Geant4  10.03.p03
 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 // 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"
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 theINCLEtaMass = 547.862;
66  const G4double theINCLOmegaMass = 782.65;
67  const G4double theINCLEtaPrimeMass = 957.78;
68  const G4double theINCLPhotonMass = 0.0;
69  G4ThreadLocal G4double protonMass = 0.0;
70  G4ThreadLocal G4double neutronMass = 0.0;
71  G4ThreadLocal G4double piPlusMass = 0.0;
72  G4ThreadLocal G4double piMinusMass = 0.0;
73  G4ThreadLocal G4double piZeroMass = 0.0;
74  G4ThreadLocal G4double etaMass = 0.0;
75  G4ThreadLocal G4double omegaMass = 0.0;
76  G4ThreadLocal G4double etaPrimeMass = 0.0;
77  G4ThreadLocal G4double photonMass = 0.0;
78 
79  // Hard-coded values of the real particle masses (MeV/c^2)
80  G4ThreadLocal G4double theRealProtonMass = 938.27203;
81  G4ThreadLocal G4double theRealNeutronMass = 939.56536;
82  G4ThreadLocal G4double theRealChargedPiMass = 139.57018;
83  G4ThreadLocal G4double theRealPiZeroMass = 134.9766;
84  G4ThreadLocal G4double theRealEtaMass = 547.862;
85  G4ThreadLocal G4double theRealOmegaMass = 782.65;
86  G4ThreadLocal G4double theRealEtaPrimeMass = 957.78;
87  G4ThreadLocal G4double theRealPhotonMass = 0.0;
88 
89  // Width (second)
90  const G4double theChargedPiWidth = 2.6033e-08;
91  const G4double thePiZeroWidth = 8.52e-17;
92  const G4double theEtaWidth = 5.025e-19; // 1.31 keV
93  const G4double theOmegaWidth = 7.7528e-23; // 8.49 MeV
94  const G4double theEtaPrimeWidth = 3.3243e-21; // 0.198 MeV
95  G4ThreadLocal G4double piPlusWidth = 0.0;
96  G4ThreadLocal G4double piMinusWidth = 0.0;
97  G4ThreadLocal G4double piZeroWidth = 0.0;
98  G4ThreadLocal G4double etaWidth = 0.0;
99  G4ThreadLocal G4double omegaWidth = 0.0;
100  G4ThreadLocal G4double etaPrimeWidth = 0.0;
101 
102 
103  const G4int mediumNucleiTableSize = 30;
104 
105  const G4double mediumDiffuseness[mediumNucleiTableSize] =
106  {0.0,0.0,0.0,0.0,0.0,1.78,1.77,1.77,1.77,1.71,
107  1.69,1.69,1.635,1.730,1.81,1.833,1.798,
108  1.841,0.567,0.571, 0.560,0.549,0.550,0.551,
109  0.580,0.575,0.569,0.537,0.0,0.0};
110  const G4double mediumRadius[mediumNucleiTableSize] =
111  {0.0,0.0,0.0,0.0,0.0,0.334,0.327,0.479,0.631,0.838,
112  0.811,1.07,1.403,1.335,1.25,1.544,1.498,1.513,
113  2.58,2.77, 2.775,2.78,2.88,2.98,3.22,3.03,2.84,
114  3.14,0.0,0.0};
115 
116  const G4double positionRMS[clusterTableZSize][clusterTableASize] = {
117  /* A= 0 1 2 3 4 5 6 7 8 9 10 11 12 */
118  /* 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},
119  /* 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},
120  /* 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},
121  /* 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},
122  /* 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},
123  /* 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},
124  /* 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},
125  /* 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},
126  /* 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}
127  };
128 
129  const G4double momentumRMS[clusterTableZSize][clusterTableASize] = {
130  /* A= 0 1 2 3 4 5 6 7 8 9 10 11 12 */
131  /* 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},
132  /* Z=1 */ {-1.0, -1.0, 77.0, 110., 153., 100., 100., 100., -1.0, -1.0, -1.0, -1.0, -1.0},
133  /* Z=2 */ {-1.0, -1.0, -1.0, 110., 153., 100., 100., 100., 100., 100., 100., -1.0, -1.0},
134  /* Z=3 */ {-1.0, -1.0, -1.0, -1.0, 153., 100., 100., 100., 100., 100., 100., 100., 100.},
135  /* Z=4 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 100., 100., 100., 100., 100., 100., 100.},
136  /* Z=5 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 100., 100., 100., 100., 100., 100., 100.},
137  /* Z=6 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 100., 100., 100., 100., 100.},
138  /* 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.},
139  /* 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.}
140  };
141 
142  const G4int elementTableSize = 113; // up to Cn
143 
145  const std::string elementTable[elementTableSize] = {
146  "",
147  "H",
148  "He",
149  "Li",
150  "Be",
151  "B",
152  "C",
153  "N",
154  "O",
155  "F",
156  "Ne",
157  "Na",
158  "Mg",
159  "Al",
160  "Si",
161  "P",
162  "S",
163  "Cl",
164  "Ar",
165  "K",
166  "Ca",
167  "Sc",
168  "Ti",
169  "V",
170  "Cr",
171  "Mn",
172  "Fe",
173  "Co",
174  "Ni",
175  "Cu",
176  "Zn",
177  "Ga",
178  "Ge",
179  "As",
180  "Se",
181  "Br",
182  "Kr",
183  "Rb",
184  "Sr",
185  "Y",
186  "Zr",
187  "Nb",
188  "Mo",
189  "Tc",
190  "Ru",
191  "Rh",
192  "Pd",
193  "Ag",
194  "Cd",
195  "In",
196  "Sn",
197  "Sb",
198  "Te",
199  "I",
200  "Xe",
201  "Cs",
202  "Ba",
203  "La",
204  "Ce",
205  "Pr",
206  "Nd",
207  "Pm",
208  "Sm",
209  "Eu",
210  "Gd",
211  "Tb",
212  "Dy",
213  "Ho",
214  "Er",
215  "Tm",
216  "Yb",
217  "Lu",
218  "Hf",
219  "Ta",
220  "W",
221  "Re",
222  "Os",
223  "Ir",
224  "Pt",
225  "Au",
226  "Hg",
227  "Tl",
228  "Pb",
229  "Bi",
230  "Po",
231  "At",
232  "Rn",
233  "Fr",
234  "Ra",
235  "Ac",
236  "Th",
237  "Pa",
238  "U",
239  "Np",
240  "Pu",
241  "Am",
242  "Cm",
243  "Bk",
244  "Cf",
245  "Es",
246  "Fm",
247  "Md",
248  "No",
249  "Lr",
250  "Rf",
251  "Db",
252  "Sg",
253  "Bh",
254  "Hs",
255  "Mt",
256  "Ds",
257  "Rg",
258  "Cn"
259  };
260 
262  const std::string elementIUPACDigits = "nubtqphsoe";
263 
264 #define INCL_DEFAULT_SEPARATION_ENERGY 6.83
265  const G4double theINCLProtonSeparationEnergy = INCL_DEFAULT_SEPARATION_ENERGY;
266  const G4double theINCLNeutronSeparationEnergy = INCL_DEFAULT_SEPARATION_ENERGY;
267  G4ThreadLocal G4double protonSeparationEnergy = INCL_DEFAULT_SEPARATION_ENERGY;
268  G4ThreadLocal G4double neutronSeparationEnergy = INCL_DEFAULT_SEPARATION_ENERGY;
269 #undef INCL_DEFAULT_SEPARATION_ENERGY
270 
271  G4ThreadLocal G4double rpCorrelationCoefficient[UnknownParticle];
272 
273  G4ThreadLocal G4double neutronSkin = 0.0;
274  G4ThreadLocal G4double neutronHalo = 0.0;
275 
276 #ifdef INCLXX_IN_GEANT4_MODE
277  G4ThreadLocal G4IonTable *theG4IonTable;
278 #endif
279 
281  G4ThreadLocal G4double constantFermiMomentum = 0.0;
282 
284  char iupacToInt(char c) {
285  return (char)(((G4int)'0')+elementIUPACDigits.find(c));
286  }
287 
289  char intToIUPAC(char n) { return elementIUPACDigits.at(n); }
290 
292  const NaturalIsotopicDistributions *getNaturalIsotopicDistributions() {
293  if(!theNaturalIsotopicDistributions)
294  theNaturalIsotopicDistributions = new NaturalIsotopicDistributions;
295  return theNaturalIsotopicDistributions;
296  }
297 
298  } // namespace
299 
300  void initialize(Config const * const theConfig /*=0*/) {
301  protonMass = theINCLNucleonMass;
302  neutronMass = theINCLNucleonMass;
303  piPlusMass = theINCLPionMass;
304  piMinusMass = theINCLPionMass;
305  piZeroMass = theINCLPionMass;
306  etaMass = theINCLEtaMass;
307  omegaMass = theINCLOmegaMass;
308  etaPrimeMass = theINCLEtaPrimeMass;
309  photonMass = theINCLPhotonMass;
310 
311  if(theConfig && theConfig->getUseRealMasses()) {
314  } else {
317  }
318 
319 #ifndef INCLXX_IN_GEANT4_MODE
320  std::string dataFilePath;
321  if(theConfig)
322  dataFilePath = theConfig->getINCLXXDataFilePath();
324 #endif
325 
326 #ifdef INCLXX_IN_GEANT4_MODE
327  G4ParticleTable *theG4ParticleTable = G4ParticleTable::GetParticleTable();
328  theG4IonTable = theG4ParticleTable->GetIonTable();
329  theRealProtonMass = theG4ParticleTable->FindParticle("proton")->GetPDGMass() / MeV;
330  theRealNeutronMass = theG4ParticleTable->FindParticle("neutron")->GetPDGMass() / MeV;
331  theRealChargedPiMass = theG4ParticleTable->FindParticle("pi+")->GetPDGMass() / MeV;
332  theRealPiZeroMass = theG4ParticleTable->FindParticle("pi0")->GetPDGMass() / MeV;
333  theRealEtaMass = theG4ParticleTable->FindParticle("eta")->GetPDGMass() / MeV;
334  theRealOmegaMass = theG4ParticleTable->FindParticle("omega")->GetPDGMass() / MeV;
335  theRealEtaPrimeMass = theG4ParticleTable->FindParticle("eta_prime")->GetPDGMass() / MeV;
336  theRealPhotonMass = theG4ParticleTable->FindParticle("gamma")->GetPDGMass() / MeV;
337 #endif
338 
339  minDeltaMass = theRealNeutronMass + theRealChargedPiMass + 0.5;
341  minDeltaMassRndm = std::atan((minDeltaMass-effectiveDeltaMass)*2./effectiveDeltaWidth);
342 
343  piPlusWidth = theChargedPiWidth;
344  piMinusWidth = theChargedPiWidth;
345  piZeroWidth = thePiZeroWidth;
346  etaWidth = theEtaWidth;
347  omegaWidth = theOmegaWidth;
348  etaPrimeWidth = theEtaPrimeWidth;
349 
350 
351  // Initialise the separation-energy function
352  if(!theConfig || theConfig->getSeparationEnergyType()==INCLSeparationEnergy)
354  else if(theConfig->getSeparationEnergyType()==RealSeparationEnergy)
358  else {
359  INCL_FATAL("Unrecognized separation-energy type in ParticleTable initialization: " << theConfig->getSeparationEnergyType() << '\n');
360  return;
361  }
362 
363  // Initialise the Fermi-momentum function
364  if(!theConfig || theConfig->getFermiMomentumType()==ConstantFermiMomentum) {
366  if(theConfig) {
367  const G4double aFermiMomentum = theConfig->getFermiMomentum();
368  if(aFermiMomentum>0.)
369  constantFermiMomentum = aFermiMomentum;
370  else
371  constantFermiMomentum = PhysicalConstants::Pf;
372  } else {
373  constantFermiMomentum = PhysicalConstants::Pf;
374  }
375  } else if(theConfig->getFermiMomentumType()==ConstantLightFermiMomentum)
377  else if(theConfig->getFermiMomentumType()==MassDependentFermiMomentum)
379  else {
380  INCL_FATAL("Unrecognized Fermi-momentum type in ParticleTable initialization: " << theConfig->getFermiMomentumType() << '\n');
381  return;
382  }
383 
384  // Initialise the r-p correlation coefficients
385  std::fill(rpCorrelationCoefficient, rpCorrelationCoefficient + UnknownParticle, 1.);
386  if(theConfig) {
387  rpCorrelationCoefficient[Proton] = theConfig->getRPCorrelationCoefficient(Proton);
388  rpCorrelationCoefficient[Neutron] = theConfig->getRPCorrelationCoefficient(Neutron);
389  }
390 
391  // Initialise the neutron-skin parameters
392  if(theConfig) {
393  neutronSkin = theConfig->getNeutronSkin();
394  neutronHalo = theConfig->getNeutronHalo();
395  }
396 
397  }
398 
400  // Actually this is the 3rd component of isospin (I_z) multiplied by 2!
401  if(t == Proton) {
402  return 1;
403  } else if(t == Neutron) {
404  return -1;
405  } else if(t == PiPlus) {
406  return 2;
407  } else if(t == PiMinus) {
408  return -2;
409  } else if(t == PiZero) {
410  return 0;
411  } else if(t == DeltaPlusPlus) {
412  return 3;
413  } else if(t == DeltaPlus) {
414  return 1;
415  } else if(t == DeltaZero) {
416  return -1;
417  } else if(t == DeltaMinus) {
418  return -3;
419  } else if(t == Eta) {
420  return 0;
421  } else if(t == Omega) {
422  return 0;
423  } else if(t == EtaPrime) {
424  return 0;
425  } else if(t == Photon) {
426  return 0;
427  }
428 
429  INCL_ERROR("Requested isospin of an unknown particle!");
430  return -10; // Unknown
431  }
432 
433  std::string getShortName(const ParticleSpecies &s) {
434  if(s.theType==Composite)
435  return getShortName(s.theA,s.theZ);
436  else
437  return getShortName(s.theType);
438  }
439 
440  std::string getName(const ParticleSpecies &s) {
441  if(s.theType==Composite)
442  return getName(s.theA,s.theZ);
443  else
444  return getName(s.theType);
445  }
446 
447  std::string getName(const G4int A, const G4int Z) {
448  std::stringstream stream;
449  stream << getElementName(Z) << "-" << A;
450  return stream.str();
451  }
452 
453  std::string getShortName(const G4int A, const G4int Z) {
454  std::stringstream stream;
455  stream << getElementName(Z);
456  if(A>0)
457  stream << A;
458  return stream.str();
459  }
460 
461  std::string getName(const ParticleType p) {
462  if(p == G4INCL::Proton) {
463  return std::string("proton");
464  } else if(p == G4INCL::Neutron) {
465  return std::string("neutron");
466  } else if(p == G4INCL::DeltaPlusPlus) {
467  return std::string("delta++");
468  } else if(p == G4INCL::DeltaPlus) {
469  return std::string("delta+");
470  } else if(p == G4INCL::DeltaZero) {
471  return std::string("delta0");
472  } else if(p == G4INCL::DeltaMinus) {
473  return std::string("delta-");
474  } else if(p == G4INCL::PiPlus) {
475  return std::string("pi+");
476  } else if(p == G4INCL::PiZero) {
477  return std::string("pi0");
478  } else if(p == G4INCL::PiMinus) {
479  return std::string("pi-");
480  } else if(p == G4INCL::Composite) {
481  return std::string("composite");
482  } else if(p == G4INCL::Eta) {
483  return std::string("eta");
484  } else if(p == G4INCL::Omega) {
485  return std::string("omega");
486  } else if(p == G4INCL::EtaPrime) {
487  return std::string("etaprime");
488  } else if(p == G4INCL::Photon) {
489  return std::string("photon");
490  }
491  return std::string("unknown");
492  }
493 
494  std::string getShortName(const ParticleType p) {
495  if(p == G4INCL::Proton) {
496  return std::string("p");
497  } else if(p == G4INCL::Neutron) {
498  return std::string("n");
499  } else if(p == G4INCL::DeltaPlusPlus) {
500  return std::string("d++");
501  } else if(p == G4INCL::DeltaPlus) {
502  return std::string("d+");
503  } else if(p == G4INCL::DeltaZero) {
504  return std::string("d0");
505  } else if(p == G4INCL::DeltaMinus) {
506  return std::string("d-");
507  } else if(p == G4INCL::PiPlus) {
508  return std::string("pi+");
509  } else if(p == G4INCL::PiZero) {
510  return std::string("pi0");
511  } else if(p == G4INCL::PiMinus) {
512  return std::string("pi-");
513  } else if(p == G4INCL::Composite) {
514  return std::string("comp");
515  } else if(p == G4INCL::Eta) {
516  return std::string("eta");
517  } else if(p == G4INCL::Omega) {
518  return std::string("omega");
519  } else if(p == G4INCL::EtaPrime) {
520  return std::string("etap");
521  } else if(p == G4INCL::Photon) {
522  return std::string("photon");
523  }
524  return std::string("unknown");
525  }
526 
528  if(pt == Proton) {
529  return protonMass;
530  } else if(pt == Neutron) {
531  return neutronMass;
532  } else if(pt == PiPlus) {
533  return piPlusMass;
534  } else if(pt == PiMinus) {
535  return piMinusMass;
536  } else if(pt == PiZero) {
537  return piZeroMass;
538  } else if(pt == Eta) {
539  return etaMass;
540  } else if(pt == Omega) {
541  return omegaMass;
542  } else if(pt == EtaPrime) {
543  return etaPrimeMass;
544  } else if(pt == Photon) {
545  return photonMass;
546  } else {
547  INCL_ERROR("getMass : Unknown particle type." << '\n');
548  return 0.0;
549  }
550  }
551 
553  switch(t) {
554  case Proton:
555  return theRealProtonMass;
556  break;
557  case Neutron:
558  return theRealNeutronMass;
559  break;
560  case PiPlus:
561  case PiMinus:
562  return theRealChargedPiMass;
563  break;
564  case PiZero:
565  return theRealPiZeroMass;
566  break;
567  case Eta:
568  return theRealEtaMass;
569  break;
570  case Omega:
571  return theRealOmegaMass;
572  break;
573  case EtaPrime:
574  return theRealEtaPrimeMass;
575  break;
576  case Photon:
577  return theRealPhotonMass;
578  break;
579  default:
580  INCL_ERROR("Particle::getRealMass : Unknown particle type." << '\n');
581  return 0.0;
582  break;
583  }
584  }
585 
586  G4double getRealMass(const G4int A, const G4int Z) {
587 // assert(A>=0);
588  // For nuclei with Z<0 or Z>A, assume that the exotic charge state is due to pions
589  if(Z<0)
590  return A*neutronMass - Z*getRealMass(PiMinus);
591  else if(Z>A)
592  return A*protonMass + (A-Z)*getRealMass(PiPlus);
593  else if(Z==0)
594  return A*getRealMass(Neutron);
595  else if(A==Z)
596  return A*getRealMass(Proton);
597  else if(A>1) {
598 #ifndef INCLXX_IN_GEANT4_MODE
599  return ::G4INCL::NuclearMassTable::getMass(A,Z);
600 #else
601  return theG4IonTable->GetNucleusMass(Z,A) / MeV;
602 #endif
603  } else
604  return 0.;
605  }
606 
607  G4double getINCLMass(const G4int A, const G4int Z) {
608 // assert(A>=0);
609  // For nuclei with Z<0 or Z>A, assume that the exotic charge state is due to pions
610  if(Z<0)
611  return A*neutronMass - Z*getINCLMass(PiMinus);
612  else if(Z>A)
613  return A*protonMass + (A-Z)*getINCLMass(PiPlus);
614  else if(A>1)
615  return Z*(protonMass - protonSeparationEnergy) + (A-Z)*(neutronMass - neutronSeparationEnergy);
616  else if(A==1 && Z==0)
617  return getINCLMass(Neutron);
618  else if(A==1 && Z==1)
619  return getINCLMass(Proton);
620  else
621  return 0.;
622  }
623 
624  G4double getTableQValue(const G4int A1, const G4int Z1, const G4int A2, const G4int Z2) {
625  return getTableMass(A1,Z1) + getTableMass(A2,Z2) - getTableMass(A1+A2,Z1+Z2);
626  }
627 
628  G4double getTableQValue(const G4int A1, const G4int Z1, const G4int A2, const G4int Z2, const G4int A3, const G4int Z3) {
629  return getTableMass(A1,Z1) + getTableMass(A2,Z2) - getTableMass(A3,Z3) - getTableMass(A1+A2-A3,Z1+Z2-Z3);
630  }
631 
633  if(p.theType == Composite)
634  return (*getTableMass)(p.theA, p.theZ);
635  else
636  return (*getTableParticleMass)(p.theType);
637  }
638 
640  switch(t) {
641  case Proton:
642  case Neutron:
643  case DeltaPlusPlus:
644  case DeltaPlus:
645  case DeltaZero:
646  case DeltaMinus:
647  return 1;
648  break;
649  case PiPlus:
650  case PiMinus:
651  case PiZero:
652  case Eta:
653  case Omega:
654  case EtaPrime:
655  case Photon:
656  return 0;
657  break;
658  default:
659  return 0;
660  break;
661  }
662  }
663 
665  switch(t) {
666  case DeltaPlusPlus:
667  return 2;
668  break;
669  case Proton:
670  case DeltaPlus:
671  case PiPlus:
672  return 1;
673  break;
674  case Neutron:
675  case DeltaZero:
676  case PiZero:
677  case Eta:
678  case Omega:
679  case EtaPrime:
680  case Photon:
681  return 0;
682  break;
683  case DeltaMinus:
684  case PiMinus:
685  return -1;
686  break;
687  default:
688  return 0;
689  break;
690  }
691  }
692 
694 // assert(A>=0);
695  if(A >= 19 || (A < 6 && A >= 2)) {
696  // For large (Woods-Saxon or Modified Harmonic Oscillator) or small
697  // (Gaussian) nuclei, the radius parameter is just the nuclear radius
698  return getRadiusParameter(t,A,Z);
699  } else if(A < clusterTableASize && Z>=0 && Z < clusterTableZSize && A >= 6) {
700  const G4double thisRMS = positionRMS[Z][A];
701  if(thisRMS>0.0)
702  return thisRMS;
703  else {
704  INCL_DEBUG("getNuclearRadius: Radius for nucleus A = " << A << " Z = " << Z << " is not available" << '\n'
705  << "returning radius for C12");
706  return positionRMS[6][12];
707  }
708  } else if(A < 19) {
709  const G4double theRadiusParameter = getRadiusParameter(t, A, Z);
710  const G4double theDiffusenessParameter = getSurfaceDiffuseness(t, A, Z);
711  // The formula yields the nuclear RMS radius based on the parameters of
712  // the nuclear-density function
713  return 1.225*theDiffusenessParameter*
714  std::sqrt((2.+5.*theRadiusParameter)/(2.+3.*theRadiusParameter));
715  } else {
716  INCL_ERROR("getNuclearRadius: No radius for nucleus A = " << A << " Z = " << Z << '\n');
717  return 0.0;
718  }
719  }
720 
723  }
724 
726 // assert(A>0);
727  if(A >= 28) {
728  // phenomenological radius fit
729  G4double r0 = (2.745e-4 * A + 1.063) * std::pow(A, 1.0/3.0);
730  if(t==Neutron)
731  r0 += neutronSkin;
732  return r0;
733  } else if(A < 6 && A >= 2) {
734  if(Z<clusterTableZSize && Z>=0) {
735  const G4double thisRMS = positionRMS[Z][A];
736  if(thisRMS>0.0)
737  return thisRMS;
738  else {
739  INCL_DEBUG("getRadiusParameter: Radius for nucleus A = " << A << " Z = " << Z << " is not available" << '\n'
740  << "returning radius for C12");
741  return positionRMS[6][12];
742  }
743  } else {
744  INCL_DEBUG("getRadiusParameter: Radius for nucleus A = " << A << " Z = " << Z << " is not available" << '\n'
745  << "returning radius for C12");
746  return positionRMS[6][12];
747  }
748  } else if(A < 28 && A >= 6) {
749  return mediumRadius[A-1];
750  // return 1.581*mediumDiffuseness[A-1]*(2.+5.*mediumRadius[A-1])/(2.+3.*mediumRadius[A-1]);
751  } else {
752  INCL_ERROR("getRadiusParameter: No radius for nucleus A = " << A << " Z = " << Z << '\n');
753  return 0.0;
754  }
755  }
756 
758  const G4double XFOISA = 8.0;
759  if(A >= 19) {
760  return getNuclearRadius(t,A,Z) + XFOISA * getSurfaceDiffuseness(t,A,Z);
761  } else if(A < 19 && A >= 6) {
762  return 5.5 + 0.3 * (G4double(A) - 6.0)/12.0;
763  } else if(A >= 2) {
764  return getNuclearRadius(t, A, Z) + 4.5;
765  } else {
766  INCL_ERROR("getMaximumNuclearRadius : No maximum radius for nucleus A = " << A << " Z = " << Z << '\n');
767  return 0.0;
768  }
769  }
770 
772  if(A >= 28) {
773  G4double a = 1.63e-4 * A + 0.510;
774  if(t==Neutron)
775  a += neutronHalo;
776  return a;
777  } else if(A < 28 && A >= 19) {
778  return mediumDiffuseness[A-1];
779  } else if(A < 19 && A >= 6) {
780  return mediumDiffuseness[A-1];
781  } else if(A < 6 && A >= 2) {
782  INCL_ERROR("getSurfaceDiffuseness: was called for A = " << A << " Z = " << Z << '\n');
783  return 0.0;
784  } else {
785  INCL_ERROR("getSurfaceDiffuseness: No diffuseness for nucleus A = " << A << " Z = " << Z << '\n');
786  return 0.0;
787  }
788  }
789 
791 // assert(Z>=0 && A>=0 && Z<=A);
793  }
794 
795  G4double getSeparationEnergyINCL(const ParticleType t, const G4int /*A*/, const G4int /*Z*/) {
796  if(t==Proton)
797  return theINCLProtonSeparationEnergy;
798  else if(t==Neutron)
799  return theINCLNeutronSeparationEnergy;
800  else {
801  INCL_ERROR("ParticleTable::getSeparationEnergyINCL : Unknown particle type." << '\n');
802  return 0.0;
803  }
804  }
805 
807  // Real separation energies for all nuclei
808  if(t==Proton)
809  return (*getTableParticleMass)(Proton) + (*getTableMass)(A-1,Z-1) - (*getTableMass)(A,Z);
810  else if(t==Neutron)
811  return (*getTableParticleMass)(Neutron) + (*getTableMass)(A-1,Z) - (*getTableMass)(A,Z);
812  else {
813  INCL_ERROR("ParticleTable::getSeparationEnergyReal : Unknown particle type." << '\n');
814  return 0.0;
815  }
816  }
817 
819  // Real separation energies for light nuclei, fixed values for heavy nuclei
821  return getSeparationEnergyReal(t, A, Z);
822  else
823  return getSeparationEnergyINCL(t, A, Z);
824  }
825 
826  G4double getProtonSeparationEnergy() { return protonSeparationEnergy; }
827 
828  G4double getNeutronSeparationEnergy() { return neutronSeparationEnergy; }
829 
830  void setProtonSeparationEnergy(const G4double s) { protonSeparationEnergy = s; }
831 
832  void setNeutronSeparationEnergy(const G4double s) { neutronSeparationEnergy = s; }
833 
834  std::string getElementName(const G4int Z) {
835  if(Z<1) {
836  INCL_WARN("getElementName called with Z<1" << '\n');
837  return elementTable[0];
838  } else if(Z<elementTableSize)
839  return elementTable[Z];
840  else
841  return getIUPACElementName(Z);
842  }
843 
844  std::string getIUPACElementName(const G4int Z) {
845  std::stringstream elementStream;
846  elementStream << Z;
847  std::string elementName = elementStream.str();
848  std::transform(elementName.begin(), elementName.end(), elementName.begin(), intToIUPAC);
849  elementName[0] = std::toupper(elementName.at(0));
850  return elementName;
851  }
852 
853  G4int parseElement(std::string pS) {
854  // Normalize the element name
855  std::transform(pS.begin(), pS.end(), pS.begin(), ::tolower);
856  pS[0] = ::toupper(pS[0]);
857 
858  const std::string *iter = std::find(elementTable, elementTable+elementTableSize, pS);
859  if(iter != elementTable+elementTableSize)
860  return iter - elementTable;
861  else
863  }
864 
865  G4int parseIUPACElement(std::string const &s) {
866  // Normalise to lower case
867  std::string elementName(s);
868  std::transform(elementName.begin(), elementName.end(), elementName.begin(), ::tolower);
869  // Return 0 if the element name contains anything but IUPAC digits
870  if(elementName.find_first_not_of(elementIUPACDigits)!=std::string::npos)
871  return 0;
872  std::transform(elementName.begin(), elementName.end(), elementName.begin(), iupacToInt);
873  std::stringstream elementStream(elementName);
874  G4int Z;
875  elementStream >> Z;
876  return Z;
877  }
878 
880  return getNaturalIsotopicDistributions()->getIsotopicDistribution(Z);
881  }
882 
884  return getNaturalIsotopicDistributions()->drawRandomIsotope(Z);
885  }
886 
887  G4double getFermiMomentumConstant(const G4int /*A*/, const G4int /*Z*/) {
888  return constantFermiMomentum;
889  }
890 
892 // assert(Z>0 && A>0 && Z<=A);
894  const G4double rms = momentumRMS[Z][A];
895  return ((rms>0.) ? rms : momentumRMS[6][12]) * Math::sqrtFiveThirds;
896  } else
897  return getFermiMomentumConstant(A,Z);
898  }
899 
901 // assert(A>0);
902  static const G4double alphaParam = 259.416; // MeV/c
903  static const G4double betaParam = 152.824; // MeV/c
904  static const G4double gammaParam = 9.5157E-2;
905  return alphaParam - betaParam*std::exp(-gammaParam*((G4double)A));
906  }
907 
909 // assert(t==Proton || t==Neutron);
910  return rpCorrelationCoefficient[t];
911  }
912 
913  G4double getNeutronSkin() { return neutronSkin; }
914 
915  G4double getNeutronHalo() { return neutronHalo; }
916 
924 
926 // assert(isosp == -2 || isosp == 0 || isosp == 2);
927  if (isosp == -2) {
928  return PiMinus;
929  }
930  else if (isosp == 0) {
931  return PiZero;
932  }
933  else {
934  return PiPlus;
935  }
936  }
937 
939 // assert(isosp == -1 || isosp == 1);
940  if (isosp == -1) {
941  return Neutron;
942  }
943  else {
944  return Proton;
945  }
946  }
947 
949 // assert(isosp == -3 || isosp == -1 || isosp == 1 || isosp == 3);
950  if (isosp == -3) {
951  return DeltaMinus;
952  }
953  else if (isosp == -1) {
954  return DeltaZero;
955  }
956  else if (isosp == 1) {
957  return DeltaPlus;
958  }
959  else {
960  return DeltaPlusPlus;
961  }
962  }
963 
965 // assert(pt == PiPlus || pt == PiMinus || pt == PiZero || pt == Eta || pt == Omega || pt == EtaPrime);
966  if(pt == PiPlus) {
967  return piPlusWidth;
968  } else if(pt == PiMinus) {
969  return piMinusWidth;
970  } else if(pt == PiZero) {
971  return piZeroWidth;
972  } else if(pt == Eta) {
973  return etaWidth;
974  } else if(pt == Omega) {
975  return omegaWidth;
976  } else if(pt == EtaPrime) {
977  return etaPrimeWidth;
978  } else {
979  INCL_ERROR("getWidth : Unknown particle type." << '\n');
980  return 0.0;
981  }
982  }
983 
984 
985  } // namespace ParticleTable
986 } // namespace G4INCL
987 
#define INCL_FATAL(x)
const G4double sqrtThreeFifths
G4double(* ParticleMassFn)(const ParticleType)
G4double getSeparationEnergyReal(const ParticleType t, const G4int A, const G4int Z)
Return the real separation energy.
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
G4ThreadLocal ParticleMassFn getTableParticleMass
Static pointer to the mass function for particles.
G4ThreadLocal G4double minDeltaMassRndm
G4double getNeutronSkin()
Get the thickness of the neutron skin.
G4double getNeutronSeparationEnergy()
Getter for neutronSeparationEnergy.
G4double getNeutronSkin() const
Get the neutron-skin thickness.
std::string const & getINCLXXDataFilePath() const
G4ThreadLocal G4double minDeltaMass2
G4double getProtonSeparationEnergy()
Getter for protonSeparationEnergy.
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
const char * p
Definition: xmltok.h:285
G4double getTableQValue(const G4int A1, const G4int Z1, const G4int A2, const G4int Z2)
Get Q-value (in MeV/c^2)
#define INCL_ERROR(x)
G4int getChargeNumber(const ParticleType t)
Get charge number from particle type.
#define INCL_WARN(x)
G4double getFermiMomentum() const
Get the Fermi momentum.
G4double(* NuclearMassFn)(const G4int, const G4int)
G4double getFermiMomentumConstant(const G4int, const G4int)
Return the constant value of the Fermi momentum.
#define G4ThreadLocal
Definition: tls.hh:89
ParticleType getNucleonType(const G4int isosp)
Get the type of nucleon.
int G4int
Definition: G4Types.hh:78
void setProtonSeparationEnergy(const G4double s)
Setter for protonSeparationEnergy.
G4double(* FermiMomentumFn)(const G4int, const G4int)
SeparationEnergyType getSeparationEnergyType() const
Get the separation-energy type.
G4double getNuclearRadius(const ParticleType t, const G4int A, const G4int Z)
ParticleType getPionType(const G4int isosp)
Get the type of pion.
const XML_Char * s
Definition: expat.h:262
G4double getRealMass(const G4INCL::ParticleType t)
Get particle mass (in MeV/c^2)
G4int getMassNumber(const ParticleType t)
Get mass number from particle type.
G4IonTable * GetIonTable() const
double A(double temperature)
G4int parseElement(std::string pS)
Get the name of the element from the atomic number.
Class that stores isotopic abundances for a given element.
std::string getIUPACElementName(const G4int Z)
Get the name of an unnamed element from the IUPAC convention.
G4double getLargestNuclearRadius(const G4int A, const G4int Z)
const G4double sqrtFiveThirds
const G4int n
G4bool getUseRealMasses() const
Whether to use real masses.
void initialize(Config const *const theConfig)
Initialize the clustering model based on the Config object.
ParticleType getDeltaType(const G4int isosp)
Get the type of delta.
G4double(* SeparationEnergyFn)(const ParticleType, const G4int, const G4int)
Functions that encapsulate a mass table.
G4double getSurfaceDiffuseness(const ParticleType t, const G4int A, const G4int Z)
G4double getRadiusParameter(const ParticleType t, const G4int A, const G4int Z)
G4ThreadLocal SeparationEnergyFn getSeparationEnergy
Static pointer to the separation-energy function.
G4double getNeutronHalo() const
Get the neutron-halo size.
G4double getWidth(const ParticleType t)
Get particle width (in s)
FermiMomentumType getFermiMomentumType() const
Get the Fermi-momentum type.
std::string getElementName(const G4int Z)
Get the name of the element from the atomic number.
G4double getMaximumNuclearRadius(const ParticleType t, const G4int A, const G4int Z)
G4double GetPDGMass() const
const G4double effectiveDeltaWidth
static G4ParticleTable * GetParticleTable()
G4double getSeparationEnergyRealForLight(const ParticleType t, const G4int A, const G4int Z)
Return the real separation energy only for light nuclei.
G4double getINCLMass(const G4int A, const G4int Z)
Get INCL nuclear mass (in MeV/c^2)
G4ThreadLocal G4double minDeltaMass
G4int getIsospin(const ParticleType t)
Get the isospin of a particle.
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4ThreadLocal NuclearMassFn getTableMass
Static pointer to the mass function for nuclei.
const G4double Pf
Fermi momentum [MeV/c].
G4int drawRandomNaturalIsotope(const G4int Z)
static constexpr double MeV
Definition: G4SIunits.hh:214
G4double getTableSpeciesMass(const ParticleSpecies &p)
G4double getMomentumRMS(const G4int A, const G4int Z)
Return the RMS of the momentum distribution (light clusters)
#define INCL_DEFAULT_SEPARATION_ENERGY
IsotopicDistribution const & getNaturalIsotopicDistribution(const G4int Z)
G4double getRPCorrelationCoefficient(const ParticleType t)
Get the value of the r-p correlation coefficient.
G4int parseIUPACElement(std::string const &pS)
Parse a IUPAC element name.
const G4double effectiveDeltaMass
double G4double
Definition: G4Types.hh:76
#define INCL_DEBUG(x)
tuple c
Definition: test.py:13
G4double getSeparationEnergyINCL(const ParticleType t, const G4int, const G4int)
Return INCL&#39;s default separation energy.
G4double getFermiMomentumMassDependent(const G4int A, const G4int)
Return the value Fermi momentum from a fit.
G4double getRPCorrelationCoefficient(const ParticleType t) const
Get the r-p correlation coefficient.
std::string getShortName(const ParticleType t)
Get the short INCL name of the particle.
void initialize(Config const *const theConfig=0)
Initialize the particle table.
void setNeutronSeparationEnergy(const G4double s)
Setter for protonSeparationEnergy.
G4double getFermiMomentumConstantLight(const G4int A, const G4int Z)
Return the constant value of the Fermi momentum - special for light.
std::string getName(const ParticleType t)
Get the native INCL name of the particle.
G4double getNeutronHalo()
Get the size of the neutron halo.
G4ThreadLocal FermiMomentumFn getFermiMomentum