Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 // Pekka Kaitaniemi, CEA and Helsinki Institute of Physics
28 // Davide Mancusi, CEA
29 // Alain Boudard, CEA
30 // Sylvie Leray, CEA
31 // Joseph Cugnon, University of Liege
32 //
33 #define INCLXX_IN_GEANT4_MODE 1
34 
35 #include "globals.hh"
36 
37 #include "G4INCLParticleTable.hh"
38 #include <algorithm>
39 // #include <cassert>
40 #include <cmath>
41 #include <cctype>
42 #include <sstream>
43 
44 #ifdef INCLXX_IN_GEANT4_MODE
45 #include "G4PhysicalConstants.hh"
46 #include "G4SystemOfUnits.hh"
47 #endif
48 
49 namespace G4INCL {
50 
52  const NaturalIsotopicDistributions *ParticleTable::theNaturalIsotopicDistributions = NULL;
53 
55  ParticleTable::NuclearMassFn ParticleTable::getTableMass;
57  ParticleTable::ParticleMassFn ParticleTable::getTableParticleMass;
59  ParticleTable::SeparationEnergyFn ParticleTable::getSeparationEnergy;
60 
61  const G4double ParticleTable::theINCLNucleonMass = 938.2796;
62  const G4double ParticleTable::theINCLPionMass = 138.0;
63  G4double ParticleTable::protonMass = 0.0;
64  G4double ParticleTable::neutronMass = 0.0;
65  G4double ParticleTable::piPlusMass = 0.0;
66  G4double ParticleTable::piMinusMass = 0.0;
67  G4double ParticleTable::piZeroMass = 0.0;
68 
69  // e^2/(4 pi epsilon_0) [MeV fm]
70  const G4double ParticleTable::eSquared = 1.439964;
71 
72  // Hard-coded values of the real particle masses (MeV/c^2)
73  G4double ParticleTable::theRealProtonMass = 938.27203;
74  G4double ParticleTable::theRealNeutronMass = 939.56536;
75  G4double ParticleTable::theRealChargedPiMass = 139.57018;
76  G4double ParticleTable::theRealPiZeroMass = 134.9766;
77 
78  const G4double ParticleTable::mediumDiffuseness[mediumNucleiTableSize] =
79  {0.0,0.0,0.0,0.0,0.0,1.78,1.77,1.77,1.77,1.71,
80  1.69,1.69,1.635,1.730,1.81,1.833,1.798,
81  1.841,0.567,0.571, 0.560,0.549,0.550,0.551,
82  0.580,0.575,0.569,0.537,0.0,0.0};
83  const G4double ParticleTable::mediumRadius[mediumNucleiTableSize] =
84  {0.0,0.0,0.0,0.0,0.0,0.334,0.327,0.479,0.631,0.838,
85  0.811,1.07,1.403,1.335,1.25,1.544,1.498,1.513,
86  2.58,2.77, 2.775,2.78,2.88,2.98,3.22,3.03,2.84,
87  3.14,0.0,0.0};
88  const G4double ParticleTable::positionRMS[clusterTableZSize][clusterTableASize] = {
89 /* A= 0 1 2 3 4 5 6 7 8 9 10 11 12 */
90 /* Z=0 */ {0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00},
91 /* Z=1 */ {0.00, 0.00, 2.10, 1.80, 1.70, 1.83, 2.60, 2.50, 0.00, 0.00, 0.00, 0.00, 0.00},
92 /* Z=2 */ {0.00, 0.00, 0.00, 1.80, 1.68, 1.70, 2.60, 2.50, 2.50, 2.50, 2.50, 0.00, 0.00},
93 /* Z=3 */ {0.00, 0.00, 0.00, 0.00, 1.70, 1.83, 2.56, 2.40, 2.50, 2.50, 2.50, 2.50, 2.50},
94 /* Z=4 */ {0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 2.60, 2.50, 2.50, 2.51, 2.50, 2.50, 2.50},
95 /* Z=5 */ {0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 2.50, 2.50, 2.50, 2.50, 2.45, 2.40, 2.50},
96 /* Z=6 */ {0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 2.50, 2.50, 2.50, 2.50, 2.47},
97 /* Z=7 */ {0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 2.50, 2.50, 2.50},
98 /* Z=8 */ {0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 2.50}
99  };
100 
101  const G4double ParticleTable::momentumRMS[clusterTableZSize][clusterTableASize] = {
102 /* A= 0 1 2 3 4 5 6 7 8 9 10 11 12 */
103 /* Z=0 */ {0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00},
104 /* Z=1 */ {0.00, 0.00, 77.0, 110., 153., 100., 100., 100., 0.00, 0.00, 0.00, 0.00, 0.00},
105 /* Z=2 */ {0.00, 0.00, 0.00, 110., 153., 100., 100., 100., 100., 100., 100., 0.00, 0.00},
106 /* Z=3 */ {0.00, 0.00, 0.00, 0.00, 153., 100., 100., 100., 100., 100., 100., 100., 100.},
107 /* Z=4 */ {0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 100., 100., 100., 100., 100., 100., 100.},
108 /* Z=5 */ {0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 100., 100., 100., 100., 100., 100., 100.},
109 /* Z=6 */ {0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 100., 100., 100., 100., 100.},
110 /* Z=7 */ {0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 100., 100., 100.},
111 /* Z=8 */ {0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 100.}
112  };
113 
115  const std::string ParticleTable::elementTable[elementTableSize] = {
116  "",
117  "H",
118  "He",
119  "Li",
120  "Be",
121  "B",
122  "C",
123  "N",
124  "O",
125  "F",
126  "Ne",
127  "Na",
128  "Mg",
129  "Al",
130  "Si",
131  "P",
132  "S",
133  "Cl",
134  "Ar",
135  "K",
136  "Ca",
137  "Sc",
138  "Ti",
139  "V",
140  "Cr",
141  "Mn",
142  "Fe",
143  "Co",
144  "Ni",
145  "Cu",
146  "Zn",
147  "Ga",
148  "Ge",
149  "As",
150  "Se",
151  "Br",
152  "Kr",
153  "Rb",
154  "Sr",
155  "Y",
156  "Zr",
157  "Nb",
158  "Mo",
159  "Tc",
160  "Ru",
161  "Rh",
162  "Pd",
163  "Ag",
164  "Cd",
165  "In",
166  "Sn",
167  "Sb",
168  "Te",
169  "I",
170  "Xe",
171  "Cs",
172  "Ba",
173  "La",
174  "Ce",
175  "Pr",
176  "Nd",
177  "Pm",
178  "Sm",
179  "Eu",
180  "Gd",
181  "Tb",
182  "Dy",
183  "Ho",
184  "Er",
185  "Tm",
186  "Yb",
187  "Lu",
188  "Hf",
189  "Ta",
190  "W",
191  "Re",
192  "Os",
193  "Ir",
194  "Pt",
195  "Au",
196  "Hg",
197  "Tl",
198  "Pb",
199  "Bi",
200  "Po",
201  "At",
202  "Rn",
203  "Fr",
204  "Ra",
205  "Ac",
206  "Th",
207  "Pa",
208  "U",
209  "Np",
210  "Pu",
211  "Am",
212  "Cm",
213  "Bk",
214  "Cf",
215  "Es",
216  "Fm",
217  "Md",
218  "No",
219  "Lr",
220  "Rf",
221  "Db",
222  "Sg",
223  "Bh",
224  "Hs",
225  "Mt",
226  "Ds",
227  "Rg",
228  "Cn"
229  };
230 
232  const std::string ParticleTable::elementIUPACDigits = "nubtqphsoe";
233 
234  /* Table for cluster decays
235  * Definition of "Stable": halflife > 1 ms
236  *
237  * These table includes decay data for clusters that INCL presently does
238  * not produce. It can't hurt.
239  *
240  * Unphysical nuclides (A<Z) are marked as stable, but should never be
241  * produced by INCL. If you find them in the output, something is fishy.
242  */
243  const ParticleTable::ClusterDecayType ParticleTable::clusterDecayMode[clusterTableZSize][clusterTableASize] =
244  {
245 /* A = 0 1 2 3 4 5 6 7 8 9 10 11 12 */
246 /* Z = 0 */ {StableCluster, StableCluster, NeutronDecay, NeutronUnbound, NeutronUnbound, NeutronUnbound, NeutronUnbound, NeutronUnbound, NeutronUnbound, NeutronUnbound, NeutronUnbound, NeutronUnbound, NeutronUnbound},
247 /* Z = 1 */ {StableCluster, StableCluster, StableCluster, StableCluster, NeutronDecay, TwoNeutronDecay, NeutronDecay, TwoNeutronDecay, NeutronUnbound, NeutronUnbound, NeutronUnbound, NeutronUnbound, NeutronUnbound},
248 /* Z = 2 */ {StableCluster, StableCluster, ProtonDecay, StableCluster, StableCluster, NeutronDecay, StableCluster, NeutronDecay, StableCluster, NeutronDecay, TwoNeutronDecay, NeutronUnbound, NeutronUnbound},
249 /* Z = 3 */ {StableCluster, StableCluster, StableCluster, ProtonUnbound, ProtonDecay, ProtonDecay, StableCluster, StableCluster, StableCluster, StableCluster, NeutronDecay, StableCluster, NeutronDecay},
250 /* Z = 4 */ {StableCluster, StableCluster, StableCluster, StableCluster, ProtonUnbound, ProtonDecay, TwoProtonDecay, StableCluster, AlphaDecay, StableCluster, StableCluster, StableCluster, StableCluster},
251 /* Z = 5 */ {StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, ProtonUnbound, TwoProtonDecay, ProtonDecay, StableCluster, ProtonDecay, StableCluster, StableCluster, StableCluster},
252 /* Z = 6 */ {StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, ProtonUnbound, ProtonUnbound, TwoProtonDecay, StableCluster, StableCluster, StableCluster, StableCluster},
253 /* Z = 7 */ {StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, ProtonUnbound, ProtonUnbound, ProtonUnbound, ProtonDecay, ProtonDecay, StableCluster},
254 /* Z = 8 */ {StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, StableCluster, ProtonUnbound, ProtonUnbound, ProtonUnbound, ProtonUnbound, ProtonDecay}
255  };
256 
260  const G4double ParticleTable::clusterPosFact[maxClusterMass+1] = {0.0, 1.0, 0.5,
261  0.33333, 0.25,
262  0.2, 0.16667,
263  0.14286, 0.125,
264  0.11111, 0.1,
265  0.09091, 0.083333};
266 
270  const G4double ParticleTable::clusterPosFact2[maxClusterMass+1] = {0.0, 1.0, 0.25,
271  0.11111, 0.0625,
272  0.04, 0.0277778,
273  0.020408, 0.015625,
274  0.012346, 0.01,
275  0.0082645, 0.0069444};
276 
277  const G4int ParticleTable::clusterZMin[maxClusterMass+1] = {0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3};
278  const G4int ParticleTable::clusterZMax[maxClusterMass+1] = {0, 0, 1, 2, 3, 3, 5, 5, 6, 6, 7, 7, 8};
279  const G4double ParticleTable::clusterPhaseSpaceCut[maxClusterMass+1] = {0.0, 70000.0, 180000.0,
280  90000.0, 90000.0,
281  128941.0 ,145607.0,
282  161365.0, 176389.0,
283  190798.0, 204681.0,
284  218109.0, 231135.0};
285  const G4double ParticleTable::effectiveNucleonMass = 938.2796;
286  const G4double ParticleTable::effectiveNucleonMass2 = 8.8036860777616e5;
287  const G4double ParticleTable::effectiveDeltaMass = 1232.0;
288  const G4double ParticleTable::effectivePionMass = 138.0;
289  const G4double ParticleTable::effectiveDeltaDecayThreshold =
290  ParticleTable::theRealNeutronMass + ParticleTable::theRealChargedPiMass
291  + 0.5;
292  const G4double ParticleTable::theINCLProtonSeparationEnergy = 6.83;
293  const G4double ParticleTable::theINCLNeutronSeparationEnergy = ParticleTable::theINCLProtonSeparationEnergy;
294  G4double ParticleTable::protonSeparationEnergy = theINCLProtonSeparationEnergy;
295  G4double ParticleTable::neutronSeparationEnergy = theINCLNeutronSeparationEnergy;
296 
297 #ifdef INCLXX_IN_GEANT4_MODE
298  G4IonTable *ParticleTable::theG4IonTable;
299 #else
300  std::vector< std::vector <G4bool> > ParticleTable::massTableMask;
301  std::vector< std::vector <G4double> > ParticleTable::massTable;
302 #endif
303 
304  void ParticleTable::initialize(Config const * const theConfig /*=0*/) {
305  protonMass = theINCLNucleonMass;
306  neutronMass = theINCLNucleonMass;
307  piPlusMass = theINCLPionMass;
308  piMinusMass = theINCLPionMass;
309  piZeroMass = theINCLPionMass;
310 
311 #ifndef INCLXX_IN_GEANT4_MODE
312  std::string dataFilePath;
313  if(theConfig)
314  dataFilePath = theConfig->getINCLXXDataFilePath();
315  readRealMasses(dataFilePath);
316 #endif
317 
318  if(theConfig && theConfig->getUseRealMasses()) {
319  getTableMass = ParticleTable::getRealMass;
320  getTableParticleMass = ParticleTable::getRealMass;
321  } else {
322  getTableMass = ParticleTable::getINCLMass;
323  getTableParticleMass = ParticleTable::getINCLMass;
324  }
325 #ifdef INCLXX_IN_GEANT4_MODE
326  G4ParticleTable *theG4ParticleTable = G4ParticleTable::GetParticleTable();
327  theG4IonTable = theG4ParticleTable->GetIonTable();
328  theRealProtonMass = theG4ParticleTable->FindParticle("proton")->GetPDGMass() / MeV;
329  theRealNeutronMass = theG4ParticleTable->FindParticle("neutron")->GetPDGMass() / MeV;
330  theRealChargedPiMass = theG4ParticleTable->FindParticle("pi+")->GetPDGMass() / MeV;
331  theRealPiZeroMass = theG4ParticleTable->FindParticle("pi0")->GetPDGMass() / MeV;
332 #endif
333 
334  // Initialise the separation-energy function
335  if(!theConfig || theConfig->getSeparationEnergyType()==INCLSeparationEnergy)
336  getSeparationEnergy = ParticleTable::getSeparationEnergyINCL;
337  else if(theConfig->getSeparationEnergyType()==RealSeparationEnergy)
338  getSeparationEnergy = ParticleTable::getSeparationEnergyReal;
339  else if(theConfig->getSeparationEnergyType()==RealForLightSeparationEnergy)
340  getSeparationEnergy = ParticleTable::getSeparationEnergyRealForLight;
341  else {
342  FATAL("Unrecognized separation-energy type in ParticleTable initialization: " << theConfig->getSeparationEnergyType() << std::endl);
343  }
344 
345  }
346 
347  G4int ParticleTable::getIsospin(const ParticleType t) {
348  // Actually this is the 3rd component of isospin (I_z) multiplied by 2!
349  if(t == Proton) {
350  return 1;
351  } else if(t == Neutron) {
352  return -1;
353  } else if(t == PiPlus) {
354  return 2;
355  } else if(t == PiMinus) {
356  return -2;
357  } else if(t == PiZero) {
358  return 0;
359  } else if(t == DeltaPlusPlus) {
360  return 3;
361  } else if(t == DeltaPlus) {
362  return 1;
363  } else if(t == DeltaZero) {
364  return -1;
365  } else if(t == DeltaMinus) {
366  return -3;
367  }
368 
369  ERROR("Requested isospin of an unknown particle!");
370  return -10; // Unknown
371  }
372 
373  std::string ParticleTable::getShortName(const ParticleSpecies s) {
374  if(s.theType==Composite)
375  return getShortName(s.theA,s.theZ);
376  else
377  return getShortName(s.theType);
378  }
379 
380  std::string ParticleTable::getName(const ParticleSpecies s) {
381  if(s.theType==Composite)
382  return getName(s.theA,s.theZ);
383  else
384  return getName(s.theType);
385  }
386 
387  std::string ParticleTable::getName(const G4int A, const G4int Z) {
388  std::stringstream stream;
389  stream << getElementName(Z) << "-" << A;
390  return stream.str();
391  }
392 
393  std::string ParticleTable::getShortName(const G4int A, const G4int Z) {
394  std::stringstream stream;
395  stream << getElementName(Z);
396  if(A>0)
397  stream << A;
398  return stream.str();
399  }
400 
401  std::string ParticleTable::getName(const ParticleType p) {
402  if(p == G4INCL::Proton) {
403  return std::string("proton");
404  } else if(p == G4INCL::Neutron) {
405  return std::string("neutron");
406  } else if(p == G4INCL::DeltaPlusPlus) {
407  return std::string("delta++");
408  } else if(p == G4INCL::DeltaPlus) {
409  return std::string("delta+");
410  } else if(p == G4INCL::DeltaZero) {
411  return std::string("delta0");
412  } else if(p == G4INCL::DeltaMinus) {
413  return std::string("delta-");
414  } else if(p == G4INCL::PiPlus) {
415  return std::string("pi+");
416  } else if(p == G4INCL::PiZero) {
417  return std::string("pi0");
418  } else if(p == G4INCL::PiMinus) {
419  return std::string("pi-");
420  } else if(p == G4INCL::Composite) {
421  return std::string("composite");
422  }
423  return std::string("unknown");
424  }
425 
426  std::string ParticleTable::getShortName(const ParticleType p) {
427  if(p == G4INCL::Proton) {
428  return std::string("p");
429  } else if(p == G4INCL::Neutron) {
430  return std::string("n");
431  } else if(p == G4INCL::DeltaPlusPlus) {
432  return std::string("d++");
433  } else if(p == G4INCL::DeltaPlus) {
434  return std::string("d+");
435  } else if(p == G4INCL::DeltaZero) {
436  return std::string("d0");
437  } else if(p == G4INCL::DeltaMinus) {
438  return std::string("d-");
439  } else if(p == G4INCL::PiPlus) {
440  return std::string("pi+");
441  } else if(p == G4INCL::PiZero) {
442  return std::string("pi0");
443  } else if(p == G4INCL::PiMinus) {
444  return std::string("pi-");
445  } else if(p == G4INCL::Composite) {
446  return std::string("comp");
447  }
448  return std::string("unknown");
449  }
450 
451  G4double ParticleTable::getINCLMass(const ParticleType pt) {
452  if(pt == Proton) {
453  return protonMass;
454  } else if(pt == Neutron) {
455  return neutronMass;
456  } else if(pt == PiPlus) {
457  return piPlusMass;
458  } else if(pt == PiMinus) {
459  return piMinusMass;
460  } else if(pt == PiZero) {
461  return piZeroMass;
462  } else {
463  ERROR("ParticleTable::getMass : Unknown particle type." << std::endl);
464  return 0.0;
465  }
466  }
467 
468  G4double ParticleTable::getRealMass(const ParticleType t) {
469  switch(t) {
470  case Proton:
471  return theRealProtonMass;
472  break;
473  case Neutron:
474  return theRealNeutronMass;
475  break;
476  case PiPlus:
477  case PiMinus:
478  return theRealChargedPiMass;
479  break;
480  case PiZero:
481  return theRealPiZeroMass;
482  break;
483  default:
484  ERROR("Particle::getRealMass : Unknown particle type." << std::endl);
485  return 0.0;
486  break;
487  }
488  }
489 
490  G4double ParticleTable::getRealMass(const G4int A, const G4int Z) {
491 // assert(A>=0);
492  // For nuclei with Z<0 or Z>A, assume that the exotic charge state is due to pions
493  if(Z<0)
494  return A*neutronMass - Z*getRealMass(PiMinus);
495  else if(Z>A)
496  return A*protonMass + (A-Z)*getRealMass(PiPlus);
497  else if(Z==0)
498  return A*getRealMass(Neutron);
499  else if(A==Z)
500  return A*getRealMass(Proton);
501  else if(A>1) {
502 #ifndef INCLXX_IN_GEANT4_MODE
503  if(ParticleTable::hasMassTable(A,Z))
504  return ParticleTable::massTable.at(Z).at(A);
505  else {
506  DEBUG("Real mass unavailable for isotope A=" << A << ", Z=" << Z
507  << ", using Weizsaecker's formula"
508  << std::endl);
509  return getWeizsaeckerMass(A,Z);
510  }
511 #else
512  return theG4IonTable->GetNucleusMass(Z,A) / MeV;
513 #endif
514  } else
515  return 0.;
516  }
517 
518  G4double ParticleTable::getINCLMass(const G4int A, const G4int Z) {
519 // assert(A>=0);
520  // For nuclei with Z<0 or Z>A, assume that the exotic charge state is due to pions
521  if(Z<0)
522  return A*neutronMass - Z*getINCLMass(PiMinus);
523  else if(Z>A)
524  return A*protonMass + (A-Z)*getINCLMass(PiPlus);
525  else if(A>1)
526  return Z*(protonMass - protonSeparationEnergy) + (A-Z)*(neutronMass - neutronSeparationEnergy);
527  else if(A==1 && Z==0)
528  return getINCLMass(Neutron);
529  else if(A==1 && Z==1)
530  return getINCLMass(Proton);
531  else
532  return 0.;
533  }
534 
535  G4double ParticleTable::getNuclearRadius(const G4int A, const G4int Z) {
536 // assert(A>0 && Z>=0 && Z<=A);
537  if(A >= 19 || (A < 6 && A >= 2)) {
538  // For large (Woods-Saxon or Modified Harmonic Oscillator) or small
539  // (Gaussian) nuclei, the radius parameter is just the nuclear radius
540  return getRadiusParameter(A,Z);
541  } else if(A < clusterTableASize && Z < clusterTableZSize && A >= 6) {
542  const G4double thisRMS = positionRMS[Z][A];
543  if(thisRMS>0.0)
544  return thisRMS;
545  else {
546  ERROR("ParticleTable::getRadiusParameter : Radius for nucleus A = " << A << " Z = " << Z << " is ==0.0" << std::endl);
547  return 0.0;
548  }
549  } else if(A < 19) {
550  const G4double theRadiusParameter = getRadiusParameter(A, Z);
551  const G4double theDiffusenessParameter = getSurfaceDiffuseness(A, Z);
552  // The formula yields the nuclear RMS radius based on the parameters of
553  // the nuclear-density function
554  return 1.581*theDiffusenessParameter*
555  (2.+5.*theRadiusParameter)/(2.+3.*theRadiusParameter);
556  } else {
557  ERROR("ParticleTable::getNuclearRadius: No radius for nucleus A = " << A << " Z = " << Z << std::endl);
558  return 0.0;
559  }
560  }
561 
562  G4double ParticleTable::getRadiusParameter(const G4int A, const G4int Z) {
563 // assert(A>0 && Z>=0 && Z<=A);
564  if(A >= 28) {
565  // phenomenological radius fit
566  return (2.745e-4 * A + 1.063) * std::pow(A, 1.0/3.0);
567  } else if(A < 6 && A >= 2) {
568  if(Z<clusterTableZSize) {
569  const G4double thisRMS = positionRMS[Z][A];
570  if(thisRMS>0.0)
571  return thisRMS;
572  else {
573  ERROR("ParticleTable::getRadiusParameter : Radius for nucleus A = " << A << " Z = " << Z << " is ==0.0" << std::endl);
574  return 0.0;
575  }
576  } else {
577  ERROR("ParticleTable::getRadiusParameter : No radius for nucleus A = " << A << " Z = " << Z << std::endl);
578  return 0.0;
579  }
580  } else if(A < 28 && A >= 6) {
581  return mediumRadius[A-1];
582  // return 1.581*mediumDiffuseness[A-1]*(2.+5.*mediumRadius[A-1])/(2.+3.*mediumRadius[A-1]);
583  } else {
584  ERROR("ParticleTable::getRadiusParameter: No radius for nucleus A = " << A << " Z = " << Z << std::endl);
585  return 0.0;
586  }
587  }
588 
589  G4double ParticleTable::getMaximumNuclearRadius(const G4int A, const G4int Z) {
590  const G4double XFOISA = 8.0;
591  if(A >= 19) {
592  return getNuclearRadius(A,Z) + XFOISA * getSurfaceDiffuseness(A,Z);
593  } else if(A < 19 && A >= 6) {
594  return 5.5 + 0.3 * (G4double(A) - 6.0)/12.0;
595  } else if(A >= 2) {
596  return ParticleTable::getNuclearRadius(A, Z) + 4.5;
597  } else {
598  ERROR("ParticleTable::getMaximumNuclearRadius : No maximum radius for nucleus A = " << A << " Z = " << Z << std::endl);
599  return 0.0;
600  }
601  }
602 
603 #ifdef INCLXX_IN_GEANT4_MODE
604  G4double ParticleTable::getSurfaceDiffuseness(const G4int A, const G4int /*Z*/ ) {
605 #else
606  G4double ParticleTable::getSurfaceDiffuseness(const G4int A, const G4int Z) {
607 #endif // INCLXX_IN_GEANT4_MODE
608 
609  if(A >= 28) {
610  return 1.63e-4 * A + 0.510;
611  } else if(A < 28 && A >= 19) {
612  return mediumDiffuseness[A-1];
613  } else if(A < 19 && A >= 6) {
614  return mediumDiffuseness[A-1];
615  } else if(A < 6 && A >= 2) {
616  ERROR("ParticleTable::getSurfaceDiffuseness: was called for A = " << A << " Z = " << Z << std::endl);
617  return 0.0;
618  } else {
619  ERROR("ParticleTable::getSurfaceDiffuseness: No diffuseness for nucleus A = " << A << " Z = " << Z << std::endl);
620  return 0.0;
621  }
622  }
623 
624  std::string ParticleTable::getElementName(const G4int Z) {
625  if(Z<1) {
626  WARN("getElementName called with Z<1" << std::endl);
627  return elementTable[0];
628  } else if(Z<elementTableSize)
629  return elementTable[Z];
630  else
631  return getIUPACElementName(Z);
632  }
633 
634 #ifndef INCLXX_IN_GEANT4_MODE
635  void ParticleTable::readRealMasses(std::string const &path) {
636  // Clear the existing tables, if any
637  massTableMask.clear();
638  massTable.clear();
639 
640  // File name
641  std::string fileName(path + "/walletlifetime.dat");
642  DEBUG("Reading real nuclear masses from file " << fileName << std::endl);
643 
644  // Open the file stream
645  std::ifstream massTableIn(fileName.c_str());
646  if(!massTableIn.good()) {
647  FATAL("Cannot open " << fileName << " data file." << std::endl);
648  return;
649  }
650 
651  // Read in the data
652  unsigned int Z, A;
653  const G4double amu = 931.494061; // atomic mass unit in MeV/c^2
654  const G4double eMass = 0.5109988; // electron mass in MeV/c^2
655  G4double excess;
656  massTableIn >> A >> Z >> excess;
657  do {
658  // Dynamically determine the table size
659  if(Z>=massTable.size()) {
660  massTable.resize(Z+1);
661  massTableMask.resize(Z+1);
662  }
663  if(A>=massTable[Z].size()) {
664  massTable[Z].resize(A+1, 0.);
665  massTableMask[Z].resize(A+1, false);
666  }
667 
668  massTable.at(Z).at(A) = A*amu + excess - Z*eMass;
669  massTableMask.at(Z).at(A) = true;
670 
671  massTableIn >> A >> Z >> excess;
672  } while(massTableIn.good());
673 
674  }
675 #endif
676 
677  std::string ParticleTable::getIUPACElementName(const G4int Z) {
678  std::stringstream elementStream;
679  elementStream << Z;
680  std::string elementName = elementStream.str();
681  std::transform(elementName.begin(), elementName.end(), elementName.begin(), ParticleTable::intToIUPAC);
682  elementName[0] = std::toupper(elementName.at(0));
683  return elementName;
684  }
685 
686  G4int ParticleTable::parseIUPACElement(std::string const &s) {
687  // Normalise to lower case
688  std::string elementName(s);
689  std::transform(elementName.begin(), elementName.end(), elementName.begin(), ::tolower);
690  // Return 0 if the element name contains anything but IUPAC digits
691  if(elementName.find_first_not_of(elementIUPACDigits)!=std::string::npos)
692  return 0;
693  std::transform(elementName.begin(), elementName.end(), elementName.begin(), ParticleTable::iupacToInt);
694  std::stringstream elementStream(elementName);
695  G4int Z;
696  elementStream >> Z;
697  return Z;
698  }
699 
700 }
Generated on Sat May 25 2013 14:34:00 for Geant4 by   doxygen 1.8.4