Geant4: G4INCLParticleTable.cc Source File

00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // INCL++ intra-nuclear cascade model
00027 // Pekka Kaitaniemi, CEA and Helsinki Institute of Physics
00028 // Davide Mancusi, CEA
00029 // Alain Boudard, CEA
00030 // Sylvie Leray, CEA
00031 // Joseph Cugnon, University of Liege
00032 //
00033 #define INCLXX_IN_GEANT4_MODE 1
00034 
00035 #include "globals.hh"
00036 
00037 #include "G4INCLParticleTable.hh"
00038 #include <algorithm>
00039 // #include <cassert>
00040 #include <cmath>
00041 #include <cctype>
00042 #include <sstream>
00043 
00044 #ifdef INCLXX_IN_GEANT4_MODE
00045 #include "G4PhysicalConstants.hh"
00046 #include "G4SystemOfUnits.hh"
00047 #endif
00048 
00049 namespace G4INCL {
00050 
00052   const NaturalIsotopicDistributions *ParticleTable::theNaturalIsotopicDistributions = NULL;
00053 
00055   ParticleTable::NuclearMassFn ParticleTable::getTableMass;
00057   ParticleTable::ParticleMassFn ParticleTable::getTableParticleMass;
00059   ParticleTable::SeparationEnergyFn ParticleTable::getSeparationEnergy;
00060 
00061   const G4double ParticleTable::theINCLNucleonMass = 938.2796;
00062   const G4double ParticleTable::theINCLPionMass = 138.0;
00063   G4double ParticleTable::protonMass = 0.0;
00064   G4double ParticleTable::neutronMass = 0.0;
00065   G4double ParticleTable::piPlusMass = 0.0;
00066   G4double ParticleTable::piMinusMass = 0.0;
00067   G4double ParticleTable::piZeroMass = 0.0;
00068 
00069   // e^2/(4 pi epsilon_0) [MeV fm]
00070   const G4double ParticleTable::eSquared = 1.439964;
00071 
00072   // Hard-coded values of the real particle masses (MeV/c^2)
00073   G4double ParticleTable::theRealProtonMass = 938.27203;
00074   G4double ParticleTable::theRealNeutronMass = 939.56536;
00075   G4double ParticleTable::theRealChargedPiMass = 139.57018;
00076   G4double ParticleTable::theRealPiZeroMass = 134.9766;
00077 
00078   const G4double ParticleTable::mediumDiffuseness[mediumNucleiTableSize] =
00079   {0.0,0.0,0.0,0.0,0.0,1.78,1.77,1.77,1.77,1.71,
00080     1.69,1.69,1.635,1.730,1.81,1.833,1.798,
00081     1.841,0.567,0.571, 0.560,0.549,0.550,0.551,
00082     0.580,0.575,0.569,0.537,0.0,0.0};
00083   const G4double ParticleTable::mediumRadius[mediumNucleiTableSize] =
00084   {0.0,0.0,0.0,0.0,0.0,0.334,0.327,0.479,0.631,0.838,
00085     0.811,1.07,1.403,1.335,1.25,1.544,1.498,1.513,
00086     2.58,2.77, 2.775,2.78,2.88,2.98,3.22,3.03,2.84,
00087     3.14,0.0,0.0};
00088   const G4double ParticleTable::positionRMS[clusterTableZSize][clusterTableASize] = {
00089 /*     A=   0     1     2     3     4     5     6     7     8     9    10    11    12 */
00090 /* 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},
00091 /* 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},
00092 /* 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},
00093 /* 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},
00094 /* 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},
00095 /* 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},
00096 /* 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},
00097 /* 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},
00098 /* 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}
00099   };
00100 
00101   const G4double ParticleTable::momentumRMS[clusterTableZSize][clusterTableASize] = {
00102 /*     A=   0     1     2     3     4     5     6     7     8     9    10    11    12 */
00103 /* 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},
00104 /* Z=1 */ {0.00, 0.00, 77.0, 110., 153., 100., 100., 100., 0.00, 0.00, 0.00, 0.00, 0.00},
00105 /* Z=2 */ {0.00, 0.00, 0.00, 110., 153., 100., 100., 100., 100., 100., 100., 0.00, 0.00},
00106 /* Z=3 */ {0.00, 0.00, 0.00, 0.00, 153., 100., 100., 100., 100., 100., 100., 100., 100.},
00107 /* Z=4 */ {0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 100., 100., 100., 100., 100., 100., 100.},
00108 /* Z=5 */ {0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 100., 100., 100., 100., 100., 100., 100.},
00109 /* Z=6 */ {0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 100., 100., 100., 100., 100.},
00110 /* 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.},
00111 /* 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.}
00112   };
00113 
00115   const std::string ParticleTable::elementTable[elementTableSize] = {
00116     "",
00117     "H",
00118     "He",
00119     "Li",
00120     "Be",
00121     "B",
00122     "C",
00123     "N",
00124     "O",
00125     "F",
00126     "Ne",
00127     "Na",
00128     "Mg",
00129     "Al",
00130     "Si",
00131     "P",
00132     "S",
00133     "Cl",
00134     "Ar",
00135     "K",
00136     "Ca",
00137     "Sc",
00138     "Ti",
00139     "V",
00140     "Cr",
00141     "Mn",
00142     "Fe",
00143     "Co",
00144     "Ni",
00145     "Cu",
00146     "Zn",
00147     "Ga",
00148     "Ge",
00149     "As",
00150     "Se",
00151     "Br",
00152     "Kr",
00153     "Rb",
00154     "Sr",
00155     "Y",
00156     "Zr",
00157     "Nb",
00158     "Mo",
00159     "Tc",
00160     "Ru",
00161     "Rh",
00162     "Pd",
00163     "Ag",
00164     "Cd",
00165     "In",
00166     "Sn",
00167     "Sb",
00168     "Te",
00169     "I",
00170     "Xe",
00171     "Cs",
00172     "Ba",
00173     "La",
00174     "Ce",
00175     "Pr",
00176     "Nd",
00177     "Pm",
00178     "Sm",
00179     "Eu",
00180     "Gd",
00181     "Tb",
00182     "Dy",
00183     "Ho",
00184     "Er",
00185     "Tm",
00186     "Yb",
00187     "Lu",
00188     "Hf",
00189     "Ta",
00190     "W",
00191     "Re",
00192     "Os",
00193     "Ir",
00194     "Pt",
00195     "Au",
00196     "Hg",
00197     "Tl",
00198     "Pb",
00199     "Bi",
00200     "Po",
00201     "At",
00202     "Rn",
00203     "Fr",
00204     "Ra",
00205     "Ac",
00206     "Th",
00207     "Pa",
00208     "U",
00209     "Np",
00210     "Pu",
00211     "Am",
00212     "Cm",
00213     "Bk",
00214     "Cf",
00215     "Es",
00216     "Fm",
00217     "Md",
00218     "No",
00219     "Lr",
00220     "Rf",
00221     "Db",
00222     "Sg",
00223     "Bh",
00224     "Hs",
00225     "Mt",
00226     "Ds",
00227     "Rg",
00228     "Cn"
00229   };
00230 
00232   const std::string ParticleTable::elementIUPACDigits = "nubtqphsoe";
00233 
00234   /* Table for cluster decays
00235    * Definition of "Stable": halflife > 1 ms
00236    *
00237    * These table includes decay data for clusters that INCL presently does
00238    * not produce. It can't hurt.
00239    *
00240    * Unphysical nuclides (A<Z) are marked as stable, but should never be
00241    * produced by INCL. If you find them in the output, something is fishy.
00242    */
00243   const ParticleTable::ClusterDecayType ParticleTable::clusterDecayMode[clusterTableZSize][clusterTableASize] =
00244   {
00245 /*                       A = 0              1               2               3               4                5               6                7               8               9             10             11             12 */
00246 /* Z =  0 */    {StableCluster, StableCluster,   NeutronDecay, NeutronUnbound, NeutronUnbound,  NeutronUnbound, NeutronUnbound,  NeutronUnbound, NeutronUnbound, NeutronUnbound,  NeutronUnbound, NeutronUnbound, NeutronUnbound},
00247 /* Z =  1 */    {StableCluster, StableCluster,  StableCluster,  StableCluster,   NeutronDecay, TwoNeutronDecay,   NeutronDecay, TwoNeutronDecay, NeutronUnbound, NeutronUnbound,  NeutronUnbound, NeutronUnbound, NeutronUnbound},
00248 /* Z =  2 */    {StableCluster, StableCluster,    ProtonDecay,  StableCluster,  StableCluster,    NeutronDecay,  StableCluster,    NeutronDecay,  StableCluster,   NeutronDecay, TwoNeutronDecay, NeutronUnbound, NeutronUnbound},
00249 /* Z =  3 */    {StableCluster, StableCluster,  StableCluster,  ProtonUnbound,    ProtonDecay,     ProtonDecay,  StableCluster,   StableCluster,  StableCluster,  StableCluster,    NeutronDecay,  StableCluster,   NeutronDecay},
00250 /* Z =  4 */    {StableCluster, StableCluster,  StableCluster,  StableCluster,  ProtonUnbound,     ProtonDecay, TwoProtonDecay,   StableCluster,     AlphaDecay,  StableCluster,   StableCluster,  StableCluster,  StableCluster},
00251 /* Z =  5 */    {StableCluster, StableCluster,  StableCluster,  StableCluster,  StableCluster,   ProtonUnbound, TwoProtonDecay,     ProtonDecay,  StableCluster,    ProtonDecay,   StableCluster,  StableCluster,  StableCluster},
00252 /* Z =  6 */    {StableCluster, StableCluster,  StableCluster,  StableCluster,  StableCluster,   StableCluster,  ProtonUnbound,   ProtonUnbound, TwoProtonDecay,  StableCluster,   StableCluster,  StableCluster,  StableCluster},
00253 /* Z =  7 */    {StableCluster, StableCluster,  StableCluster,  StableCluster,  StableCluster,   StableCluster,  StableCluster,   ProtonUnbound,  ProtonUnbound,  ProtonUnbound,     ProtonDecay,    ProtonDecay,  StableCluster},
00254 /* Z =  8 */    {StableCluster, StableCluster,  StableCluster,  StableCluster,  StableCluster,   StableCluster,  StableCluster,   StableCluster,  ProtonUnbound,  ProtonUnbound,   ProtonUnbound,  ProtonUnbound,    ProtonDecay}
00255   };
00256 
00260   const G4double ParticleTable::clusterPosFact[maxClusterMass+1] = {0.0, 1.0, 0.5,
00261                                                       0.33333, 0.25,
00262                                                       0.2, 0.16667,
00263                                                       0.14286, 0.125,
00264                                                       0.11111, 0.1,
00265                                                       0.09091, 0.083333};
00266 
00270   const G4double ParticleTable::clusterPosFact2[maxClusterMass+1] = {0.0, 1.0, 0.25,
00271                                                        0.11111, 0.0625,
00272                                                        0.04, 0.0277778,
00273                                                        0.020408, 0.015625,
00274                                                        0.012346, 0.01,
00275                                                        0.0082645, 0.0069444};
00276 
00277   const G4int ParticleTable::clusterZMin[maxClusterMass+1] = {0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3};
00278   const G4int ParticleTable::clusterZMax[maxClusterMass+1] = {0, 0, 1, 2, 3, 3, 5, 5, 6, 6, 7, 7, 8};
00279   const G4double ParticleTable::clusterPhaseSpaceCut[maxClusterMass+1] = {0.0, 70000.0, 180000.0,
00280                                                             90000.0, 90000.0,
00281                                                             128941.0 ,145607.0,
00282                                                             161365.0, 176389.0,
00283                                                             190798.0, 204681.0,
00284                                                             218109.0, 231135.0};
00285   const G4double ParticleTable::effectiveNucleonMass = 938.2796;
00286   const G4double ParticleTable::effectiveNucleonMass2 = 8.8036860777616e5;
00287   const G4double ParticleTable::effectiveDeltaMass = 1232.0;
00288   const G4double ParticleTable::effectivePionMass = 138.0;
00289   const G4double ParticleTable::effectiveDeltaDecayThreshold =
00290     ParticleTable::theRealNeutronMass + ParticleTable::theRealChargedPiMass
00291     + 0.5;
00292   const G4double ParticleTable::theINCLProtonSeparationEnergy = 6.83;
00293   const G4double ParticleTable::theINCLNeutronSeparationEnergy = ParticleTable::theINCLProtonSeparationEnergy;
00294   G4double ParticleTable::protonSeparationEnergy = theINCLProtonSeparationEnergy;
00295   G4double ParticleTable::neutronSeparationEnergy = theINCLNeutronSeparationEnergy;
00296 
00297 #ifdef INCLXX_IN_GEANT4_MODE
00298   G4IonTable *ParticleTable::theG4IonTable;
00299 #else
00300   std::vector< std::vector <G4bool> > ParticleTable::massTableMask;
00301   std::vector< std::vector <G4double> > ParticleTable::massTable;
00302 #endif
00303 
00304   void ParticleTable::initialize(Config const * const theConfig /*=0*/) {
00305     protonMass = theINCLNucleonMass;
00306     neutronMass = theINCLNucleonMass;
00307     piPlusMass = theINCLPionMass;
00308     piMinusMass = theINCLPionMass;
00309     piZeroMass = theINCLPionMass;
00310 
00311 #ifndef INCLXX_IN_GEANT4_MODE
00312     std::string dataFilePath;
00313     if(theConfig)
00314       dataFilePath = theConfig->getINCLXXDataFilePath();
00315     readRealMasses(dataFilePath);
00316 #endif
00317 
00318     if(theConfig && theConfig->getUseRealMasses()) {
00319       getTableMass = ParticleTable::getRealMass;
00320       getTableParticleMass = ParticleTable::getRealMass;
00321     } else {
00322       getTableMass = ParticleTable::getINCLMass;
00323       getTableParticleMass = ParticleTable::getINCLMass;
00324     }
00325 #ifdef INCLXX_IN_GEANT4_MODE
00326     G4ParticleTable *theG4ParticleTable = G4ParticleTable::GetParticleTable();
00327     theG4IonTable = theG4ParticleTable->GetIonTable();
00328     theRealProtonMass = theG4ParticleTable->FindParticle("proton")->GetPDGMass() / MeV;
00329     theRealNeutronMass = theG4ParticleTable->FindParticle("neutron")->GetPDGMass() / MeV;
00330     theRealChargedPiMass = theG4ParticleTable->FindParticle("pi+")->GetPDGMass() / MeV;
00331     theRealPiZeroMass = theG4ParticleTable->FindParticle("pi0")->GetPDGMass() / MeV;
00332 #endif
00333 
00334     // Initialise the separation-energy function
00335     if(!theConfig || theConfig->getSeparationEnergyType()==INCLSeparationEnergy)
00336       getSeparationEnergy = ParticleTable::getSeparationEnergyINCL;
00337     else if(theConfig->getSeparationEnergyType()==RealSeparationEnergy)
00338       getSeparationEnergy = ParticleTable::getSeparationEnergyReal;
00339     else if(theConfig->getSeparationEnergyType()==RealForLightSeparationEnergy)
00340       getSeparationEnergy = ParticleTable::getSeparationEnergyRealForLight;
00341     else {
00342       FATAL("Unrecognized separation-energy type in ParticleTable initialization: " << theConfig->getSeparationEnergyType() << std::endl);
00343     }
00344 
00345   }
00346 
00347   G4int ParticleTable::getIsospin(const ParticleType t) {
00348     // Actually this is the 3rd component of isospin (I_z) multiplied by 2!
00349     if(t == Proton) {
00350       return 1;
00351     } else if(t == Neutron) {
00352       return -1;
00353     } else if(t == PiPlus) {
00354       return 2;
00355     } else if(t == PiMinus) {
00356       return -2;
00357     } else if(t == PiZero) {
00358       return 0;
00359     } else if(t == DeltaPlusPlus) {
00360       return 3;
00361     } else if(t == DeltaPlus) {
00362       return 1;
00363     } else if(t == DeltaZero) {
00364       return -1;
00365     } else if(t == DeltaMinus) {
00366       return -3;
00367     }
00368 
00369     ERROR("Requested isospin of an unknown particle!");
00370     return -10; // Unknown
00371   }
00372 
00373   std::string ParticleTable::getShortName(const ParticleSpecies s) {
00374     if(s.theType==Composite)
00375       return getShortName(s.theA,s.theZ);
00376     else
00377       return getShortName(s.theType);
00378   }
00379 
00380   std::string ParticleTable::getName(const ParticleSpecies s) {
00381     if(s.theType==Composite)
00382       return getName(s.theA,s.theZ);
00383     else
00384       return getName(s.theType);
00385   }
00386 
00387   std::string ParticleTable::getName(const G4int A, const G4int Z) {
00388     std::stringstream stream;
00389     stream << getElementName(Z) << "-" << A;
00390     return stream.str();
00391   }
00392 
00393   std::string ParticleTable::getShortName(const G4int A, const G4int Z) {
00394     std::stringstream stream;
00395     stream << getElementName(Z);
00396     if(A>0)
00397       stream << A;
00398     return stream.str();
00399   }
00400 
00401   std::string ParticleTable::getName(const ParticleType p) {
00402     if(p == G4INCL::Proton) {
00403       return std::string("proton");
00404     } else if(p == G4INCL::Neutron) {
00405       return std::string("neutron");
00406     } else if(p == G4INCL::DeltaPlusPlus) {
00407       return std::string("delta++");
00408     } else if(p == G4INCL::DeltaPlus) {
00409       return std::string("delta+");
00410     } else if(p == G4INCL::DeltaZero) {
00411       return std::string("delta0");
00412     } else if(p == G4INCL::DeltaMinus) {
00413       return std::string("delta-");
00414     } else if(p == G4INCL::PiPlus) {
00415       return std::string("pi+");
00416     } else if(p == G4INCL::PiZero) {
00417       return std::string("pi0");
00418     } else if(p == G4INCL::PiMinus) {
00419       return std::string("pi-");
00420     } else if(p == G4INCL::Composite) {
00421       return std::string("composite");
00422     }
00423     return std::string("unknown");
00424   }
00425 
00426   std::string ParticleTable::getShortName(const ParticleType p) {
00427     if(p == G4INCL::Proton) {
00428       return std::string("p");
00429     } else if(p == G4INCL::Neutron) {
00430       return std::string("n");
00431     } else if(p == G4INCL::DeltaPlusPlus) {
00432       return std::string("d++");
00433     } else if(p == G4INCL::DeltaPlus) {
00434       return std::string("d+");
00435     } else if(p == G4INCL::DeltaZero) {
00436       return std::string("d0");
00437     } else if(p == G4INCL::DeltaMinus) {
00438       return std::string("d-");
00439     } else if(p == G4INCL::PiPlus) {
00440       return std::string("pi+");
00441     } else if(p == G4INCL::PiZero) {
00442       return std::string("pi0");
00443     } else if(p == G4INCL::PiMinus) {
00444       return std::string("pi-");
00445     } else if(p == G4INCL::Composite) {
00446       return std::string("comp");
00447     }
00448     return std::string("unknown");
00449   }
00450 
00451   G4double ParticleTable::getINCLMass(const ParticleType pt) {
00452     if(pt == Proton) {
00453       return protonMass;
00454     } else if(pt == Neutron) {
00455       return neutronMass;
00456     } else if(pt == PiPlus) {
00457       return piPlusMass;
00458     } else if(pt == PiMinus) {
00459       return piMinusMass;
00460     } else if(pt == PiZero) {
00461       return piZeroMass;
00462     } else {
00463       ERROR("ParticleTable::getMass : Unknown particle type." << std::endl);
00464       return 0.0;
00465     }
00466   }
00467 
00468   G4double ParticleTable::getRealMass(const ParticleType t) {
00469     switch(t) {
00470       case Proton:
00471         return theRealProtonMass;
00472         break;
00473       case Neutron:
00474         return theRealNeutronMass;
00475         break;
00476       case PiPlus:
00477       case PiMinus:
00478         return theRealChargedPiMass;
00479         break;
00480       case PiZero:
00481         return theRealPiZeroMass;
00482         break;
00483       default:
00484         ERROR("Particle::getRealMass : Unknown particle type." << std::endl);
00485         return 0.0;
00486         break;
00487     }
00488   }
00489 
00490   G4double ParticleTable::getRealMass(const G4int A, const G4int Z) {
00491 // assert(A>=0);
00492     // For nuclei with Z<0 or Z>A, assume that the exotic charge state is due to pions
00493     if(Z<0)
00494       return A*neutronMass - Z*getRealMass(PiMinus);
00495     else if(Z>A)
00496       return A*protonMass + (A-Z)*getRealMass(PiPlus);
00497     else if(Z==0)
00498       return A*getRealMass(Neutron);
00499     else if(A==Z)
00500       return A*getRealMass(Proton);
00501     else if(A>1) {
00502 #ifndef INCLXX_IN_GEANT4_MODE
00503       if(ParticleTable::hasMassTable(A,Z))
00504         return ParticleTable::massTable.at(Z).at(A);
00505       else {
00506         DEBUG("Real mass unavailable for isotope A=" << A << ", Z=" << Z
00507             << ", using Weizsaecker's formula"
00508             << std::endl);
00509         return getWeizsaeckerMass(A,Z);
00510       }
00511 #else
00512       return theG4IonTable->GetNucleusMass(Z,A) / MeV;
00513 #endif
00514     } else
00515       return 0.;
00516   }
00517 
00518   G4double ParticleTable::getINCLMass(const G4int A, const G4int Z) {
00519 // assert(A>=0);
00520     // For nuclei with Z<0 or Z>A, assume that the exotic charge state is due to pions
00521     if(Z<0)
00522       return A*neutronMass - Z*getINCLMass(PiMinus);
00523     else if(Z>A)
00524       return A*protonMass + (A-Z)*getINCLMass(PiPlus);
00525     else if(A>1)
00526       return Z*(protonMass - protonSeparationEnergy) + (A-Z)*(neutronMass - neutronSeparationEnergy);
00527     else if(A==1 && Z==0)
00528       return getINCLMass(Neutron);
00529     else if(A==1 && Z==1)
00530       return getINCLMass(Proton);
00531     else
00532       return 0.;
00533   }
00534 
00535   G4double ParticleTable::getNuclearRadius(const G4int A, const G4int Z) {
00536 // assert(A>0 && Z>=0 && Z<=A);
00537     if(A >= 19 || (A < 6 && A >= 2)) {
00538       // For large (Woods-Saxon or Modified Harmonic Oscillator) or small
00539       // (Gaussian) nuclei, the radius parameter is just the nuclear radius
00540       return getRadiusParameter(A,Z);
00541     } else if(A < clusterTableASize && Z < clusterTableZSize && A >= 6) {
00542       const G4double thisRMS = positionRMS[Z][A];
00543       if(thisRMS>0.0)
00544         return thisRMS;
00545       else {
00546         ERROR("ParticleTable::getRadiusParameter : Radius for nucleus A = " << A << " Z = " << Z << " is ==0.0" << std::endl);
00547         return 0.0;
00548       }
00549     } else if(A < 19) {
00550       const G4double theRadiusParameter = getRadiusParameter(A, Z);
00551       const G4double theDiffusenessParameter = getSurfaceDiffuseness(A, Z);
00552       // The formula yields the nuclear RMS radius based on the parameters of
00553       // the nuclear-density function
00554       return 1.581*theDiffusenessParameter*
00555         (2.+5.*theRadiusParameter)/(2.+3.*theRadiusParameter);
00556     } else {
00557       ERROR("ParticleTable::getNuclearRadius: No radius for nucleus A = " << A << " Z = " << Z << std::endl);
00558       return 0.0;
00559     }
00560   }
00561 
00562   G4double ParticleTable::getRadiusParameter(const G4int A, const G4int Z) {
00563 // assert(A>0 && Z>=0 && Z<=A);
00564     if(A >= 28) {
00565       // phenomenological radius fit
00566       return (2.745e-4 * A + 1.063) * std::pow(A, 1.0/3.0);
00567     } else if(A < 6 && A >= 2) {
00568       if(Z<clusterTableZSize) {
00569         const G4double thisRMS = positionRMS[Z][A];
00570         if(thisRMS>0.0)
00571           return thisRMS;
00572         else {
00573           ERROR("ParticleTable::getRadiusParameter : Radius for nucleus A = " << A << " Z = " << Z << " is ==0.0" << std::endl);
00574           return 0.0;
00575         }
00576       } else {
00577         ERROR("ParticleTable::getRadiusParameter : No radius for nucleus A = " << A << " Z = " << Z << std::endl);
00578         return 0.0;
00579       }
00580     } else if(A < 28 && A >= 6) {
00581       return mediumRadius[A-1];
00582       //      return 1.581*mediumDiffuseness[A-1]*(2.+5.*mediumRadius[A-1])/(2.+3.*mediumRadius[A-1]);
00583     } else {
00584       ERROR("ParticleTable::getRadiusParameter: No radius for nucleus A = " << A << " Z = " << Z << std::endl);
00585       return 0.0;
00586     }
00587   }
00588 
00589   G4double ParticleTable::getMaximumNuclearRadius(const G4int A, const G4int Z) {
00590     const G4double XFOISA = 8.0;
00591     if(A >= 19) {
00592       return getNuclearRadius(A,Z) + XFOISA * getSurfaceDiffuseness(A,Z);
00593     } else if(A < 19 && A >= 6) {
00594       return 5.5 + 0.3 * (G4double(A) - 6.0)/12.0;
00595     } else if(A >= 2) {
00596       return ParticleTable::getNuclearRadius(A, Z) + 4.5;
00597     } else {
00598       ERROR("ParticleTable::getMaximumNuclearRadius : No maximum radius for nucleus A = " << A << " Z = " << Z << std::endl);
00599       return 0.0;
00600     }
00601   }
00602 
00603 #ifdef INCLXX_IN_GEANT4_MODE
00604   G4double ParticleTable::getSurfaceDiffuseness(const G4int A, const G4int /*Z*/ ) {
00605 #else
00606   G4double ParticleTable::getSurfaceDiffuseness(const G4int A, const G4int Z) {
00607 #endif // INCLXX_IN_GEANT4_MODE
00608 
00609     if(A >= 28) {
00610       return 1.63e-4 * A + 0.510;
00611     } else if(A < 28 && A >= 19) {
00612       return mediumDiffuseness[A-1];
00613     } else if(A < 19 && A >= 6) {
00614       return mediumDiffuseness[A-1];
00615     } else if(A < 6 && A >= 2) {
00616       ERROR("ParticleTable::getSurfaceDiffuseness: was called for A = " << A << " Z = " << Z << std::endl);
00617       return 0.0;
00618     } else {
00619       ERROR("ParticleTable::getSurfaceDiffuseness: No diffuseness for nucleus A = " << A << " Z = " << Z << std::endl);
00620       return 0.0;
00621     }
00622   }
00623 
00624   std::string ParticleTable::getElementName(const G4int Z) {
00625     if(Z<1) {
00626       WARN("getElementName called with Z<1" << std::endl);
00627       return elementTable[0];
00628     } else if(Z<elementTableSize)
00629       return elementTable[Z];
00630     else
00631       return getIUPACElementName(Z);
00632   }
00633 
00634 #ifndef INCLXX_IN_GEANT4_MODE
00635   void ParticleTable::readRealMasses(std::string const &path) {
00636     // Clear the existing tables, if any
00637     massTableMask.clear();
00638     massTable.clear();
00639 
00640     // File name
00641     std::string fileName(path + "/walletlifetime.dat");
00642     DEBUG("Reading real nuclear masses from file " << fileName << std::endl);
00643 
00644     // Open the file stream
00645     std::ifstream massTableIn(fileName.c_str());
00646     if(!massTableIn.good()) {
00647       FATAL("Cannot open " << fileName << " data file." << std::endl);
00648       return;
00649     }
00650 
00651     // Read in the data
00652     unsigned int Z, A;
00653     const G4double amu = 931.494061; // atomic mass unit in MeV/c^2
00654     const G4double eMass = 0.5109988; // electron mass in MeV/c^2
00655     G4double excess;
00656     massTableIn >> A >> Z >> excess;
00657     do {
00658       // Dynamically determine the table size
00659       if(Z>=massTable.size()) {
00660         massTable.resize(Z+1);
00661         massTableMask.resize(Z+1);
00662       }
00663       if(A>=massTable[Z].size()) {
00664         massTable[Z].resize(A+1, 0.);
00665         massTableMask[Z].resize(A+1, false);
00666       }
00667 
00668       massTable.at(Z).at(A) = A*amu + excess - Z*eMass;
00669       massTableMask.at(Z).at(A) = true;
00670 
00671       massTableIn >> A >> Z >> excess;
00672     } while(massTableIn.good());
00673 
00674   }
00675 #endif
00676 
00677   std::string ParticleTable::getIUPACElementName(const G4int Z) {
00678     std::stringstream elementStream;
00679     elementStream << Z;
00680     std::string elementName = elementStream.str();
00681     std::transform(elementName.begin(), elementName.end(), elementName.begin(), ParticleTable::intToIUPAC);
00682     elementName[0] = std::toupper(elementName.at(0));
00683     return elementName;
00684   }
00685 
00686   G4int ParticleTable::parseIUPACElement(std::string const &s) {
00687     // Normalise to lower case
00688     std::string elementName(s);
00689     std::transform(elementName.begin(), elementName.end(), elementName.begin(), ::tolower);
00690     // Return 0 if the element name contains anything but IUPAC digits
00691     if(elementName.find_first_not_of(elementIUPACDigits)!=std::string::npos)
00692       return 0;
00693     std::transform(elementName.begin(), elementName.end(), elementName.begin(), ParticleTable::iupacToInt);
00694     std::stringstream elementStream(elementName);
00695     G4int Z;
00696     elementStream >> Z;
00697     return Z;
00698   }
00699 
00700 }