Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // INCL++ intra-nuclear cascade model 26 // INCL++ intra-nuclear cascade model 27 // Alain Boudard, CEA-Saclay, France << 27 // Pekka Kaitaniemi, CEA and Helsinki Institute of Physics 28 // Joseph Cugnon, University of Liege, Belgium << 28 // Davide Mancusi, CEA 29 // Jean-Christophe David, CEA-Saclay, France << 29 // Alain Boudard, CEA 30 // Pekka Kaitaniemi, CEA-Saclay, France, and H << 30 // Sylvie Leray, CEA 31 // Sylvie Leray, CEA-Saclay, France << 31 // Joseph Cugnon, University of Liege 32 // Davide Mancusi, CEA-Saclay, France << 33 // 32 // 34 #define INCLXX_IN_GEANT4_MODE 1 33 #define INCLXX_IN_GEANT4_MODE 1 35 34 36 #include "globals.hh" 35 #include "globals.hh" 37 36 38 #include "G4INCLParticleTable.hh" 37 #include "G4INCLParticleTable.hh" 39 #include "G4INCLNuclearMassTable.hh" << 40 #include <algorithm> 38 #include <algorithm> 41 // #include <cassert> 39 // #include <cassert> 42 #include <cmath> 40 #include <cmath> 43 #include <cctype> 41 #include <cctype> 44 #include <sstream> 42 #include <sstream> 45 #ifdef INCLXX_IN_GEANT4_MODE << 46 #include "G4SystemOfUnits.hh" << 47 #endif << 48 43 49 #ifdef INCLXX_IN_GEANT4_MODE 44 #ifdef INCLXX_IN_GEANT4_MODE 50 #include "G4PhysicalConstants.hh" 45 #include "G4PhysicalConstants.hh" 51 #include "G4SystemOfUnits.hh" 46 #include "G4SystemOfUnits.hh" 52 #endif 47 #endif 53 48 54 namespace G4INCL { 49 namespace G4INCL { 55 50 56 namespace ParticleTable { << 51 /// \brief Static instance of the NaturalIsotopicAbundances class 57 << 52 const NaturalIsotopicDistributions *ParticleTable::theNaturalIsotopicDistributions = NULL; 58 namespace { << 59 << 60 /// \brief Static instance of the Natura << 61 const NaturalIsotopicDistributions *theN << 62 << 63 const G4double theINCLNucleonMass = 938. << 64 const G4double theINCLPionMass = 138.0; << 65 const G4double theINCLLambdaMass = 1115. << 66 //const G4double theINCLSigmaMass = 1197 << 67 //const G4double theINCLKaonMass = 497.6 << 68 const G4double theINCLEtaMass = 547.862; << 69 const G4double theINCLOmegaMass = 782.65 << 70 const G4double theINCLEtaPrimeMass = 957 << 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 << 78 G4ThreadLocal G4double SigmaZeroMass = 0 << 79 G4ThreadLocal G4double SigmaMinusMass = << 80 G4ThreadLocal G4double LambdaMass = 0.0; << 81 G4ThreadLocal G4double XiMinusMass = 0.0 << 82 G4ThreadLocal G4double XiZeroMass = 0.0; << 83 G4ThreadLocal G4double antiProtonMass = << 84 G4ThreadLocal G4double antiNeutronMass = << 85 G4ThreadLocal G4double antiSigmaPlusMass << 86 G4ThreadLocal G4double antiSigmaZeroMass << 87 G4ThreadLocal G4double antiSigmaMinusMas << 88 G4ThreadLocal G4double antiLambdaMass = << 89 G4ThreadLocal G4double antiXiMinusMass = << 90 G4ThreadLocal G4double antiXiZeroMass = << 91 G4ThreadLocal G4double KPlusMass = 0.0; << 92 G4ThreadLocal G4double KZeroMass = 0.0; << 93 G4ThreadLocal G4double KZeroBarMass = 0. << 94 G4ThreadLocal G4double KShortMass = 0.0; << 95 G4ThreadLocal G4double KLongMass = 0.0; << 96 G4ThreadLocal G4double KMinusMass = 0.0; << 97 G4ThreadLocal G4double etaMass = 0.0; << 98 G4ThreadLocal G4double omegaMass = 0.0; << 99 G4ThreadLocal G4double etaPrimeMass = 0. << 100 G4ThreadLocal G4double photonMass = 0.0; << 101 << 102 // Hard-coded values of the real particl << 103 G4ThreadLocal G4double theRealProtonMass << 104 G4ThreadLocal G4double theRealNeutronMas << 105 G4ThreadLocal G4double theRealChargedPiM << 106 G4ThreadLocal G4double theRealPiZeroMass << 107 G4ThreadLocal G4double theRealLambdaMass << 108 G4ThreadLocal G4double theRealSigmaPlusM << 109 G4ThreadLocal G4double theRealSigmaZeroM << 110 G4ThreadLocal G4double theRealSigmaMinus << 111 G4ThreadLocal G4double theRealAntiProton << 112 G4ThreadLocal G4double theRealXiMinusMas << 113 G4ThreadLocal G4double theRealXiZeroMass << 114 G4ThreadLocal G4double theRealAntiNeutro << 115 G4ThreadLocal G4double theRealAntiLambda << 116 G4ThreadLocal G4double theRealAntiSigmaP << 117 G4ThreadLocal G4double theRealAntiSigmaZ << 118 G4ThreadLocal G4double theRealAntiSigmaM << 119 G4ThreadLocal G4double theRealAntiXiMinu << 120 G4ThreadLocal G4double theRealAntiXiZero << 121 G4ThreadLocal G4double theRealChargedKao << 122 G4ThreadLocal G4double theRealNeutralKao << 123 G4ThreadLocal G4double theRealEtaMass = << 124 G4ThreadLocal G4double theRealOmegaMass << 125 G4ThreadLocal G4double theRealEtaPrimeMa << 126 G4ThreadLocal G4double theRealPhotonMass << 127 << 128 // Width (second) << 129 const G4double theChargedPiWidth = 2.603 << 130 const G4double thePiZeroWidth = 8.52e-17 << 131 const G4double theEtaWidth = 5.025e-19; << 132 const G4double theOmegaWidth = 7.7528e-2 << 133 const G4double theEtaPrimeWidth = 3.3243 << 134 const G4double theChargedKaonWidth = 1.2 << 135 const G4double theKShortWidth = 8.954e-1 << 136 const G4double theKLongWidth = 5.116e-08 << 137 const G4double theLambdaWidth = 2.632e-1 << 138 const G4double theSigmaPlusWidth = 8.018 << 139 const G4double theSigmaZeroWidth = 7.4e- << 140 const G4double theSigmaMinusWidth = 1.47 << 141 //const G4double theXiMinusWidth = 1.639 << 142 //const G4double theXiZeroWidth = 2.90e- << 143 //const G4double theAntiLambdaWidth = 2. << 144 //const G4double theAntiSigmaPlusWidth = << 145 //const G4double theAntiSigmaZeroWidth = << 146 //const G4double theAntiSigmaMinusWidth << 147 //const G4double theAntiXiMinusWidth = 1 << 148 //const G4double theAntiXiZeroWidth = 2. << 149 G4ThreadLocal G4double piPlusWidth = 0.0 << 150 G4ThreadLocal G4double piMinusWidth = 0. << 151 G4ThreadLocal G4double piZeroWidth = 0.0 << 152 G4ThreadLocal G4double etaWidth = 0.0; << 153 G4ThreadLocal G4double omegaWidth = 0.0; << 154 G4ThreadLocal G4double etaPrimeWidth = 0 << 155 G4ThreadLocal G4double LambdaWidth = 0.0 << 156 G4ThreadLocal G4double SigmaPlusWidth = << 157 G4ThreadLocal G4double SigmaZeroWidth = << 158 G4ThreadLocal G4double SigmaMinusWidth = << 159 G4ThreadLocal G4double KPlusWidth = 0.0; << 160 G4ThreadLocal G4double KMinusWidth = 0.0 << 161 G4ThreadLocal G4double KShortWidth = 0.0 << 162 G4ThreadLocal G4double KLongWidth = 0.0; << 163 G4ThreadLocal G4double XiMinusWidth = 0. << 164 G4ThreadLocal G4double XiZeroWidth = 0.0 << 165 G4ThreadLocal G4double antiLambdaWidth = << 166 G4ThreadLocal G4double antiSigmaZeroWidt << 167 G4ThreadLocal G4double antiSigmaMinusWid << 168 G4ThreadLocal G4double antiSigmaPlusWidt << 169 G4ThreadLocal G4double antiXiZeroWidth = << 170 G4ThreadLocal G4double antiXiMinusWidth << 171 << 172 const G4int mediumNucleiTableSize = 30; << 173 << 174 const G4double mediumDiffuseness[mediumN << 175 {0.0,0.0,0.0,0.0,0.0,1.78,1.77,1.77,1.69 << 176 1.69,1.72,1.635,1.730,1.81,1.833,1.798 << 177 1.93,0.567,0.571, 0.560,0.549,0.550,0. << 178 0.580,0.575,0.569,0.537,0.0,0.0}; << 179 const G4double mediumRadius[mediumNuclei << 180 {0.0,0.0,0.0,0.0,0.0,0.334,0.327,0.479,0 << 181 0.811,0.84,1.403,1.335,1.25,1.544,1.49 << 182 2.58,2.77, 2.775,2.78,2.88,2.98,3.22,3 << 183 3.14,0.0,0.0}; << 184 << 185 const G4double positionRMS[clusterTableZ << 186 /* A= 0 1 2 3 4 << 187 /* Z=0 */ {-1.0, -1.0, -1.0, -1.0, -1. << 188 /* Z=1 */ {-1.0, -1.0, 2.10, 1.80, 1.7 << 189 /* Z=2 */ {-1.0, -1.0, -1.0, 1.80, 1.6 << 190 /* Z=3 */ {-1.0, -1.0, -1.0, -1.0, 1.7 << 191 /* Z=4 */ {-1.0, -1.0, -1.0, -1.0, -1. << 192 /* Z=5 */ {-1.0, -1.0, -1.0, -1.0, -1. << 193 /* Z=6 */ {-1.0, -1.0, -1.0, -1.0, -1. << 194 /* Z=7 */ {-1.0, -1.0, -1.0, -1.0, -1. << 195 /* Z=8 */ {-1.0, -1.0, -1.0, -1.0, -1. << 196 }; << 197 << 198 const G4double momentumRMS[clusterTableZ << 199 /* A= 0 1 2 3 4 << 200 /* Z=0 */ {-1.0, -1.0, -1.0, -1.0, -1. << 201 /* Z=1 */ {-1.0, -1.0, 77.0, 110., 153 << 202 /* Z=2 */ {-1.0, -1.0, -1.0, 110., 153 << 203 /* Z=3 */ {-1.0, -1.0, -1.0, -1.0, 153 << 204 /* Z=4 */ {-1.0, -1.0, -1.0, -1.0, -1. << 205 /* Z=5 */ {-1.0, -1.0, -1.0, -1.0, -1. << 206 /* Z=6 */ {-1.0, -1.0, -1.0, -1.0, -1. << 207 /* Z=7 */ {-1.0, -1.0, -1.0, -1.0, -1. << 208 /* Z=8 */ {-1.0, -1.0, -1.0, -1.0, -1. << 209 }; << 210 << 211 const G4int elementTableSize = 113; // u << 212 << 213 /// \brief Table of chemical element nam << 214 const std::string elementTable[elementTa << 215 "", << 216 "H", << 217 "He", << 218 "Li", << 219 "Be", << 220 "B", << 221 "C", << 222 "N", << 223 "O", << 224 "F", << 225 "Ne", << 226 "Na", << 227 "Mg", << 228 "Al", << 229 "Si", << 230 "P", << 231 "S", << 232 "Cl", << 233 "Ar", << 234 "K", << 235 "Ca", << 236 "Sc", << 237 "Ti", << 238 "V", << 239 "Cr", << 240 "Mn", << 241 "Fe", << 242 "Co", << 243 "Ni", << 244 "Cu", << 245 "Zn", << 246 "Ga", << 247 "Ge", << 248 "As", << 249 "Se", << 250 "Br", << 251 "Kr", << 252 "Rb", << 253 "Sr", << 254 "Y", << 255 "Zr", << 256 "Nb", << 257 "Mo", << 258 "Tc", << 259 "Ru", << 260 "Rh", << 261 "Pd", << 262 "Ag", << 263 "Cd", << 264 "In", << 265 "Sn", << 266 "Sb", << 267 "Te", << 268 "I", << 269 "Xe", << 270 "Cs", << 271 "Ba", << 272 "La", << 273 "Ce", << 274 "Pr", << 275 "Nd", << 276 "Pm", << 277 "Sm", << 278 "Eu", << 279 "Gd", << 280 "Tb", << 281 "Dy", << 282 "Ho", << 283 "Er", << 284 "Tm", << 285 "Yb", << 286 "Lu", << 287 "Hf", << 288 "Ta", << 289 "W", << 290 "Re", << 291 "Os", << 292 "Ir", << 293 "Pt", << 294 "Au", << 295 "Hg", << 296 "Tl", << 297 "Pb", << 298 "Bi", << 299 "Po", << 300 "At", << 301 "Rn", << 302 "Fr", << 303 "Ra", << 304 "Ac", << 305 "Th", << 306 "Pa", << 307 "U", << 308 "Np", << 309 "Pu", << 310 "Am", << 311 "Cm", << 312 "Bk", << 313 "Cf", << 314 "Es", << 315 "Fm", << 316 "Md", << 317 "No", << 318 "Lr", << 319 "Rf", << 320 "Db", << 321 "Sg", << 322 "Bh", << 323 "Hs", << 324 "Mt", << 325 "Ds", << 326 "Rg", << 327 "Cn" << 328 }; << 329 << 330 /// \brief Digit names to compose IUPAC << 331 const std::string elementIUPACDigits = " << 332 << 333 #define INCL_DEFAULT_SEPARATION_ENERGY 6.83 << 334 const G4double theINCLProtonSeparationEn << 335 const G4double theINCLNeutronSeparationE << 336 const G4double theINCLLambdaSeparationEn << 337 //const G4double theINCLantiProtonSepara << 338 const G4double theINCLantiProtonSeparati << 339 G4ThreadLocal G4double protonSeparationE << 340 G4ThreadLocal G4double neutronSeparation << 341 G4ThreadLocal G4double lambdaSeparationE << 342 //G4ThreadLocal G4double antiprotonSepar << 343 //G4ThreadLocal G4double antiprotonSepar << 344 #undef INCL_DEFAULT_SEPARATION_ENERGY << 345 << 346 G4ThreadLocal G4double rpCorrelationCoef << 347 53 348 G4ThreadLocal G4double neutronSkin = 0.0 << 54 /// \brief Static pointer to the mass function for nuclei 349 G4ThreadLocal G4double neutronHalo = 0.0 << 55 ParticleTable::NuclearMassFn ParticleTable::getTableMass; >> 56 /// \brief Static pointer to the mass function for particles >> 57 ParticleTable::ParticleMassFn ParticleTable::getTableParticleMass; >> 58 /// \brief Static pointer to the separation-energy function >> 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 >> 114 /// \brief Table of chemical element names >> 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 >> 231 /// \brief Digit names to compose IUPAC element names >> 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 >> 257 /** >> 258 * Precomputed factor 1.0/A >> 259 */ >> 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 >> 267 /** >> 268 * Precomputed factor (1.0/A)^2 >> 269 */ >> 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; 350 296 351 #ifdef INCLXX_IN_GEANT4_MODE 297 #ifdef INCLXX_IN_GEANT4_MODE 352 G4ThreadLocal G4IonTable *theG4IonTable; << 298 G4IonTable *ParticleTable::theG4IonTable; >> 299 #else >> 300 std::vector< std::vector <G4bool> > ParticleTable::massTableMask; >> 301 std::vector< std::vector <G4double> > ParticleTable::massTable; 353 #endif 302 #endif 354 303 355 /// \brief Default value for constant Fe << 304 void ParticleTable::initialize(Config const * const theConfig /*=0*/) { 356 G4ThreadLocal G4double constantFermiMome << 305 protonMass = theINCLNucleonMass; 357 << 306 neutronMass = theINCLNucleonMass; 358 /// \brief Transform a IUPAC char to an << 307 piPlusMass = theINCLPionMass; 359 char iupacToInt(char c) { << 308 piMinusMass = theINCLPionMass; 360 return (char)(((G4int)'0')+elementIUPA << 309 piZeroMass = theINCLPionMass; 361 } << 362 << 363 /// \brief Transform an integer digit (r << 364 char intToIUPAC(char n) { return element << 365 << 366 /// \brief Get the singleton instance of << 367 const NaturalIsotopicDistributions *getN << 368 if(!theNaturalIsotopicDistributions) << 369 theNaturalIsotopicDistributions = ne << 370 return theNaturalIsotopicDistributions << 371 } << 372 << 373 } // namespace << 374 << 375 void initialize(Config const * const theCo << 376 protonMass = theINCLNucleonMass; << 377 neutronMass = theINCLNucleonMass; << 378 piPlusMass = theINCLPionMass; << 379 piMinusMass = theINCLPionMass; << 380 piZeroMass = theINCLPionMass; << 381 << 382 etaMass = theINCLEtaMass; << 383 omegaMass = theINCLOmegaMass; << 384 etaPrimeMass = theINCLEtaPrimeMass; << 385 photonMass = theINCLPhotonMass; << 386 << 387 SigmaPlusMass = theRealSigmaPlusMass; << 388 SigmaMinusMass = theRealSigmaMinusMass; << 389 SigmaZeroMass = theRealSigmaZeroMass; << 390 LambdaMass = theINCLLambdaMass; << 391 KPlusMass = theRealChargedKaonMass; << 392 KZeroMass = theRealNeutralKaonMass; << 393 KZeroBarMass = theRealNeutralKaonMass; << 394 KShortMass = theRealNeutralKaonMass; << 395 KLongMass = theRealNeutralKaonMass; << 396 KMinusMass = theRealChargedKaonMass; << 397 << 398 antiProtonMass = theRealAntiProtonMass; << 399 XiZeroMass = theRealXiZeroMass; << 400 XiMinusMass = theRealXiMinusMass; << 401 antiNeutronMass = theRealAntiNeutronMass << 402 antiSigmaPlusMass = theRealAntiSigmaPlus << 403 antiSigmaMinusMass = theRealAntiSigmaMin << 404 antiSigmaZeroMass = theRealAntiSigmaZero << 405 antiLambdaMass = theRealAntiLambdaMass; << 406 antiXiZeroMass = theRealAntiXiZeroMass; << 407 antiXiMinusMass = theRealAntiXiMinusMass << 408 << 409 if(theConfig && theConfig->getUseRealMas << 410 getTableMass = getRealMass; << 411 getTableParticleMass = getRealMass; << 412 } else { << 413 getTableMass = getINCLMass; << 414 getTableParticleMass = getINCLMass; << 415 } << 416 310 417 #ifndef INCLXX_IN_GEANT4_MODE 311 #ifndef INCLXX_IN_GEANT4_MODE 418 std::string dataFilePath; << 312 std::string dataFilePath; 419 if(theConfig) << 313 if(theConfig) 420 dataFilePath = theConfig->getINCLXXDat << 314 dataFilePath = theConfig->getINCLXXDataFilePath(); 421 NuclearMassTable::initialize(dataFilePat << 315 readRealMasses(dataFilePath); 422 #endif << 423 << 424 #ifdef INCLXX_IN_GEANT4_MODE << 425 G4ParticleTable *theG4ParticleTable = G4 << 426 theG4IonTable = theG4ParticleTable->GetI << 427 theRealProtonMass = theG4ParticleTable-> << 428 theRealNeutronMass = theG4ParticleTable- << 429 theRealChargedPiMass = theG4ParticleTabl << 430 theRealPiZeroMass = theG4ParticleTable-> << 431 << 432 theRealEtaMass = theG4ParticleTable->Fin << 433 theRealOmegaMass = theG4ParticleTable->F << 434 theRealEtaPrimeMass = theG4ParticleTable << 435 theRealPhotonMass = theG4ParticleTable-> << 436 << 437 theRealSigmaPlusMass = theG4ParticleTabl << 438 theRealSigmaZeroMass = theG4ParticleTabl << 439 theRealSigmaMinusMass = theG4ParticleTab << 440 theRealLambdaMass = theG4ParticleTable-> << 441 theRealChargedKaonMass = theG4ParticleTa << 442 theRealNeutralKaonMass = theG4ParticleTa << 443 << 444 theRealAntiProtonMass = theG4ParticleTab << 445 theRealAntiNeutronMass = theG4ParticleTa << 446 theRealXiZeroMass = theG4ParticleTable-> << 447 theRealXiMinusMass = theG4ParticleTable- << 448 theRealAntiSigmaPlusMass = theG4Particle << 449 theRealAntiSigmaZeroMass = theG4Particle << 450 theRealAntiSigmaMinusMass = theG4Particl << 451 theRealAntiLambdaMass = theG4ParticleTab << 452 theRealAntiXiZeroMass = theG4ParticleTab << 453 theRealAntiXiMinusMass = theG4ParticleTa << 454 #endif 316 #endif 455 317 456 minDeltaMass = theRealNeutronMass + theR << 318 if(theConfig && theConfig->getUseRealMasses()) { 457 minDeltaMass2 = minDeltaMass*minDeltaMas << 319 getTableMass = ParticleTable::getRealMass; 458 minDeltaMassRndm = std::atan((minDeltaMa << 320 getTableParticleMass = ParticleTable::getRealMass; 459 << 321 } else { 460 piPlusWidth = theChargedPiWidth; << 322 getTableMass = ParticleTable::getINCLMass; 461 piMinusWidth = theChargedPiWidth; << 323 getTableParticleMass = ParticleTable::getINCLMass; 462 piZeroWidth = thePiZeroWidth; << 324 } 463 etaWidth = theEtaWidth; << 464 omegaWidth = theOmegaWidth; << 465 etaPrimeWidth = theEtaPrimeWidth; << 466 << 467 SigmaMinusWidth = theSigmaMinusWidth; << 468 SigmaPlusWidth = theSigmaPlusWidth; << 469 SigmaZeroWidth = theSigmaZeroWidth; << 470 LambdaWidth = theLambdaWidth; << 471 KPlusWidth = theChargedKaonWidth; << 472 KMinusWidth = theChargedKaonWidth; << 473 KShortWidth = theKShortWidth; << 474 KLongWidth = theKLongWidth; << 475 << 476 // Initialise HFB tables << 477 #ifdef INCLXX_IN_GEANT4_MODE 325 #ifdef INCLXX_IN_GEANT4_MODE 478 HFB::initialize(); << 326 G4ParticleTable *theG4ParticleTable = G4ParticleTable::GetParticleTable(); 479 #else << 327 theG4IonTable = theG4ParticleTable->GetIonTable(); 480 HFB::initialize(dataFilePath); << 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; 481 #endif 332 #endif 482 333 483 // Initialise the separation-energy func << 334 // Initialise the separation-energy function 484 if(!theConfig || theConfig->getSeparatio << 335 if(!theConfig || theConfig->getSeparationEnergyType()==INCLSeparationEnergy) 485 getSeparationEnergy = getSeparationEne << 336 getSeparationEnergy = ParticleTable::getSeparationEnergyINCL; 486 else if(theConfig->getSeparationEnergyTy << 337 else if(theConfig->getSeparationEnergyType()==RealSeparationEnergy) 487 getSeparationEnergy = getSeparationEne << 338 getSeparationEnergy = ParticleTable::getSeparationEnergyReal; 488 else if(theConfig->getSeparationEnergyTy << 339 else if(theConfig->getSeparationEnergyType()==RealForLightSeparationEnergy) 489 getSeparationEnergy = getSeparationEne << 340 getSeparationEnergy = ParticleTable::getSeparationEnergyRealForLight; 490 else { << 341 else { 491 INCL_FATAL("Unrecognized separation-en << 342 FATAL("Unrecognized separation-energy type in ParticleTable initialization: " << theConfig->getSeparationEnergyType() << std::endl); 492 return; << 343 } 493 } << 344 494 << 345 } 495 // Initialise the Fermi-momentum functio << 346 496 if(!theConfig || theConfig->getFermiMome << 347 G4int ParticleTable::getIsospin(const ParticleType t) { 497 getFermiMomentum = ParticleTable::getF << 348 // Actually this is the 3rd component of isospin (I_z) multiplied by 2! 498 if(theConfig) { << 349 if(t == Proton) { 499 const G4double aFermiMomentum = theC << 350 return 1; 500 if(aFermiMomentum>0.) << 351 } else if(t == Neutron) { 501 constantFermiMomentum = aFermiMome << 352 return -1; 502 else << 353 } else if(t == PiPlus) { 503 constantFermiMomentum = PhysicalCo << 354 return 2; 504 } else { << 355 } else if(t == PiMinus) { 505 constantFermiMomentum = PhysicalCons << 356 return -2; 506 } << 357 } else if(t == PiZero) { 507 } else if(theConfig->getFermiMomentumTyp << 358 return 0; 508 getFermiMomentum = ParticleTable::getF << 359 } else if(t == DeltaPlusPlus) { 509 else if(theConfig->getFermiMomentumType( << 360 return 3; 510 getFermiMomentum = ParticleTable::getF << 361 } else if(t == DeltaPlus) { 511 else { << 362 return 1; 512 INCL_FATAL("Unrecognized Fermi-momentu << 363 } else if(t == DeltaZero) { 513 return; << 364 return -1; 514 } << 365 } else if(t == DeltaMinus) { 515 << 366 return -3; 516 // Initialise the r-p correlation coeffi << 367 } 517 std::fill(rpCorrelationCoefficient, rpCo << 368 518 if(theConfig) { << 369 ERROR("Requested isospin of an unknown particle!"); 519 rpCorrelationCoefficient[Proton] = the << 370 return -10; // Unknown 520 rpCorrelationCoefficient[Neutron] = th << 371 } 521 } << 372 522 << 373 std::string ParticleTable::getShortName(const ParticleSpecies s) { 523 // Initialise the neutron-skin parameter << 374 if(s.theType==Composite) 524 if(theConfig) { << 375 return getShortName(s.theA,s.theZ); 525 neutronSkin = theConfig->getNeutronSki << 376 else 526 neutronHalo = theConfig->getNeutronHal << 377 return getShortName(s.theType); 527 } << 378 } 528 << 379 529 } << 380 std::string ParticleTable::getName(const ParticleSpecies s) { 530 << 381 if(s.theType==Composite) 531 G4int getIsospin(const ParticleType t) { << 382 return getName(s.theA,s.theZ); 532 // Actually this is the 3rd component of << 383 else 533 if(t == Proton) { << 384 return getName(s.theType); 534 return 1; << 385 } 535 } else if(t == Neutron) { << 386 536 return -1; << 387 std::string ParticleTable::getName(const G4int A, const G4int Z) { 537 } else if(t == PiPlus) { << 388 std::stringstream stream; 538 return 2; << 389 stream << getElementName(Z) << "-" << A; 539 } else if(t == PiMinus) { << 390 return stream.str(); 540 return -2; << 391 } 541 } else if(t == PiZero) { << 392 542 return 0; << 393 std::string ParticleTable::getShortName(const G4int A, const G4int Z) { 543 } else if(t == DeltaPlusPlus) { << 394 std::stringstream stream; 544 return 3; << 395 stream << getElementName(Z); 545 } else if(t == DeltaPlus) { << 396 if(A>0) 546 return 1; << 397 stream << A; 547 } else if(t == DeltaZero) { << 398 return stream.str(); 548 return -1; << 399 } 549 } else if(t == DeltaMinus) { << 400 550 return -3; << 401 std::string ParticleTable::getName(const ParticleType p) { 551 } else if(t == Lambda) { << 402 if(p == G4INCL::Proton) { 552 return 0; << 403 return std::string("proton"); 553 } else if(t == SigmaPlus) { << 404 } else if(p == G4INCL::Neutron) { 554 return 2; << 405 return std::string("neutron"); 555 } else if(t == SigmaZero) { << 406 } else if(p == G4INCL::DeltaPlusPlus) { 556 return 0; << 407 return std::string("delta++"); 557 } else if(t == SigmaMinus) { << 408 } else if(p == G4INCL::DeltaPlus) { 558 return -2; << 409 return std::string("delta+"); 559 } else if(t == KPlus) { << 410 } else if(p == G4INCL::DeltaZero) { 560 return 1; << 411 return std::string("delta0"); 561 } else if(t == KZero) { << 412 } else if(p == G4INCL::DeltaMinus) { 562 return -1; << 413 return std::string("delta-"); 563 } else if(t == KZeroBar) { << 414 } else if(p == G4INCL::PiPlus) { 564 return 1; << 415 return std::string("pi+"); 565 } else if(t == KShort) { << 416 } else if(p == G4INCL::PiZero) { 566 return 0; << 417 return std::string("pi0"); 567 } else if(t == KLong) { << 418 } else if(p == G4INCL::PiMinus) { 568 return 0; << 419 return std::string("pi-"); 569 } else if(t == KMinus) { << 420 } else if(p == G4INCL::Composite) { 570 return -1; << 421 return std::string("composite"); 571 } else if(t == Eta) { << 422 } 572 return 0; << 423 return std::string("unknown"); 573 } else if(t == Omega) { << 424 } 574 return 0; << 425 575 } else if(t == EtaPrime) { << 426 std::string ParticleTable::getShortName(const ParticleType p) { 576 return 0; << 427 if(p == G4INCL::Proton) { 577 } else if(t == Photon) { << 428 return std::string("p"); 578 return 0; << 429 } else if(p == G4INCL::Neutron) { 579 } else if(t == antiProton) { << 430 return std::string("n"); 580 return -1; << 431 } else if(p == G4INCL::DeltaPlusPlus) { 581 } else if(t == XiMinus) { << 432 return std::string("d++"); 582 return -1; << 433 } else if(p == G4INCL::DeltaPlus) { 583 } else if(t == XiZero) { << 434 return std::string("d+"); 584 return 1; << 435 } else if(p == G4INCL::DeltaZero) { 585 } else if(t == antiNeutron) { << 436 return std::string("d0"); 586 return 1; << 437 } else if(p == G4INCL::DeltaMinus) { 587 } else if(t == antiLambda) { << 438 return std::string("d-"); 588 return 0; << 439 } else if(p == G4INCL::PiPlus) { 589 } else if(t == antiSigmaPlus) { << 440 return std::string("pi+"); 590 return -2; << 441 } else if(p == G4INCL::PiZero) { 591 } else if(t == antiSigmaZero) { << 442 return std::string("pi0"); 592 return 0; << 443 } else if(p == G4INCL::PiMinus) { 593 } else if(t == antiSigmaMinus) { << 444 return std::string("pi-"); 594 return 2; << 445 } else if(p == G4INCL::Composite) { 595 } else if(t == antiXiMinus) { << 446 return std::string("comp"); 596 return 1; << 447 } 597 } else if(t == antiXiZero) { << 448 return std::string("unknown"); 598 return -1; << 449 } 599 } << 450 600 INCL_ERROR("Requested isospin of an unkn << 451 G4double ParticleTable::getINCLMass(const ParticleType pt) { 601 return -10; // Unknown << 452 if(pt == Proton) { 602 } << 453 return protonMass; 603 << 454 } else if(pt == Neutron) { 604 std::string getShortName(const ParticleSpe << 455 return neutronMass; 605 if(sp.theType==Composite && sp.theS == 0 << 456 } else if(pt == PiPlus) { 606 return getShortName(sp.theA,sp.theZ); << 457 return piPlusMass; 607 else if(sp.theType==Composite) << 458 } else if(pt == PiMinus) { 608 return getName(sp.theA,sp.theZ,sp.theS << 459 return piMinusMass; 609 else << 460 } else if(pt == PiZero) { 610 return getShortName(sp.theType); << 461 return piZeroMass; 611 } << 462 } else { 612 << 463 ERROR("ParticleTable::getMass : Unknown particle type." << std::endl); 613 std::string getName(const ParticleSpecies << 464 return 0.0; 614 if(sp.theType==Composite && sp.theS == 0 << 465 } 615 return getName(sp.theA,sp.theZ); << 466 } 616 else if(sp.theType==Composite) << 467 617 return getName(sp.theA,sp.theZ,sp.theS << 468 G4double ParticleTable::getRealMass(const ParticleType t) { 618 else << 469 switch(t) { 619 return getName(sp.theType); << 470 case Proton: 620 } << 471 return theRealProtonMass; 621 << 472 break; 622 std::string getName(const G4int A, const G << 473 case Neutron: 623 std::stringstream stream; << 474 return theRealNeutronMass; 624 stream << getElementName(Z) << "-" << A; << 475 break; 625 return stream.str(); << 476 case PiPlus: 626 } << 477 case PiMinus: 627 << 478 return theRealChargedPiMass; 628 std::string getName(const G4int A, const G << 479 break; 629 std::stringstream stream; << 480 case PiZero: 630 if(S >= 0) // S < 0 for hypernuclei << 481 return theRealPiZeroMass; 631 return getName(A, Z); << 482 break; 632 else if(S == -1) << 483 default: 633 stream << getElementName(Z) << "-" << << 484 ERROR("Particle::getRealMass : Unknown particle type." << std::endl); 634 else << 635 stream << getElementName(Z) << "-" << << 636 return stream.str(); << 637 } << 638 << 639 std::string getShortName(const G4int A, co << 640 std::stringstream stream; << 641 stream << getElementName(Z); << 642 if(A>0) << 643 stream << A; << 644 return stream.str(); << 645 } << 646 << 647 std::string getName(const ParticleType p) << 648 if(p == G4INCL::Proton) { << 649 return std::string("proton"); << 650 } else if(p == G4INCL::Neutron) { << 651 return std::string("neutron"); << 652 } else if(p == G4INCL::DeltaPlusPlus) { << 653 return std::string("delta++"); << 654 } else if(p == G4INCL::DeltaPlus) { << 655 return std::string("delta+"); << 656 } else if(p == G4INCL::DeltaZero) { << 657 return std::string("delta0"); << 658 } else if(p == G4INCL::DeltaMinus) { << 659 return std::string("delta-"); << 660 } else if(p == G4INCL::PiPlus) { << 661 return std::string("pi+"); << 662 } else if(p == G4INCL::PiZero) { << 663 return std::string("pi0"); << 664 } else if(p == G4INCL::PiMinus) { << 665 return std::string("pi-"); << 666 } else if(p == G4INCL::Lambda) { << 667 return std::string("lambda"); << 668 } else if(p == G4INCL::SigmaPlus) { << 669 return std::string("sigma+"); << 670 } else if(p == G4INCL::SigmaZero) { << 671 return std::string("sigma0"); << 672 } else if(p == G4INCL::SigmaMinus) { << 673 return std::string("sigma-"); << 674 } else if(p == G4INCL::antiProton) { << 675 return std::string("antiproton"); << 676 } else if(p == G4INCL::XiMinus) { << 677 return std::string("xi-"); << 678 } else if(p == G4INCL::XiZero) { << 679 return std::string("xi0"); << 680 } else if(p == G4INCL::antiNeutron) { << 681 return std::string("antineutron"); << 682 } else if(p == G4INCL::antiSigmaPlus) { << 683 return std::string("antisigma+"); << 684 } else if(p == G4INCL::antiSigmaZero) { << 685 return std::string("antisigma0"); << 686 } else if(p == G4INCL::antiSigmaMinus) { << 687 return std::string("antisigma-"); << 688 } else if(p == G4INCL::antiLambda) { << 689 return std::string("antilambda"); << 690 } else if(p == G4INCL::antiXiMinus) { << 691 return std::string("antixi-"); << 692 } else if(p == G4INCL::antiXiZero) { << 693 return std::string("antixi0"); << 694 } else if(p == G4INCL::KPlus) { << 695 return std::string("kaon+"); << 696 } else if(p == G4INCL::KZero) { << 697 return std::string("kaon0"); << 698 } else if(p == G4INCL::KZeroBar) { << 699 return std::string("kaon0bar"); << 700 } else if(p == G4INCL::KMinus) { << 701 return std::string("kaon-"); << 702 } else if(p == G4INCL::KShort) { << 703 return std::string("kaonshort"); << 704 } else if(p == G4INCL::KLong) { << 705 return std::string("kaonlong"); << 706 } else if(p == G4INCL::Composite) { << 707 return std::string("composite"); << 708 } else if(p == G4INCL::Eta) { << 709 return std::string("eta"); << 710 } else if(p == G4INCL::Omega) { << 711 return std::string("omega"); << 712 } else if(p == G4INCL::EtaPrime) { << 713 return std::string("etaprime"); << 714 } else if(p == G4INCL::Photon) { << 715 return std::string("photon"); << 716 } << 717 return std::string("unknown"); << 718 } << 719 << 720 std::string getShortName(const ParticleTyp << 721 if(p == G4INCL::Proton) { << 722 return std::string("p"); << 723 } else if(p == G4INCL::Neutron) { << 724 return std::string("n"); << 725 } else if(p == G4INCL::DeltaPlusPlus) { << 726 return std::string("d++"); << 727 } else if(p == G4INCL::DeltaPlus) { << 728 return std::string("d+"); << 729 } else if(p == G4INCL::DeltaZero) { << 730 return std::string("d0"); << 731 } else if(p == G4INCL::DeltaMinus) { << 732 return std::string("d-"); << 733 } else if(p == G4INCL::PiPlus) { << 734 return std::string("pi+"); << 735 } else if(p == G4INCL::PiZero) { << 736 return std::string("pi0"); << 737 } else if(p == G4INCL::PiMinus) { << 738 return std::string("pi-"); << 739 } else if(p == G4INCL::Lambda) { << 740 return std::string("l"); << 741 } else if(p == G4INCL::SigmaPlus) { << 742 return std::string("s+"); << 743 } else if(p == G4INCL::SigmaZero) { << 744 return std::string("s0"); << 745 } else if(p == G4INCL::SigmaMinus) { << 746 return std::string("s-"); << 747 } else if(p == G4INCL::antiProton) { << 748 return std::string("pb"); << 749 } else if(p == G4INCL::XiMinus) { << 750 return std::string("x-"); << 751 } else if(p == G4INCL::XiZero) { << 752 return std::string("x0"); << 753 } else if(p == G4INCL::antiNeutron) { << 754 return std::string("nb"); << 755 } else if(p == G4INCL::antiSigmaPlus) { << 756 return std::string("s+b"); << 757 } else if(p == G4INCL::antiSigmaZero) { << 758 return std::string("s0b"); << 759 } else if(p == G4INCL::antiSigmaMinus) { << 760 return std::string("s-b"); << 761 } else if(p == G4INCL::antiLambda) { << 762 return std::string("lb"); << 763 } else if(p == G4INCL::antiXiMinus) { << 764 return std::string("x-b"); << 765 } else if(p == G4INCL::antiXiZero) { << 766 return std::string("x0b"); << 767 } else if(p == G4INCL::KPlus) { << 768 return std::string("k+"); << 769 } else if(p == G4INCL::KZero) { << 770 return std::string("k0"); << 771 } else if(p == G4INCL::KZeroBar) { << 772 return std::string("k0b"); << 773 } else if(p == G4INCL::KMinus) { << 774 return std::string("k-"); << 775 } else if(p == G4INCL::KShort) { << 776 return std::string("ks"); << 777 } else if(p == G4INCL::KLong) { << 778 return std::string("kl"); << 779 } else if(p == G4INCL::Composite) { << 780 return std::string("comp"); << 781 } else if(p == G4INCL::Eta) { << 782 return std::string("eta"); << 783 } else if(p == G4INCL::Omega) { << 784 return std::string("omega"); << 785 } else if(p == G4INCL::EtaPrime) { << 786 return std::string("etap"); << 787 } else if(p == G4INCL::Photon) { << 788 return std::string("photon"); << 789 } << 790 return std::string("unknown"); << 791 } << 792 << 793 G4double getINCLMass(const ParticleType pt << 794 if(pt == Proton) { << 795 return protonMass; << 796 } else if(pt == Neutron) { << 797 return neutronMass; << 798 } else if(pt == PiPlus) { << 799 return piPlusMass; << 800 } else if(pt == PiMinus) { << 801 return piMinusMass; << 802 } else if(pt == PiZero) { << 803 return piZeroMass; << 804 } else if(pt == SigmaPlus) { << 805 return SigmaPlusMass; << 806 } else if(pt == SigmaMinus) { << 807 return SigmaMinusMass; << 808 } else if(pt == SigmaZero) { << 809 return SigmaZeroMass; << 810 } else if(pt == Lambda) { << 811 return LambdaMass; << 812 } else if(pt == antiProton) { << 813 return antiProtonMass; << 814 } else if(pt == XiMinus) { << 815 return XiMinusMass; << 816 } else if(pt == XiZero) { << 817 return XiZeroMass; << 818 } else if(pt == antiNeutron) { << 819 return antiNeutronMass; << 820 } else if(pt == antiSigmaPlus) { << 821 return antiSigmaPlusMass; << 822 } else if(pt == antiSigmaMinus) { << 823 return antiSigmaMinusMass; << 824 } else if(pt == antiSigmaZero) { << 825 return antiSigmaZeroMass; << 826 } else if(pt == antiLambda) { << 827 return antiLambdaMass; << 828 } else if(pt == antiXiMinus) { << 829 return antiXiMinusMass; << 830 } else if(pt == antiXiZero) { << 831 return antiXiZeroMass; << 832 } else if(pt == KPlus) { << 833 return KPlusMass; << 834 } else if(pt == KZero) { << 835 return KZeroMass; << 836 } else if(pt == KZeroBar) { << 837 return KZeroBarMass; << 838 } else if(pt == KMinus) { << 839 return KMinusMass; << 840 } else if(pt == KShort) { << 841 return KShortMass; << 842 } else if(pt == KLong) { << 843 return KLongMass; << 844 } else if(pt == Eta) { << 845 return etaMass; << 846 } else if(pt == Omega) { << 847 return omegaMass; << 848 } else if(pt == EtaPrime) { << 849 return etaPrimeMass; << 850 } else if(pt == Photon) { << 851 return photonMass; << 852 } else { << 853 INCL_ERROR("getMass : Unknown particle << 854 return 0.0; 485 return 0.0; 855 } << 486 break; 856 } 487 } 857 << 488 } 858 G4double getRealMass(const ParticleType t) << 489 859 switch(t) { << 490 G4double ParticleTable::getRealMass(const G4int A, const G4int Z) { 860 case Proton: << 861 return theRealProtonMass; << 862 break; << 863 case Neutron: << 864 return theRealNeutronMass; << 865 break; << 866 case PiPlus: << 867 case PiMinus: << 868 return theRealChargedPiMass; << 869 break; << 870 case PiZero: << 871 return theRealPiZeroMass; << 872 break; << 873 case Eta: << 874 return theRealEtaMass; << 875 break; << 876 case Omega: << 877 return theRealOmegaMass; << 878 break; << 879 case EtaPrime: << 880 return theRealEtaPrimeMass; << 881 break; << 882 case Photon: << 883 return theRealPhotonMass; << 884 break; << 885 case Lambda: << 886 return theRealLambdaMass; << 887 break; << 888 case KPlus: << 889 case KMinus: << 890 return theRealChargedKaonMass; << 891 break; << 892 case KZero: << 893 case KZeroBar: << 894 case KShort: << 895 case KLong: << 896 return theRealNeutralKaonMass; << 897 break; << 898 case SigmaPlus: << 899 return theRealSigmaPlusMass; << 900 break; << 901 case SigmaZero: << 902 return theRealSigmaZeroMass; << 903 break; << 904 case SigmaMinus: << 905 return theRealSigmaMinusMass; << 906 break; << 907 case antiProton: << 908 return theRealAntiProtonMass; << 909 break; << 910 case XiMinus: << 911 return theRealXiMinusMass; << 912 break; << 913 case XiZero: << 914 return theRealXiZeroMass; << 915 break; << 916 case antiNeutron: << 917 return theRealAntiNeutronMass; << 918 break; << 919 case antiSigmaPlus: << 920 return theRealAntiSigmaPlusMass; << 921 break; << 922 case antiSigmaZero: << 923 return theRealAntiSigmaZeroMass; << 924 break; << 925 case antiSigmaMinus: << 926 return theRealAntiSigmaMinusMass; << 927 break; << 928 case antiXiMinus: << 929 return theRealAntiXiMinusMass; << 930 break; << 931 case antiXiZero: << 932 return theRealAntiXiZeroMass; << 933 break; << 934 case antiLambda: << 935 return theRealAntiLambdaMass; << 936 break; << 937 default: << 938 INCL_ERROR("Particle::getRealMass : << 939 return 0.0; << 940 break; << 941 } << 942 } << 943 << 944 G4double getRealMass(const G4int A, const << 945 // assert(A>=0); 491 // assert(A>=0); 946 // For nuclei with Z<0 or Z>A, assume th << 492 // For nuclei with Z<0 or Z>A, assume that the exotic charge state is due to pions 947 if(Z<0 && S<0) << 493 if(Z<0) 948 return (A+S)*theRealNeutronMass - S*La << 494 return A*neutronMass - Z*getRealMass(PiMinus); 949 else if(Z>A && S<0) << 495 else if(Z>A) 950 return (A+S)*theRealProtonMass - S*Lam << 496 return A*protonMass + (A-Z)*getRealMass(PiPlus); 951 if(Z<0) << 497 else if(Z==0) 952 return (A)*theRealNeutronMass - Z*getR << 498 return A*getRealMass(Neutron); 953 else if(Z>A) << 499 else if(A==Z) 954 return (A)*theRealProtonMass + (A-Z)*g << 500 return A*getRealMass(Proton); 955 else if(Z==0 && S==0) << 501 else if(A>1) { 956 return A*theRealNeutronMass; << 957 else if(A==Z) << 958 return A*theRealProtonMass; << 959 else if(Z==0 && S<0) << 960 return (A+S)*theRealNeutronMass-S*Lamb << 961 else if(A>1) { << 962 #ifndef INCLXX_IN_GEANT4_MODE 502 #ifndef INCLXX_IN_GEANT4_MODE 963 return ::G4INCL::NuclearMassTable::get << 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 } 964 #else 511 #else 965 if(S<0) return theG4IonTable->GetNucle << 512 return theG4IonTable->GetNucleusMass(Z,A) / MeV; 966 else return theG4IonTable->GetNucle << 967 #endif 513 #endif 968 } else << 514 } else 969 return 0.; << 515 return 0.; 970 } << 516 } 971 517 972 G4double getINCLMass(const G4int A, const << 518 G4double ParticleTable::getINCLMass(const G4int A, const G4int Z) { 973 // assert(A>=0); 519 // assert(A>=0); 974 // For nuclei with Z<0 or Z>A, assume th << 520 // For nuclei with Z<0 or Z>A, assume that the exotic charge state is due to pions 975 // Note that S<0 for lambda << 521 if(Z<0) 976 if(Z<0 && S<0) << 522 return A*neutronMass - Z*getINCLMass(PiMinus); 977 return (A+S)*neutronMass - S*LambdaMas << 523 else if(Z>A) 978 else if(Z>A && S<0) << 524 return A*protonMass + (A-Z)*getINCLMass(PiPlus); 979 return (A+S)*protonMass - S*LambdaMass << 525 else if(A>1) 980 else if(Z<0) << 526 return Z*(protonMass - protonSeparationEnergy) + (A-Z)*(neutronMass - neutronSeparationEnergy); 981 return (A)*neutronMass - Z*getINCLMass << 527 else if(A==1 && Z==0) 982 else if(Z>A) << 528 return getINCLMass(Neutron); 983 return (A)*protonMass + (A-Z)*getINCLM << 529 else if(A==1 && Z==1) 984 else if(A>1 && S<0) << 530 return getINCLMass(Proton); 985 return Z*(protonMass - protonSeparatio << 531 else 986 else if(A>1) << 532 return 0.; 987 return Z*(protonMass - protonSeparatio << 533 } 988 else if(A==1 && Z==0 && S==0) << 534 989 return getINCLMass(Neutron); << 535 G4double ParticleTable::getNuclearRadius(const G4int A, const G4int Z) { 990 else if(A==1 && Z==1 && S==0) << 536 // assert(A>0 && Z>=0 && Z<=A); 991 return getINCLMass(Proton); << 537 if(A >= 19 || (A < 6 && A >= 2)) { 992 else if(A==1 && Z==0 && S==-1) << 538 // For large (Woods-Saxon or Modified Harmonic Oscillator) or small 993 return getINCLMass(Lambda); << 539 // (Gaussian) nuclei, the radius parameter is just the nuclear radius 994 else << 540 return getRadiusParameter(A,Z); 995 return 0.; << 541 } else if(A < clusterTableASize && Z < clusterTableZSize && A >= 6) { 996 } << 542 const G4double thisRMS = positionRMS[Z][A]; 997 << 543 if(thisRMS>0.0) 998 G4double getTableQValue(const G4int A1, co << 544 return thisRMS; 999 return getTableMass(A1,Z1,S1) + getTable << 545 else { 1000 } << 546 ERROR("ParticleTable::getRadiusParameter : Radius for nucleus A = " << A << " Z = " << Z << " is ==0.0" << std::endl); 1001 << 547 return 0.0; 1002 G4double getTableQValue(const G4int A1, c << 1003 return getTableMass(A1,Z1,S1) + getTabl << 1004 } << 1005 << 1006 G4double getTableSpeciesMass(const Partic << 1007 if(p.theType == Composite) << 1008 return (*getTableMass)(p.theA, p.theZ << 1009 else << 1010 return (*getTableParticleMass)(p.theT << 1011 } << 1012 << 1013 G4int getMassNumber(const ParticleType t) << 1014 << 1015 switch(t) { << 1016 case Proton: << 1017 case Neutron: << 1018 case DeltaPlusPlus: << 1019 case DeltaPlus: << 1020 case DeltaZero: << 1021 case DeltaMinus: << 1022 case SigmaPlus: << 1023 case SigmaZero: << 1024 case SigmaMinus: << 1025 case Lambda: << 1026 case XiZero: << 1027 case XiMinus: << 1028 return 1; << 1029 break; << 1030 case antiProton: << 1031 case antiNeutron: << 1032 case antiSigmaPlus: << 1033 case antiSigmaZero: << 1034 case antiSigmaMinus: << 1035 case antiLambda: << 1036 case antiXiZero: << 1037 case antiXiMinus: << 1038 return -1; << 1039 break; << 1040 case PiPlus: << 1041 case PiMinus: << 1042 case PiZero: << 1043 case KPlus: << 1044 case KZero: << 1045 case KZeroBar: << 1046 case KShort: << 1047 case KLong: << 1048 case KMinus: << 1049 case Eta: << 1050 case Omega: << 1051 case EtaPrime: << 1052 case Photon: << 1053 return 0; << 1054 break; << 1055 default: << 1056 return 0; << 1057 break; << 1058 } << 1059 } << 1060 << 1061 G4int getChargeNumber(const ParticleType << 1062 switch(t) { << 1063 case DeltaPlusPlus: << 1064 return 2; << 1065 break; << 1066 case Proton: << 1067 case DeltaPlus: << 1068 case PiPlus: << 1069 case SigmaPlus: << 1070 case KPlus: << 1071 case antiSigmaMinus: << 1072 case antiXiMinus: << 1073 return 1; << 1074 break; << 1075 case Neutron: << 1076 case DeltaZero: << 1077 case PiZero: << 1078 case SigmaZero: << 1079 case Lambda: << 1080 case KZero: << 1081 case KZeroBar: << 1082 case KShort: << 1083 case KLong: << 1084 case Eta: << 1085 case Omega: << 1086 case EtaPrime: << 1087 case Photon: << 1088 case XiZero: << 1089 case antiNeutron: << 1090 case antiLambda: << 1091 case antiSigmaZero: << 1092 case antiXiZero: << 1093 return 0; << 1094 break; << 1095 case DeltaMinus: << 1096 case PiMinus: << 1097 case SigmaMinus: << 1098 case KMinus: << 1099 case antiProton: << 1100 case XiMinus: << 1101 case antiSigmaPlus: << 1102 return -1; << 1103 break; << 1104 default: << 1105 return 0; << 1106 break; << 1107 } << 1108 } << 1109 << 1110 G4int getStrangenessNumber(const Particle << 1111 switch(t) { << 1112 case DeltaPlusPlus: << 1113 case DeltaPlus: << 1114 case DeltaZero: << 1115 case DeltaMinus: << 1116 case Proton: << 1117 case Neutron: << 1118 case PiPlus: << 1119 case PiZero: << 1120 case PiMinus: << 1121 case Eta: << 1122 case Omega: << 1123 case EtaPrime: << 1124 case Photon: << 1125 case antiProton: << 1126 case antiNeutron: << 1127 return 0; << 1128 break; << 1129 case XiMinus: << 1130 case XiZero: << 1131 case antiXiMinus: << 1132 case antiXiZero: << 1133 return 2; << 1134 break; << 1135 case antiLambda: << 1136 case antiSigmaPlus: << 1137 case antiSigmaZero: << 1138 case antiSigmaMinus: << 1139 return 1; << 1140 break; << 1141 case Lambda: << 1142 case SigmaPlus: << 1143 case SigmaZero: << 1144 case SigmaMinus: << 1145 case KZeroBar: << 1146 case KMinus: << 1147 return -1; << 1148 break; << 1149 case KPlus: << 1150 case KZero: << 1151 return 1; << 1152 break; << 1153 case KShort: << 1154 return 0; << 1155 break; << 1156 case KLong: << 1157 return 0; << 1158 break; << 1159 default: << 1160 return 0; << 1161 break; << 1162 } 548 } 1163 } << 549 } else if(A < 19) { 1164 << 550 const G4double theRadiusParameter = getRadiusParameter(A, Z); 1165 G4double getNuclearRadius(const ParticleT << 551 const G4double theDiffusenessParameter = getSurfaceDiffuseness(A, Z); 1166 // assert(A>=0); << 552 // The formula yields the nuclear RMS radius based on the parameters of 1167 if(A > 19 || (A < 6 && A >= 2)) { << 553 // the nuclear-density function 1168 // For large (Woods-Saxon or Modified << 554 return 1.581*theDiffusenessParameter* 1169 // (Gaussian) nuclei, the radius para << 555 (2.+5.*theRadiusParameter)/(2.+3.*theRadiusParameter); 1170 return getRadiusParameter(t,A,Z); << 556 } else { 1171 } else if(A < clusterTableASize && Z>=0 << 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) { 1172 const G4double thisRMS = positionRMS[ 569 const G4double thisRMS = positionRMS[Z][A]; 1173 if(thisRMS>0.0) 570 if(thisRMS>0.0) 1174 return thisRMS; 571 return thisRMS; 1175 else { 572 else { 1176 INCL_DEBUG("getNuclearRadius: Radiu << 573 ERROR("ParticleTable::getRadiusParameter : Radius for nucleus A = " << A << " Z = " << Z << " is ==0.0" << std::endl); 1177 << "returning radius for << 574 return 0.0; 1178 return positionRMS[6][12]; << 1179 } << 1180 } else if(A <= 19) { << 1181 const G4double theRadiusParameter = g << 1182 const G4double theDiffusenessParamete << 1183 // The formula yields the nuclear RMS << 1184 // the nuclear-density function << 1185 return 1.225*theDiffusenessParameter* << 1186 std::sqrt((2.+5.*theRadiusParameter << 1187 } else { << 1188 INCL_ERROR("getNuclearRadius: No radi << 1189 return 0.0; << 1190 } << 1191 } << 1192 << 1193 G4double getLargestNuclearRadius(const G4 << 1194 return Math::max(getNuclearRadius(Proto << 1195 } << 1196 << 1197 G4double getRadiusParameter(const Particl << 1198 // assert(A>0); << 1199 if(A > 19) { << 1200 // radius fit for lambdas << 1201 if(t==Lambda){ << 1202 G4double r0 = (1.128+0.439*std::pow( << 1203 return r0; << 1204 } << 1205 // phenomenological radius fit << 1206 G4double r0 = (2.745e-4 * A + 1.063) << 1207 // HFB calculations << 1208 if(getRPCorrelationCoefficient(t)<1.) << 1209 G4double r0hfb = HFB::getRadiusParam << 1210 if(r0hfb>0.)r0 = r0hfb; << 1211 } << 1212 // << 1213 if(t==Neutron) << 1214 r0 += neutronSkin; << 1215 return r0; << 1216 } else if(A < 6 && A >= 2) { << 1217 if(Z<clusterTableZSize && Z>=0) { << 1218 const G4double thisRMS = positionRM << 1219 if(thisRMS>0.0) << 1220 return thisRMS; << 1221 else { << 1222 INCL_DEBUG("getRadiusParameter: R << 1223 << "returning radius f << 1224 return positionRMS[6][12]; << 1225 } << 1226 } else { << 1227 INCL_DEBUG("getRadiusParameter: Rad << 1228 << "returning radius for << 1229 return positionRMS[6][12]; << 1230 } << 1231 } else if(A <= 19 && A >= 6) { << 1232 if(t==Lambda){ << 1233 G4double r0 = (1.128+0.439*std::pow( << 1234 return r0; << 1235 } << 1236 // HFB calculations << 1237 if(getRPCorrelationCoefficient(t)<1.) << 1238 G4double r0hfb = HFB::getSurfaceDiff << 1239 if(r0hfb>0.)return r0hfb; << 1240 } << 1241 return mediumRadius[A-1]; << 1242 // return 1.581*mediumDiffusenes << 1243 } else { << 1244 INCL_ERROR("getRadiusParameter: No ra << 1245 return 0.0; << 1246 } << 1247 } << 1248 << 1249 G4double getMaximumNuclearRadius(const Pa << 1250 const G4double XFOISA = 8.0; << 1251 if(A > 19) { << 1252 return getNuclearRadius(t,A,Z) + XFOI << 1253 } else if(A <= 19 && A >= 6) { << 1254 return 5.5 + 0.3 * (G4double(A) - 6.0 << 1255 } else if(A >= 2) { << 1256 return getNuclearRadius(t, A, Z) + 4. << 1257 } else { << 1258 INCL_ERROR("getMaximumNuclearRadius : << 1259 return 0.0; << 1260 } << 1261 } << 1262 << 1263 G4double getSurfaceDiffuseness(const Part << 1264 if(A > 19) { << 1265 // phenomenological fit << 1266 G4double a = 1.63e-4 * A + 0.510; << 1267 // HFB calculations << 1268 if(getRPCorrelationCoefficient(t)<1.) << 1269 G4double ahfb = HFB::getSurfaceDiff << 1270 if(ahfb>0.)a=ahfb; << 1271 } << 1272 // << 1273 if(t==Lambda){ << 1274 // Like for neutrons << 1275 G4double ahfb = HFB::getSurfaceDiff << 1276 if(ahfb>0.)a=ahfb; << 1277 } << 1278 if(t==Neutron) << 1279 a += neutronHalo; << 1280 return a; << 1281 } else if(A <= 19 && A >= 6) { << 1282 // HFB calculations << 1283 if(getRPCorrelationCoefficient(t)<1.) << 1284 G4double ahfb = HFB::getRadiusParam << 1285 if(ahfb>0.)return ahfb; << 1286 } 575 } 1287 return mediumDiffuseness[A-1]; << 1288 } else if(A < 6 && A >= 2) { << 1289 INCL_ERROR("getSurfaceDiffuseness: wa << 1290 return 0.0; << 1291 } else { 576 } else { 1292 INCL_ERROR("getSurfaceDiffuseness: No << 577 ERROR("ParticleTable::getRadiusParameter : No radius for nucleus A = " << A << " Z = " << Z << std::endl); 1293 return 0.0; 578 return 0.0; 1294 } 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; 1295 } 600 } >> 601 } 1296 602 1297 G4double getMomentumRMS(const G4int A, co << 603 #ifdef INCLXX_IN_GEANT4_MODE 1298 // assert(Z>=0 && A>=0 && Z<=A); << 604 G4double ParticleTable::getSurfaceDiffuseness(const G4int A, const G4int /*Z*/ ) { 1299 return getFermiMomentum(A,Z) * Math::sq << 605 #else 1300 } << 606 G4double ParticleTable::getSurfaceDiffuseness(const G4int A, const G4int Z) { >> 607 #endif // INCLXX_IN_GEANT4_MODE 1301 608 1302 G4double getSeparationEnergyINCL(const Pa << 609 if(A >= 28) { 1303 if(t==Proton) << 610 return 1.63e-4 * A + 0.510; 1304 return theINCLProtonSeparationEnergy; << 611 } else if(A < 28 && A >= 19) { 1305 else if(t==Neutron) << 612 return mediumDiffuseness[A-1]; 1306 return theINCLNeutronSeparationEnergy << 613 } else if(A < 19 && A >= 6) { 1307 else if(t==Lambda) << 614 return mediumDiffuseness[A-1]; 1308 return theINCLLambdaSeparationEnergy; << 615 } else if(A < 6 && A >= 2) { 1309 else if(t==antiProton) << 616 ERROR("ParticleTable::getSurfaceDiffuseness: was called for A = " << A << " Z = " << Z << std::endl); 1310 return theINCLantiProtonSeparationEne << 617 return 0.0; 1311 else { << 618 } else { 1312 INCL_ERROR("ParticleTable::getSeparat << 619 ERROR("ParticleTable::getSurfaceDiffuseness: No diffuseness for nucleus A = " << A << " Z = " << Z << std::endl); 1313 return 0.0; << 620 return 0.0; 1314 } << 621 } 1315 } << 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 } 1316 633 1317 G4double getSeparationEnergyReal(const Pa << 634 #ifndef INCLXX_IN_GEANT4_MODE 1318 // Real separation energies for all nuc << 635 void ParticleTable::readRealMasses(std::string const &path) { 1319 if(t==Proton) << 636 // Clear the existing tables, if any 1320 return (*getTableParticleMass)(Proton << 637 massTableMask.clear(); 1321 else if(t==Neutron) << 638 massTable.clear(); 1322 return (*getTableParticleMass)(Neutro << 639 1323 else if(t==Lambda) << 640 // File name 1324 return (*getTableParticleMass)(Lambda << 641 std::string fileName(path + "/walletlifetime.dat"); 1325 else { << 642 DEBUG("Reading real nuclear masses from file " << fileName << std::endl); 1326 INCL_ERROR("ParticleTable::getSeparat << 643 1327 return 0.0; << 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); 1328 } 666 } 1329 } << 1330 << 1331 G4double getSeparationEnergyRealForLight( << 1332 // Real separation energies for light n << 1333 if(Z<clusterTableZSize && A<clusterTabl << 1334 return getSeparationEnergyReal(t, A, << 1335 else << 1336 return getSeparationEnergyINCL(t, A, << 1337 } << 1338 << 1339 G4double getProtonSeparationEnergy() { re << 1340 << 1341 G4double getNeutronSeparationEnergy() { r << 1342 << 1343 G4double getLambdaSeparationEnergy() { re << 1344 << 1345 void setProtonSeparationEnergy(const G4do << 1346 << 1347 void setNeutronSeparationEnergy(const G4d << 1348 << 1349 void setLambdaSeparationEnergy(const G4do << 1350 << 1351 std::string getElementName(const G4int Z) << 1352 if(Z<1) { << 1353 INCL_WARN("getElementName called with << 1354 return elementTable[0]; << 1355 } else if(Z<elementTableSize) << 1356 return elementTable[Z]; << 1357 else << 1358 return getIUPACElementName(Z); << 1359 } << 1360 << 1361 std::string getIUPACElementName(const G4i << 1362 std::stringstream elementStream; << 1363 elementStream << Z; << 1364 std::string elementName = elementStream << 1365 std::transform(elementName.begin(), ele << 1366 elementName[0] = (char)std::toupper(ele << 1367 return elementName; << 1368 } << 1369 << 1370 G4int parseElement(std::string pS) { << 1371 // Normalize the element name << 1372 std::transform(pS.begin(), pS.end(), pS << 1373 pS[0] = (char)std::toupper(pS[0]); << 1374 << 1375 const std::string *iter = std::find(ele << 1376 if(iter != elementTable+elementTableSiz << 1377 return G4int(iter - elementTable); << 1378 else << 1379 return ParticleTable::parseIUPACEleme << 1380 } << 1381 << 1382 G4int parseIUPACElement(std::string const << 1383 // Normalise to lower case << 1384 std::string elementName(sel); << 1385 std::transform(elementName.begin(), ele << 1386 // Return 0 if the element name contain << 1387 if(elementName.find_first_not_of(elemen << 1388 return 0; << 1389 std::transform(elementName.begin(), ele << 1390 std::stringstream elementStream(element << 1391 G4int Z; << 1392 elementStream >> Z; << 1393 return Z; << 1394 } << 1395 << 1396 IsotopicDistribution const &getNaturalIso << 1397 return getNaturalIsotopicDistributions( << 1398 } << 1399 << 1400 G4int drawRandomNaturalIsotope(const G4in << 1401 return getNaturalIsotopicDistributions( << 1402 } << 1403 << 1404 G4double getFermiMomentumConstant(const G << 1405 return constantFermiMomentum; << 1406 } << 1407 << 1408 G4double getFermiMomentumConstantLight(co << 1409 // assert(Z>0 && A>0 && Z<=A); << 1410 if(Z<clusterTableZSize && A<clusterTabl << 1411 const G4double rms = momentumRMS[Z][A << 1412 return ((rms>0.) ? rms : momentumRMS[ << 1413 } else << 1414 return getFermiMomentumConstant(A,Z); << 1415 } << 1416 << 1417 G4double getFermiMomentumMassDependent(co << 1418 // assert(A>0); << 1419 static const G4double alphaParam = 259. << 1420 static const G4double betaParam = 152. << 1421 static const G4double gammaParam = 9.51 << 1422 return alphaParam - betaParam*std::exp( << 1423 } << 1424 << 1425 G4double getRPCorrelationCoefficient(cons << 1426 // assert(t==Proton || t==Neutron || t==Lambd << 1427 return rpCorrelationCoefficient[t]; << 1428 } << 1429 << 1430 G4double getNeutronSkin() { return neutro << 1431 667 1432 G4double getNeutronHalo() { return neutro << 668 massTable.at(Z).at(A) = A*amu + excess - Z*eMass; >> 669 massTableMask.at(Z).at(A) = true; 1433 670 1434 G4ThreadLocal G4double minDeltaMass = 0.; << 671 massTableIn >> A >> Z >> excess; 1435 G4ThreadLocal G4double minDeltaMass2 = 0. << 672 } while(massTableIn.good()); 1436 G4ThreadLocal G4double minDeltaMassRndm = << 1437 G4ThreadLocal NuclearMassFn getTableMass << 1438 G4ThreadLocal ParticleMassFn getTablePart << 1439 G4ThreadLocal SeparationEnergyFn getSepar << 1440 G4ThreadLocal FermiMomentumFn getFermiMom << 1441 << 1442 ParticleType getPionType(const G4int isos << 1443 // assert(isosp == -2 || isosp == 0 || isosp << 1444 if (isosp == -2) { << 1445 return PiMinus; << 1446 } << 1447 else if (isosp == 0) { << 1448 return PiZero; << 1449 } << 1450 else { << 1451 return PiPlus; << 1452 } << 1453 } << 1454 673 1455 ParticleType getNucleonType(const G4int i << 674 } 1456 // assert(isosp == -1 || isosp == 1); << 675 #endif 1457 if (isosp == -1) { << 1458 return Neutron; << 1459 } << 1460 else { << 1461 return Proton; << 1462 } << 1463 } << 1464 << 1465 ParticleType getDeltaType(const G4int iso << 1466 // assert(isosp == -3 || isosp == -1 || isosp << 1467 if (isosp == -3) { << 1468 return DeltaMinus; << 1469 } << 1470 else if (isosp == -1) { << 1471 return DeltaZero; << 1472 } << 1473 else if (isosp == 1) { << 1474 return DeltaPlus; << 1475 } << 1476 else { << 1477 return DeltaPlusPlus; << 1478 } << 1479 } << 1480 << 1481 ParticleType getSigmaType(const G4int iso << 1482 // assert(isosp == -2 || isosp == 0 || isosp << 1483 if (isosp == -2) { << 1484 return SigmaMinus; << 1485 } << 1486 else if (isosp == 0) { << 1487 return SigmaZero; << 1488 } << 1489 else { << 1490 return SigmaPlus; << 1491 } << 1492 } << 1493 << 1494 ParticleType getXiType(const G4int isosp) << 1495 // assert(isosp == -1 || isosp == 1); << 1496 if (isosp == -1) { << 1497 return XiMinus; << 1498 } << 1499 else { << 1500 return XiZero; << 1501 } << 1502 } << 1503 << 1504 /*ParticleType getAntiNucleonType(const G4i << 1505 // assert(isosp == -1); //|| isosp == 1 << 1506 if (isosp == -1) { << 1507 return antiProton; << 1508 } << 1509 else { << 1510 return antiNeutron; << 1511 } << 1512 }*/ << 1513 << 1514 ParticleType getAntiSigmaType(const G4int << 1515 // assert(isosp == -2 || isosp == 0 || isosp << 1516 if (isosp == -2) { << 1517 return antiSigmaPlus; << 1518 } << 1519 else if (isosp == 0) { << 1520 return antiSigmaZero; << 1521 } << 1522 else { << 1523 return antiSigmaMinus; << 1524 } << 1525 } << 1526 << 1527 ParticleType getAntiXiType(const G4int is << 1528 // assert(isosp == -1 || isosp == 1); << 1529 if (isosp == -1) { << 1530 return antiXiZero; << 1531 } << 1532 else { << 1533 return antiXiMinus; << 1534 } << 1535 } << 1536 << 1537 ParticleType getKaonType(const G4int isos << 1538 // assert(isosp == -1 || isosp == 1); << 1539 if (isosp == -1) { << 1540 return KZero; << 1541 } << 1542 else { << 1543 return KPlus; << 1544 } << 1545 } << 1546 << 1547 ParticleType getAntiKaonType(const G4int << 1548 // assert(isosp == -1 || isosp == 1); << 1549 if (isosp == -1) { << 1550 return KMinus; << 1551 } << 1552 else { << 1553 return KZeroBar; << 1554 } << 1555 } << 1556 676 1557 G4double getWidth(const ParticleType pt) << 677 std::string ParticleTable::getIUPACElementName(const G4int Z) { 1558 // assert(pt == PiPlus || pt == PiMinus || pt << 678 std::stringstream elementStream; 1559 if(pt == PiPlus) { << 679 elementStream << Z; 1560 return piPlusWidth; << 680 std::string elementName = elementStream.str(); 1561 } else if(pt == PiMinus) { << 681 std::transform(elementName.begin(), elementName.end(), elementName.begin(), ParticleTable::intToIUPAC); 1562 return piMinusWidth; << 682 elementName[0] = std::toupper(elementName.at(0)); 1563 } else if(pt == PiZero) { << 683 return elementName; 1564 return piZeroWidth; << 684 } 1565 } else if(pt == Eta) { << 685 1566 return etaWidth; << 686 G4int ParticleTable::parseIUPACElement(std::string const &s) { 1567 } else if(pt == Omega) { << 687 // Normalise to lower case 1568 return omegaWidth; << 688 std::string elementName(s); 1569 } else if(pt == EtaPrime) { << 689 std::transform(elementName.begin(), elementName.end(), elementName.begin(), ::tolower); 1570 return etaPrimeWidth; << 690 // Return 0 if the element name contains anything but IUPAC digits 1571 } else if(pt == SigmaPlus) { << 691 if(elementName.find_first_not_of(elementIUPACDigits)!=std::string::npos) 1572 return SigmaPlusWidth; << 692 return 0; 1573 } else if(pt == SigmaZero) { << 693 std::transform(elementName.begin(), elementName.end(), elementName.begin(), ParticleTable::iupacToInt); 1574 return SigmaZeroWidth; << 694 std::stringstream elementStream(elementName); 1575 } else if(pt == SigmaMinus) { << 695 G4int Z; 1576 return SigmaMinusWidth; << 696 elementStream >> Z; 1577 } else if(pt == KPlus) { << 697 return Z; 1578 return KPlusWidth; << 698 } 1579 } else if(pt == KMinus) { << 1580 return KMinusWidth; << 1581 } else if(pt == KShort) { << 1582 return KShortWidth; << 1583 } else if(pt == KLong) { << 1584 return KLongWidth; << 1585 } else if(pt == Lambda) { << 1586 return LambdaWidth; << 1587 } else if(pt == XiMinus) { << 1588 return XiMinusWidth; << 1589 } else if(pt == XiZero) { << 1590 return XiZeroWidth; << 1591 } else if(pt == antiSigmaPlus) { << 1592 return antiSigmaPlusWidth; << 1593 } else if(pt == antiSigmaZero) { << 1594 return antiSigmaZeroWidth; << 1595 } else if(pt == antiSigmaMinus) { << 1596 return antiSigmaMinusWidth; << 1597 } else if(pt == antiLambda) { << 1598 return antiLambdaWidth; << 1599 } else if(pt == antiXiMinus) { << 1600 return antiXiMinusWidth; << 1601 } else if(pt == antiXiZero) { << 1602 return antiXiZeroWidth; << 1603 } else { << 1604 INCL_ERROR("getWidth : Unknown << 1605 return 0.0; << 1606 } << 1607 } << 1608 << 1609 } // namespace ParticleTable << 1610 } // namespace G4INCL << 1611 699 >> 700 } 1612 701