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 // Alain Boudard, CEA-Saclay, France 28 // Joseph Cugnon, University of Liege, Belgium 28 // Joseph Cugnon, University of Liege, Belgium 29 // Jean-Christophe David, CEA-Saclay, France 29 // Jean-Christophe David, CEA-Saclay, France 30 // Pekka Kaitaniemi, CEA-Saclay, France, and H 30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland 31 // Sylvie Leray, CEA-Saclay, France 31 // Sylvie Leray, CEA-Saclay, France 32 // Davide Mancusi, CEA-Saclay, France 32 // Davide Mancusi, CEA-Saclay, France 33 // 33 // 34 #define INCLXX_IN_GEANT4_MODE 1 34 #define INCLXX_IN_GEANT4_MODE 1 35 35 36 #include "globals.hh" 36 #include "globals.hh" 37 37 38 #include "G4INCLParticleTable.hh" 38 #include "G4INCLParticleTable.hh" 39 #include "G4INCLNuclearMassTable.hh" 39 #include "G4INCLNuclearMassTable.hh" 40 #include <algorithm> 40 #include <algorithm> 41 // #include <cassert> 41 // #include <cassert> 42 #include <cmath> 42 #include <cmath> 43 #include <cctype> 43 #include <cctype> 44 #include <sstream> 44 #include <sstream> 45 #ifdef INCLXX_IN_GEANT4_MODE 45 #ifdef INCLXX_IN_GEANT4_MODE 46 #include "G4SystemOfUnits.hh" 46 #include "G4SystemOfUnits.hh" 47 #endif 47 #endif 48 48 49 #ifdef INCLXX_IN_GEANT4_MODE 49 #ifdef INCLXX_IN_GEANT4_MODE 50 #include "G4PhysicalConstants.hh" 50 #include "G4PhysicalConstants.hh" 51 #include "G4SystemOfUnits.hh" 51 #include "G4SystemOfUnits.hh" 52 #endif 52 #endif 53 53 54 namespace G4INCL { 54 namespace G4INCL { 55 55 56 namespace ParticleTable { 56 namespace ParticleTable { 57 57 58 namespace { 58 namespace { 59 59 60 /// \brief Static instance of the Natura 60 /// \brief Static instance of the NaturalIsotopicAbundances class 61 const NaturalIsotopicDistributions *theN 61 const NaturalIsotopicDistributions *theNaturalIsotopicDistributions = NULL; 62 62 63 const G4double theINCLNucleonMass = 938. 63 const G4double theINCLNucleonMass = 938.2796; 64 const G4double theINCLPionMass = 138.0; 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; 65 const G4double theINCLEtaMass = 547.862; 69 const G4double theINCLOmegaMass = 782.65 66 const G4double theINCLOmegaMass = 782.65; 70 const G4double theINCLEtaPrimeMass = 957 67 const G4double theINCLEtaPrimeMass = 957.78; 71 const G4double theINCLPhotonMass = 0.0; 68 const G4double theINCLPhotonMass = 0.0; 72 G4ThreadLocal G4double protonMass = 0.0; 69 G4ThreadLocal G4double protonMass = 0.0; 73 G4ThreadLocal G4double neutronMass = 0.0 70 G4ThreadLocal G4double neutronMass = 0.0; 74 G4ThreadLocal G4double piPlusMass = 0.0; 71 G4ThreadLocal G4double piPlusMass = 0.0; 75 G4ThreadLocal G4double piMinusMass = 0.0 72 G4ThreadLocal G4double piMinusMass = 0.0; 76 G4ThreadLocal G4double piZeroMass = 0.0; 73 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; 74 G4ThreadLocal G4double etaMass = 0.0; 98 G4ThreadLocal G4double omegaMass = 0.0; 75 G4ThreadLocal G4double omegaMass = 0.0; 99 G4ThreadLocal G4double etaPrimeMass = 0. 76 G4ThreadLocal G4double etaPrimeMass = 0.0; 100 G4ThreadLocal G4double photonMass = 0.0; 77 G4ThreadLocal G4double photonMass = 0.0; 101 78 102 // Hard-coded values of the real particl 79 // Hard-coded values of the real particle masses (MeV/c^2) 103 G4ThreadLocal G4double theRealProtonMass 80 G4ThreadLocal G4double theRealProtonMass = 938.27203; 104 G4ThreadLocal G4double theRealNeutronMas 81 G4ThreadLocal G4double theRealNeutronMass = 939.56536; 105 G4ThreadLocal G4double theRealChargedPiM 82 G4ThreadLocal G4double theRealChargedPiMass = 139.57018; 106 G4ThreadLocal G4double theRealPiZeroMass 83 G4ThreadLocal G4double theRealPiZeroMass = 134.9766; 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 = 84 G4ThreadLocal G4double theRealEtaMass = 547.862; 124 G4ThreadLocal G4double theRealOmegaMass 85 G4ThreadLocal G4double theRealOmegaMass = 782.65; 125 G4ThreadLocal G4double theRealEtaPrimeMa 86 G4ThreadLocal G4double theRealEtaPrimeMass = 957.78; 126 G4ThreadLocal G4double theRealPhotonMass 87 G4ThreadLocal G4double theRealPhotonMass = 0.0; 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 88 >> 89 // Width (second) >> 90 const G4double theChargedPiWidth = 2.6033e-08; >> 91 const G4double thePiZeroWidth = 8.52e-17; >> 92 const G4double theEtaWidth = 5.025e-19; // 1.31 keV >> 93 const G4double theOmegaWidth = 7.7528e-23; // 8.49 MeV >> 94 const G4double theEtaPrimeWidth = 3.3243e-21; // 0.198 MeV >> 95 G4ThreadLocal G4double piPlusWidth = 0.0; >> 96 G4ThreadLocal G4double piMinusWidth = 0.0; >> 97 G4ThreadLocal G4double piZeroWidth = 0.0; >> 98 G4ThreadLocal G4double etaWidth = 0.0; >> 99 G4ThreadLocal G4double omegaWidth = 0.0; >> 100 G4ThreadLocal G4double etaPrimeWidth = 0.0; >> 101 >> 102 172 const G4int mediumNucleiTableSize = 30; 103 const G4int mediumNucleiTableSize = 30; 173 104 174 const G4double mediumDiffuseness[mediumN 105 const G4double mediumDiffuseness[mediumNucleiTableSize] = 175 {0.0,0.0,0.0,0.0,0.0,1.78,1.77,1.77,1.69 << 106 {0.0,0.0,0.0,0.0,0.0,1.78,1.77,1.77,1.77,1.71, 176 1.69,1.72,1.635,1.730,1.81,1.833,1.798 << 107 1.69,1.69,1.635,1.730,1.81,1.833,1.798, 177 1.93,0.567,0.571, 0.560,0.549,0.550,0. << 108 1.841,0.567,0.571, 0.560,0.549,0.550,0.551, 178 0.580,0.575,0.569,0.537,0.0,0.0}; 109 0.580,0.575,0.569,0.537,0.0,0.0}; 179 const G4double mediumRadius[mediumNuclei 110 const G4double mediumRadius[mediumNucleiTableSize] = 180 {0.0,0.0,0.0,0.0,0.0,0.334,0.327,0.479,0 111 {0.0,0.0,0.0,0.0,0.0,0.334,0.327,0.479,0.631,0.838, 181 0.811,0.84,1.403,1.335,1.25,1.544,1.49 << 112 0.811,1.07,1.403,1.335,1.25,1.544,1.498,1.513, 182 2.58,2.77, 2.775,2.78,2.88,2.98,3.22,3 113 2.58,2.77, 2.775,2.78,2.88,2.98,3.22,3.03,2.84, 183 3.14,0.0,0.0}; 114 3.14,0.0,0.0}; 184 115 185 const G4double positionRMS[clusterTableZ 116 const G4double positionRMS[clusterTableZSize][clusterTableASize] = { 186 /* A= 0 1 2 3 4 117 /* A= 0 1 2 3 4 5 6 7 8 9 10 11 12 */ 187 /* Z=0 */ {-1.0, -1.0, -1.0, -1.0, -1. 118 /* Z=0 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0}, 188 /* Z=1 */ {-1.0, -1.0, 2.10, 1.80, 1.7 119 /* Z=1 */ {-1.0, -1.0, 2.10, 1.80, 1.70, 1.83, 2.60, 2.50, -1.0, -1.0, -1.0, -1.0, -1.0}, 189 /* Z=2 */ {-1.0, -1.0, -1.0, 1.80, 1.6 120 /* Z=2 */ {-1.0, -1.0, -1.0, 1.80, 1.68, 1.70, 2.60, 2.50, 2.50, 2.50, 2.50, -1.0, -1.0}, 190 /* Z=3 */ {-1.0, -1.0, -1.0, -1.0, 1.7 121 /* Z=3 */ {-1.0, -1.0, -1.0, -1.0, 1.70, 1.83, 2.56, 2.40, 2.50, 2.50, 2.50, 2.50, 2.50}, 191 /* Z=4 */ {-1.0, -1.0, -1.0, -1.0, -1. 122 /* Z=4 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 2.60, 2.50, 2.50, 2.51, 2.50, 2.50, 2.50}, 192 /* Z=5 */ {-1.0, -1.0, -1.0, -1.0, -1. 123 /* Z=5 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 2.50, 2.50, 2.50, 2.50, 2.45, 2.40, 2.50}, 193 /* Z=6 */ {-1.0, -1.0, -1.0, -1.0, -1. 124 /* Z=6 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 2.50, 2.50, 2.50, 2.50, 2.47}, 194 /* Z=7 */ {-1.0, -1.0, -1.0, -1.0, -1. 125 /* Z=7 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 2.50, 2.50, 2.50}, 195 /* Z=8 */ {-1.0, -1.0, -1.0, -1.0, -1. 126 /* Z=8 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 2.50} 196 }; 127 }; 197 128 198 const G4double momentumRMS[clusterTableZ 129 const G4double momentumRMS[clusterTableZSize][clusterTableASize] = { 199 /* A= 0 1 2 3 4 130 /* A= 0 1 2 3 4 5 6 7 8 9 10 11 12 */ 200 /* Z=0 */ {-1.0, -1.0, -1.0, -1.0, -1. 131 /* Z=0 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0}, 201 /* Z=1 */ {-1.0, -1.0, 77.0, 110., 153 132 /* Z=1 */ {-1.0, -1.0, 77.0, 110., 153., 100., 100., 100., -1.0, -1.0, -1.0, -1.0, -1.0}, 202 /* Z=2 */ {-1.0, -1.0, -1.0, 110., 153 133 /* Z=2 */ {-1.0, -1.0, -1.0, 110., 153., 100., 100., 100., 100., 100., 100., -1.0, -1.0}, 203 /* Z=3 */ {-1.0, -1.0, -1.0, -1.0, 153 134 /* Z=3 */ {-1.0, -1.0, -1.0, -1.0, 153., 100., 100., 100., 100., 100., 100., 100., 100.}, 204 /* Z=4 */ {-1.0, -1.0, -1.0, -1.0, -1. 135 /* Z=4 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 100., 100., 100., 100., 100., 100., 100.}, 205 /* Z=5 */ {-1.0, -1.0, -1.0, -1.0, -1. 136 /* Z=5 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 100., 100., 100., 100., 100., 100., 100.}, 206 /* Z=6 */ {-1.0, -1.0, -1.0, -1.0, -1. 137 /* Z=6 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 100., 100., 100., 100., 100.}, 207 /* Z=7 */ {-1.0, -1.0, -1.0, -1.0, -1. 138 /* Z=7 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 100., 100., 100.}, 208 /* Z=8 */ {-1.0, -1.0, -1.0, -1.0, -1. 139 /* Z=8 */ {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 100.} 209 }; 140 }; 210 141 211 const G4int elementTableSize = 113; // u 142 const G4int elementTableSize = 113; // up to Cn 212 143 213 /// \brief Table of chemical element nam 144 /// \brief Table of chemical element names 214 const std::string elementTable[elementTa 145 const std::string elementTable[elementTableSize] = { 215 "", 146 "", 216 "H", 147 "H", 217 "He", 148 "He", 218 "Li", 149 "Li", 219 "Be", 150 "Be", 220 "B", 151 "B", 221 "C", 152 "C", 222 "N", 153 "N", 223 "O", 154 "O", 224 "F", 155 "F", 225 "Ne", 156 "Ne", 226 "Na", 157 "Na", 227 "Mg", 158 "Mg", 228 "Al", 159 "Al", 229 "Si", 160 "Si", 230 "P", 161 "P", 231 "S", 162 "S", 232 "Cl", 163 "Cl", 233 "Ar", 164 "Ar", 234 "K", 165 "K", 235 "Ca", 166 "Ca", 236 "Sc", 167 "Sc", 237 "Ti", 168 "Ti", 238 "V", 169 "V", 239 "Cr", 170 "Cr", 240 "Mn", 171 "Mn", 241 "Fe", 172 "Fe", 242 "Co", 173 "Co", 243 "Ni", 174 "Ni", 244 "Cu", 175 "Cu", 245 "Zn", 176 "Zn", 246 "Ga", 177 "Ga", 247 "Ge", 178 "Ge", 248 "As", 179 "As", 249 "Se", 180 "Se", 250 "Br", 181 "Br", 251 "Kr", 182 "Kr", 252 "Rb", 183 "Rb", 253 "Sr", 184 "Sr", 254 "Y", 185 "Y", 255 "Zr", 186 "Zr", 256 "Nb", 187 "Nb", 257 "Mo", 188 "Mo", 258 "Tc", 189 "Tc", 259 "Ru", 190 "Ru", 260 "Rh", 191 "Rh", 261 "Pd", 192 "Pd", 262 "Ag", 193 "Ag", 263 "Cd", 194 "Cd", 264 "In", 195 "In", 265 "Sn", 196 "Sn", 266 "Sb", 197 "Sb", 267 "Te", 198 "Te", 268 "I", 199 "I", 269 "Xe", 200 "Xe", 270 "Cs", 201 "Cs", 271 "Ba", 202 "Ba", 272 "La", 203 "La", 273 "Ce", 204 "Ce", 274 "Pr", 205 "Pr", 275 "Nd", 206 "Nd", 276 "Pm", 207 "Pm", 277 "Sm", 208 "Sm", 278 "Eu", 209 "Eu", 279 "Gd", 210 "Gd", 280 "Tb", 211 "Tb", 281 "Dy", 212 "Dy", 282 "Ho", 213 "Ho", 283 "Er", 214 "Er", 284 "Tm", 215 "Tm", 285 "Yb", 216 "Yb", 286 "Lu", 217 "Lu", 287 "Hf", 218 "Hf", 288 "Ta", 219 "Ta", 289 "W", 220 "W", 290 "Re", 221 "Re", 291 "Os", 222 "Os", 292 "Ir", 223 "Ir", 293 "Pt", 224 "Pt", 294 "Au", 225 "Au", 295 "Hg", 226 "Hg", 296 "Tl", 227 "Tl", 297 "Pb", 228 "Pb", 298 "Bi", 229 "Bi", 299 "Po", 230 "Po", 300 "At", 231 "At", 301 "Rn", 232 "Rn", 302 "Fr", 233 "Fr", 303 "Ra", 234 "Ra", 304 "Ac", 235 "Ac", 305 "Th", 236 "Th", 306 "Pa", 237 "Pa", 307 "U", 238 "U", 308 "Np", 239 "Np", 309 "Pu", 240 "Pu", 310 "Am", 241 "Am", 311 "Cm", 242 "Cm", 312 "Bk", 243 "Bk", 313 "Cf", 244 "Cf", 314 "Es", 245 "Es", 315 "Fm", 246 "Fm", 316 "Md", 247 "Md", 317 "No", 248 "No", 318 "Lr", 249 "Lr", 319 "Rf", 250 "Rf", 320 "Db", 251 "Db", 321 "Sg", 252 "Sg", 322 "Bh", 253 "Bh", 323 "Hs", 254 "Hs", 324 "Mt", 255 "Mt", 325 "Ds", 256 "Ds", 326 "Rg", 257 "Rg", 327 "Cn" 258 "Cn" 328 }; 259 }; 329 260 330 /// \brief Digit names to compose IUPAC 261 /// \brief Digit names to compose IUPAC element names 331 const std::string elementIUPACDigits = " 262 const std::string elementIUPACDigits = "nubtqphsoe"; 332 263 333 #define INCL_DEFAULT_SEPARATION_ENERGY 6.83 264 #define INCL_DEFAULT_SEPARATION_ENERGY 6.83 334 const G4double theINCLProtonSeparationEn 265 const G4double theINCLProtonSeparationEnergy = INCL_DEFAULT_SEPARATION_ENERGY; 335 const G4double theINCLNeutronSeparationE 266 const G4double theINCLNeutronSeparationEnergy = INCL_DEFAULT_SEPARATION_ENERGY; 336 const G4double theINCLLambdaSeparationEn << 337 //const G4double theINCLantiProtonSepara << 338 const G4double theINCLantiProtonSeparati << 339 G4ThreadLocal G4double protonSeparationE 267 G4ThreadLocal G4double protonSeparationEnergy = INCL_DEFAULT_SEPARATION_ENERGY; 340 G4ThreadLocal G4double neutronSeparation 268 G4ThreadLocal G4double neutronSeparationEnergy = INCL_DEFAULT_SEPARATION_ENERGY; 341 G4ThreadLocal G4double lambdaSeparationE << 342 //G4ThreadLocal G4double antiprotonSepar << 343 //G4ThreadLocal G4double antiprotonSepar << 344 #undef INCL_DEFAULT_SEPARATION_ENERGY 269 #undef INCL_DEFAULT_SEPARATION_ENERGY 345 270 346 G4ThreadLocal G4double rpCorrelationCoef 271 G4ThreadLocal G4double rpCorrelationCoefficient[UnknownParticle]; 347 272 348 G4ThreadLocal G4double neutronSkin = 0.0 273 G4ThreadLocal G4double neutronSkin = 0.0; 349 G4ThreadLocal G4double neutronHalo = 0.0 274 G4ThreadLocal G4double neutronHalo = 0.0; 350 275 351 #ifdef INCLXX_IN_GEANT4_MODE 276 #ifdef INCLXX_IN_GEANT4_MODE 352 G4ThreadLocal G4IonTable *theG4IonTable; 277 G4ThreadLocal G4IonTable *theG4IonTable; 353 #endif 278 #endif 354 279 355 /// \brief Default value for constant Fe 280 /// \brief Default value for constant Fermi momentum 356 G4ThreadLocal G4double constantFermiMome 281 G4ThreadLocal G4double constantFermiMomentum = 0.0; 357 282 358 /// \brief Transform a IUPAC char to an 283 /// \brief Transform a IUPAC char to an char representing an integer digit 359 char iupacToInt(char c) { 284 char iupacToInt(char c) { 360 return (char)(((G4int)'0')+elementIUPA 285 return (char)(((G4int)'0')+elementIUPACDigits.find(c)); 361 } 286 } 362 287 363 /// \brief Transform an integer digit (r 288 /// \brief Transform an integer digit (represented by a char) to a IUPAC char 364 char intToIUPAC(char n) { return element 289 char intToIUPAC(char n) { return elementIUPACDigits.at(n); } 365 290 366 /// \brief Get the singleton instance of 291 /// \brief Get the singleton instance of the natural isotopic distributions 367 const NaturalIsotopicDistributions *getN 292 const NaturalIsotopicDistributions *getNaturalIsotopicDistributions() { 368 if(!theNaturalIsotopicDistributions) 293 if(!theNaturalIsotopicDistributions) 369 theNaturalIsotopicDistributions = ne 294 theNaturalIsotopicDistributions = new NaturalIsotopicDistributions; 370 return theNaturalIsotopicDistributions 295 return theNaturalIsotopicDistributions; 371 } 296 } 372 297 373 } // namespace 298 } // namespace 374 299 375 void initialize(Config const * const theCo 300 void initialize(Config const * const theConfig /*=0*/) { 376 protonMass = theINCLNucleonMass; 301 protonMass = theINCLNucleonMass; 377 neutronMass = theINCLNucleonMass; << 302 neutronMass = theINCLNucleonMass; 378 piPlusMass = theINCLPionMass; 303 piPlusMass = theINCLPionMass; 379 piMinusMass = theINCLPionMass; 304 piMinusMass = theINCLPionMass; 380 piZeroMass = theINCLPionMass; 305 piZeroMass = theINCLPionMass; 381 << 382 etaMass = theINCLEtaMass; 306 etaMass = theINCLEtaMass; 383 omegaMass = theINCLOmegaMass; 307 omegaMass = theINCLOmegaMass; 384 etaPrimeMass = theINCLEtaPrimeMass; 308 etaPrimeMass = theINCLEtaPrimeMass; 385 photonMass = theINCLPhotonMass; 309 photonMass = theINCLPhotonMass; 386 310 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 311 if(theConfig && theConfig->getUseRealMasses()) { 410 getTableMass = getRealMass; 312 getTableMass = getRealMass; 411 getTableParticleMass = getRealMass; 313 getTableParticleMass = getRealMass; 412 } else { 314 } else { 413 getTableMass = getINCLMass; 315 getTableMass = getINCLMass; 414 getTableParticleMass = getINCLMass; 316 getTableParticleMass = getINCLMass; 415 } 317 } 416 318 417 #ifndef INCLXX_IN_GEANT4_MODE 319 #ifndef INCLXX_IN_GEANT4_MODE 418 std::string dataFilePath; 320 std::string dataFilePath; 419 if(theConfig) 321 if(theConfig) 420 dataFilePath = theConfig->getINCLXXDat 322 dataFilePath = theConfig->getINCLXXDataFilePath(); 421 NuclearMassTable::initialize(dataFilePat 323 NuclearMassTable::initialize(dataFilePath, getRealMass(Proton), getRealMass(Neutron)); 422 #endif 324 #endif 423 325 424 #ifdef INCLXX_IN_GEANT4_MODE 326 #ifdef INCLXX_IN_GEANT4_MODE 425 G4ParticleTable *theG4ParticleTable = G4 327 G4ParticleTable *theG4ParticleTable = G4ParticleTable::GetParticleTable(); 426 theG4IonTable = theG4ParticleTable->GetI 328 theG4IonTable = theG4ParticleTable->GetIonTable(); 427 theRealProtonMass = theG4ParticleTable-> 329 theRealProtonMass = theG4ParticleTable->FindParticle("proton")->GetPDGMass() / MeV; 428 theRealNeutronMass = theG4ParticleTable- << 330 theRealNeutronMass = theG4ParticleTable->FindParticle("neutron")->GetPDGMass() / MeV; 429 theRealChargedPiMass = theG4ParticleTabl 331 theRealChargedPiMass = theG4ParticleTable->FindParticle("pi+")->GetPDGMass() / MeV; 430 theRealPiZeroMass = theG4ParticleTable-> 332 theRealPiZeroMass = theG4ParticleTable->FindParticle("pi0")->GetPDGMass() / MeV; 431 << 432 theRealEtaMass = theG4ParticleTable->Fin 333 theRealEtaMass = theG4ParticleTable->FindParticle("eta")->GetPDGMass() / MeV; 433 theRealOmegaMass = theG4ParticleTable->F 334 theRealOmegaMass = theG4ParticleTable->FindParticle("omega")->GetPDGMass() / MeV; 434 theRealEtaPrimeMass = theG4ParticleTable 335 theRealEtaPrimeMass = theG4ParticleTable->FindParticle("eta_prime")->GetPDGMass() / MeV; 435 theRealPhotonMass = theG4ParticleTable-> 336 theRealPhotonMass = theG4ParticleTable->FindParticle("gamma")->GetPDGMass() / MeV; 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 337 #endif 455 338 456 minDeltaMass = theRealNeutronMass + theR 339 minDeltaMass = theRealNeutronMass + theRealChargedPiMass + 0.5; 457 minDeltaMass2 = minDeltaMass*minDeltaMas 340 minDeltaMass2 = minDeltaMass*minDeltaMass; 458 minDeltaMassRndm = std::atan((minDeltaMa 341 minDeltaMassRndm = std::atan((minDeltaMass-effectiveDeltaMass)*2./effectiveDeltaWidth); 459 342 460 piPlusWidth = theChargedPiWidth; << 343 piPlusWidth = theChargedPiWidth; 461 piMinusWidth = theChargedPiWidth; << 344 piMinusWidth = theChargedPiWidth; 462 piZeroWidth = thePiZeroWidth; << 345 piZeroWidth = thePiZeroWidth; 463 etaWidth = theEtaWidth; << 346 etaWidth = theEtaWidth; 464 omegaWidth = theOmegaWidth; << 347 omegaWidth = theOmegaWidth; 465 etaPrimeWidth = theEtaPrimeWidth; << 348 etaPrimeWidth = theEtaPrimeWidth; 466 << 349 467 SigmaMinusWidth = theSigmaMinusWidth; << 350 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 << 478 HFB::initialize(); << 479 #else << 480 HFB::initialize(dataFilePath); << 481 #endif << 482 << 483 // Initialise the separation-energy func 351 // Initialise the separation-energy function 484 if(!theConfig || theConfig->getSeparatio 352 if(!theConfig || theConfig->getSeparationEnergyType()==INCLSeparationEnergy) 485 getSeparationEnergy = getSeparationEne 353 getSeparationEnergy = getSeparationEnergyINCL; 486 else if(theConfig->getSeparationEnergyTy 354 else if(theConfig->getSeparationEnergyType()==RealSeparationEnergy) 487 getSeparationEnergy = getSeparationEne 355 getSeparationEnergy = getSeparationEnergyReal; 488 else if(theConfig->getSeparationEnergyTy 356 else if(theConfig->getSeparationEnergyType()==RealForLightSeparationEnergy) 489 getSeparationEnergy = getSeparationEne 357 getSeparationEnergy = getSeparationEnergyRealForLight; 490 else { 358 else { 491 INCL_FATAL("Unrecognized separation-en 359 INCL_FATAL("Unrecognized separation-energy type in ParticleTable initialization: " << theConfig->getSeparationEnergyType() << '\n'); 492 return; 360 return; 493 } 361 } 494 362 495 // Initialise the Fermi-momentum functio 363 // Initialise the Fermi-momentum function 496 if(!theConfig || theConfig->getFermiMome 364 if(!theConfig || theConfig->getFermiMomentumType()==ConstantFermiMomentum) { 497 getFermiMomentum = ParticleTable::getF 365 getFermiMomentum = ParticleTable::getFermiMomentumConstant; 498 if(theConfig) { 366 if(theConfig) { 499 const G4double aFermiMomentum = theC 367 const G4double aFermiMomentum = theConfig->getFermiMomentum(); 500 if(aFermiMomentum>0.) 368 if(aFermiMomentum>0.) 501 constantFermiMomentum = aFermiMome 369 constantFermiMomentum = aFermiMomentum; 502 else 370 else 503 constantFermiMomentum = PhysicalCo 371 constantFermiMomentum = PhysicalConstants::Pf; 504 } else { 372 } else { 505 constantFermiMomentum = PhysicalCons 373 constantFermiMomentum = PhysicalConstants::Pf; 506 } 374 } 507 } else if(theConfig->getFermiMomentumTyp 375 } else if(theConfig->getFermiMomentumType()==ConstantLightFermiMomentum) 508 getFermiMomentum = ParticleTable::getF 376 getFermiMomentum = ParticleTable::getFermiMomentumConstantLight; 509 else if(theConfig->getFermiMomentumType( 377 else if(theConfig->getFermiMomentumType()==MassDependentFermiMomentum) 510 getFermiMomentum = ParticleTable::getF 378 getFermiMomentum = ParticleTable::getFermiMomentumMassDependent; 511 else { 379 else { 512 INCL_FATAL("Unrecognized Fermi-momentu 380 INCL_FATAL("Unrecognized Fermi-momentum type in ParticleTable initialization: " << theConfig->getFermiMomentumType() << '\n'); 513 return; 381 return; 514 } 382 } 515 383 516 // Initialise the r-p correlation coeffi 384 // Initialise the r-p correlation coefficients 517 std::fill(rpCorrelationCoefficient, rpCo 385 std::fill(rpCorrelationCoefficient, rpCorrelationCoefficient + UnknownParticle, 1.); 518 if(theConfig) { 386 if(theConfig) { 519 rpCorrelationCoefficient[Proton] = the 387 rpCorrelationCoefficient[Proton] = theConfig->getRPCorrelationCoefficient(Proton); 520 rpCorrelationCoefficient[Neutron] = th 388 rpCorrelationCoefficient[Neutron] = theConfig->getRPCorrelationCoefficient(Neutron); 521 } 389 } 522 390 523 // Initialise the neutron-skin parameter 391 // Initialise the neutron-skin parameters 524 if(theConfig) { 392 if(theConfig) { 525 neutronSkin = theConfig->getNeutronSki 393 neutronSkin = theConfig->getNeutronSkin(); 526 neutronHalo = theConfig->getNeutronHal 394 neutronHalo = theConfig->getNeutronHalo(); 527 } 395 } 528 396 529 } 397 } 530 398 531 G4int getIsospin(const ParticleType t) { 399 G4int getIsospin(const ParticleType t) { 532 // Actually this is the 3rd component of 400 // Actually this is the 3rd component of isospin (I_z) multiplied by 2! 533 if(t == Proton) { 401 if(t == Proton) { 534 return 1; 402 return 1; 535 } else if(t == Neutron) { 403 } else if(t == Neutron) { 536 return -1; 404 return -1; 537 } else if(t == PiPlus) { 405 } else if(t == PiPlus) { 538 return 2; 406 return 2; 539 } else if(t == PiMinus) { 407 } else if(t == PiMinus) { 540 return -2; 408 return -2; 541 } else if(t == PiZero) { 409 } else if(t == PiZero) { 542 return 0; 410 return 0; 543 } else if(t == DeltaPlusPlus) { 411 } else if(t == DeltaPlusPlus) { 544 return 3; 412 return 3; 545 } else if(t == DeltaPlus) { 413 } else if(t == DeltaPlus) { 546 return 1; 414 return 1; 547 } else if(t == DeltaZero) { 415 } else if(t == DeltaZero) { 548 return -1; 416 return -1; 549 } else if(t == DeltaMinus) { 417 } else if(t == DeltaMinus) { 550 return -3; << 418 return -3; 551 } else if(t == Lambda) { << 552 return 0; << 553 } else if(t == SigmaPlus) { << 554 return 2; << 555 } else if(t == SigmaZero) { << 556 return 0; << 557 } else if(t == SigmaMinus) { << 558 return -2; << 559 } else if(t == KPlus) { << 560 return 1; << 561 } else if(t == KZero) { << 562 return -1; << 563 } else if(t == KZeroBar) { << 564 return 1; << 565 } else if(t == KShort) { << 566 return 0; << 567 } else if(t == KLong) { << 568 return 0; << 569 } else if(t == KMinus) { << 570 return -1; << 571 } else if(t == Eta) { 419 } else if(t == Eta) { 572 return 0; << 420 return 0; 573 } else if(t == Omega) { 421 } else if(t == Omega) { 574 return 0; << 422 return 0; 575 } else if(t == EtaPrime) { 423 } else if(t == EtaPrime) { 576 return 0; << 424 return 0; 577 } else if(t == Photon) { 425 } else if(t == Photon) { 578 return 0; << 426 return 0; 579 } else if(t == antiProton) { << 580 return -1; << 581 } else if(t == XiMinus) { << 582 return -1; << 583 } else if(t == XiZero) { << 584 return 1; << 585 } else if(t == antiNeutron) { << 586 return 1; << 587 } else if(t == antiLambda) { << 588 return 0; << 589 } else if(t == antiSigmaPlus) { << 590 return -2; << 591 } else if(t == antiSigmaZero) { << 592 return 0; << 593 } else if(t == antiSigmaMinus) { << 594 return 2; << 595 } else if(t == antiXiMinus) { << 596 return 1; << 597 } else if(t == antiXiZero) { << 598 return -1; << 599 } 427 } >> 428 600 INCL_ERROR("Requested isospin of an unkn 429 INCL_ERROR("Requested isospin of an unknown particle!"); 601 return -10; // Unknown 430 return -10; // Unknown 602 } 431 } 603 432 604 std::string getShortName(const ParticleSpe << 433 std::string getShortName(const ParticleSpecies &s) { 605 if(sp.theType==Composite && sp.theS == 0 << 434 if(s.theType==Composite) 606 return getShortName(sp.theA,sp.theZ); << 435 return getShortName(s.theA,s.theZ); 607 else if(sp.theType==Composite) << 608 return getName(sp.theA,sp.theZ,sp.theS << 609 else 436 else 610 return getShortName(sp.theType); << 437 return getShortName(s.theType); 611 } 438 } 612 << 439 613 std::string getName(const ParticleSpecies << 440 std::string getName(const ParticleSpecies &s) { 614 if(sp.theType==Composite && sp.theS == 0 << 441 if(s.theType==Composite) 615 return getName(sp.theA,sp.theZ); << 442 return getName(s.theA,s.theZ); 616 else if(sp.theType==Composite) << 617 return getName(sp.theA,sp.theZ,sp.theS << 618 else 443 else 619 return getName(sp.theType); << 444 return getName(s.theType); 620 } 445 } 621 446 622 std::string getName(const G4int A, const G 447 std::string getName(const G4int A, const G4int Z) { 623 std::stringstream stream; 448 std::stringstream stream; 624 stream << getElementName(Z) << "-" << A; 449 stream << getElementName(Z) << "-" << A; 625 return stream.str(); 450 return stream.str(); 626 } 451 } 627 452 628 std::string getName(const G4int A, const G << 629 std::stringstream stream; << 630 if(S >= 0) // S < 0 for hypernuclei << 631 return getName(A, Z); << 632 else if(S == -1) << 633 stream << getElementName(Z) << "-" << << 634 else << 635 stream << getElementName(Z) << "-" << << 636 return stream.str(); << 637 } << 638 << 639 std::string getShortName(const G4int A, co 453 std::string getShortName(const G4int A, const G4int Z) { 640 std::stringstream stream; 454 std::stringstream stream; 641 stream << getElementName(Z); 455 stream << getElementName(Z); 642 if(A>0) 456 if(A>0) 643 stream << A; 457 stream << A; 644 return stream.str(); 458 return stream.str(); 645 } 459 } 646 460 647 std::string getName(const ParticleType p) 461 std::string getName(const ParticleType p) { 648 if(p == G4INCL::Proton) { 462 if(p == G4INCL::Proton) { 649 return std::string("proton"); 463 return std::string("proton"); 650 } else if(p == G4INCL::Neutron) { 464 } else if(p == G4INCL::Neutron) { 651 return std::string("neutron"); 465 return std::string("neutron"); 652 } else if(p == G4INCL::DeltaPlusPlus) { 466 } else if(p == G4INCL::DeltaPlusPlus) { 653 return std::string("delta++"); 467 return std::string("delta++"); 654 } else if(p == G4INCL::DeltaPlus) { 468 } else if(p == G4INCL::DeltaPlus) { 655 return std::string("delta+"); 469 return std::string("delta+"); 656 } else if(p == G4INCL::DeltaZero) { 470 } else if(p == G4INCL::DeltaZero) { 657 return std::string("delta0"); 471 return std::string("delta0"); 658 } else if(p == G4INCL::DeltaMinus) { 472 } else if(p == G4INCL::DeltaMinus) { 659 return std::string("delta-"); 473 return std::string("delta-"); 660 } else if(p == G4INCL::PiPlus) { 474 } else if(p == G4INCL::PiPlus) { 661 return std::string("pi+"); 475 return std::string("pi+"); 662 } else if(p == G4INCL::PiZero) { 476 } else if(p == G4INCL::PiZero) { 663 return std::string("pi0"); 477 return std::string("pi0"); 664 } else if(p == G4INCL::PiMinus) { 478 } else if(p == G4INCL::PiMinus) { 665 return std::string("pi-"); 479 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) { 480 } else if(p == G4INCL::Composite) { 707 return std::string("composite"); 481 return std::string("composite"); 708 } else if(p == G4INCL::Eta) { 482 } else if(p == G4INCL::Eta) { 709 return std::string("eta"); 483 return std::string("eta"); 710 } else if(p == G4INCL::Omega) { 484 } else if(p == G4INCL::Omega) { 711 return std::string("omega"); 485 return std::string("omega"); 712 } else if(p == G4INCL::EtaPrime) { 486 } else if(p == G4INCL::EtaPrime) { 713 return std::string("etaprime"); 487 return std::string("etaprime"); 714 } else if(p == G4INCL::Photon) { 488 } else if(p == G4INCL::Photon) { 715 return std::string("photon"); << 489 return std::string("photon"); 716 } 490 } 717 return std::string("unknown"); 491 return std::string("unknown"); 718 } 492 } 719 493 720 std::string getShortName(const ParticleTyp 494 std::string getShortName(const ParticleType p) { 721 if(p == G4INCL::Proton) { 495 if(p == G4INCL::Proton) { 722 return std::string("p"); 496 return std::string("p"); 723 } else if(p == G4INCL::Neutron) { 497 } else if(p == G4INCL::Neutron) { 724 return std::string("n"); 498 return std::string("n"); 725 } else if(p == G4INCL::DeltaPlusPlus) { 499 } else if(p == G4INCL::DeltaPlusPlus) { 726 return std::string("d++"); 500 return std::string("d++"); 727 } else if(p == G4INCL::DeltaPlus) { 501 } else if(p == G4INCL::DeltaPlus) { 728 return std::string("d+"); 502 return std::string("d+"); 729 } else if(p == G4INCL::DeltaZero) { 503 } else if(p == G4INCL::DeltaZero) { 730 return std::string("d0"); 504 return std::string("d0"); 731 } else if(p == G4INCL::DeltaMinus) { 505 } else if(p == G4INCL::DeltaMinus) { 732 return std::string("d-"); 506 return std::string("d-"); 733 } else if(p == G4INCL::PiPlus) { 507 } else if(p == G4INCL::PiPlus) { 734 return std::string("pi+"); 508 return std::string("pi+"); 735 } else if(p == G4INCL::PiZero) { 509 } else if(p == G4INCL::PiZero) { 736 return std::string("pi0"); 510 return std::string("pi0"); 737 } else if(p == G4INCL::PiMinus) { 511 } else if(p == G4INCL::PiMinus) { 738 return std::string("pi-"); 512 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) { 513 } else if(p == G4INCL::Composite) { 780 return std::string("comp"); 514 return std::string("comp"); 781 } else if(p == G4INCL::Eta) { 515 } else if(p == G4INCL::Eta) { 782 return std::string("eta"); 516 return std::string("eta"); 783 } else if(p == G4INCL::Omega) { 517 } else if(p == G4INCL::Omega) { 784 return std::string("omega"); 518 return std::string("omega"); 785 } else if(p == G4INCL::EtaPrime) { 519 } else if(p == G4INCL::EtaPrime) { 786 return std::string("etap"); 520 return std::string("etap"); 787 } else if(p == G4INCL::Photon) { 521 } else if(p == G4INCL::Photon) { 788 return std::string("photon"); 522 return std::string("photon"); 789 } 523 } 790 return std::string("unknown"); 524 return std::string("unknown"); 791 } 525 } 792 526 793 G4double getINCLMass(const ParticleType pt 527 G4double getINCLMass(const ParticleType pt) { 794 if(pt == Proton) { 528 if(pt == Proton) { 795 return protonMass; 529 return protonMass; 796 } else if(pt == Neutron) { 530 } else if(pt == Neutron) { 797 return neutronMass; 531 return neutronMass; 798 } else if(pt == PiPlus) { 532 } else if(pt == PiPlus) { 799 return piPlusMass; 533 return piPlusMass; 800 } else if(pt == PiMinus) { 534 } else if(pt == PiMinus) { 801 return piMinusMass; 535 return piMinusMass; 802 } else if(pt == PiZero) { 536 } else if(pt == PiZero) { 803 return piZeroMass; 537 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) { 538 } else if(pt == Eta) { 845 return etaMass; 539 return etaMass; 846 } else if(pt == Omega) { 540 } else if(pt == Omega) { 847 return omegaMass; 541 return omegaMass; 848 } else if(pt == EtaPrime) { 542 } else if(pt == EtaPrime) { 849 return etaPrimeMass; 543 return etaPrimeMass; 850 } else if(pt == Photon) { 544 } else if(pt == Photon) { 851 return photonMass; << 545 return photonMass; 852 } else { 546 } else { 853 INCL_ERROR("getMass : Unknown particle 547 INCL_ERROR("getMass : Unknown particle type." << '\n'); 854 return 0.0; 548 return 0.0; 855 } 549 } 856 } 550 } 857 << 551 858 G4double getRealMass(const ParticleType t) 552 G4double getRealMass(const ParticleType t) { 859 switch(t) { 553 switch(t) { 860 case Proton: 554 case Proton: 861 return theRealProtonMass; 555 return theRealProtonMass; 862 break; 556 break; 863 case Neutron: 557 case Neutron: 864 return theRealNeutronMass; 558 return theRealNeutronMass; 865 break; 559 break; 866 case PiPlus: 560 case PiPlus: 867 case PiMinus: 561 case PiMinus: 868 return theRealChargedPiMass; 562 return theRealChargedPiMass; 869 break; 563 break; 870 case PiZero: 564 case PiZero: 871 return theRealPiZeroMass; 565 return theRealPiZeroMass; 872 break; 566 break; 873 case Eta: 567 case Eta: 874 return theRealEtaMass; 568 return theRealEtaMass; 875 break; 569 break; 876 case Omega: 570 case Omega: 877 return theRealOmegaMass; 571 return theRealOmegaMass; 878 break; 572 break; 879 case EtaPrime: 573 case EtaPrime: 880 return theRealEtaPrimeMass; 574 return theRealEtaPrimeMass; 881 break; 575 break; 882 case Photon: 576 case Photon: 883 return theRealPhotonMass; 577 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; 578 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: 579 default: 938 INCL_ERROR("Particle::getRealMass : 580 INCL_ERROR("Particle::getRealMass : Unknown particle type." << '\n'); 939 return 0.0; 581 return 0.0; 940 break; 582 break; 941 } 583 } 942 } 584 } 943 << 585 944 G4double getRealMass(const G4int A, const << 586 G4double getRealMass(const G4int A, const G4int Z) { 945 // assert(A>=0); 587 // assert(A>=0); 946 // For nuclei with Z<0 or Z>A, assume th 588 // For nuclei with Z<0 or Z>A, assume that the exotic charge state is due to pions 947 if(Z<0 && S<0) << 948 return (A+S)*theRealNeutronMass - S*La << 949 else if(Z>A && S<0) << 950 return (A+S)*theRealProtonMass - S*Lam << 951 if(Z<0) 589 if(Z<0) 952 return (A)*theRealNeutronMass - Z*getR << 590 return A*neutronMass - Z*getRealMass(PiMinus); 953 else if(Z>A) 591 else if(Z>A) 954 return (A)*theRealProtonMass + (A-Z)*g << 592 return A*protonMass + (A-Z)*getRealMass(PiPlus); 955 else if(Z==0 && S==0) << 593 else if(Z==0) 956 return A*theRealNeutronMass; << 594 return A*getRealMass(Neutron); 957 else if(A==Z) 595 else if(A==Z) 958 return A*theRealProtonMass; << 596 return A*getRealMass(Proton); 959 else if(Z==0 && S<0) << 960 return (A+S)*theRealNeutronMass-S*Lamb << 961 else if(A>1) { 597 else if(A>1) { 962 #ifndef INCLXX_IN_GEANT4_MODE 598 #ifndef INCLXX_IN_GEANT4_MODE 963 return ::G4INCL::NuclearMassTable::get << 599 return ::G4INCL::NuclearMassTable::getMass(A,Z); 964 #else 600 #else 965 if(S<0) return theG4IonTable->GetNucle << 601 return theG4IonTable->GetNucleusMass(Z,A) / MeV; 966 else return theG4IonTable->GetNucle << 967 #endif 602 #endif 968 } else 603 } else 969 return 0.; 604 return 0.; 970 } 605 } 971 606 972 G4double getINCLMass(const G4int A, const << 607 G4double getINCLMass(const G4int A, const G4int Z) { 973 // assert(A>=0); 608 // assert(A>=0); 974 // For nuclei with Z<0 or Z>A, assume th 609 // 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 << 610 if(Z<0) 976 if(Z<0 && S<0) << 611 return A*neutronMass - Z*getINCLMass(PiMinus); 977 return (A+S)*neutronMass - S*LambdaMas << 978 else if(Z>A && S<0) << 979 return (A+S)*protonMass - S*LambdaMass << 980 else if(Z<0) << 981 return (A)*neutronMass - Z*getINCLMass << 982 else if(Z>A) 612 else if(Z>A) 983 return (A)*protonMass + (A-Z)*getINCLM << 613 return A*protonMass + (A-Z)*getINCLMass(PiPlus); 984 else if(A>1 && S<0) << 985 return Z*(protonMass - protonSeparatio << 986 else if(A>1) 614 else if(A>1) 987 return Z*(protonMass - protonSeparatio 615 return Z*(protonMass - protonSeparationEnergy) + (A-Z)*(neutronMass - neutronSeparationEnergy); 988 else if(A==1 && Z==0 && S==0) << 616 else if(A==1 && Z==0) 989 return getINCLMass(Neutron); 617 return getINCLMass(Neutron); 990 else if(A==1 && Z==1 && S==0) << 618 else if(A==1 && Z==1) 991 return getINCLMass(Proton); 619 return getINCLMass(Proton); 992 else if(A==1 && Z==0 && S==-1) << 993 return getINCLMass(Lambda); << 994 else 620 else 995 return 0.; 621 return 0.; 996 } 622 } 997 623 998 G4double getTableQValue(const G4int A1, co << 624 G4double getTableQValue(const G4int A1, const G4int Z1, const G4int A2, const G4int Z2) { 999 return getTableMass(A1,Z1,S1) + getTable << 625 return getTableMass(A1,Z1) + getTableMass(A2,Z2) - getTableMass(A1+A2,Z1+Z2); 1000 } 626 } 1001 627 1002 G4double getTableQValue(const G4int A1, c << 628 G4double getTableQValue(const G4int A1, const G4int Z1, const G4int A2, const G4int Z2, const G4int A3, const G4int Z3) { 1003 return getTableMass(A1,Z1,S1) + getTabl << 629 return getTableMass(A1,Z1) + getTableMass(A2,Z2) - getTableMass(A3,Z3) - getTableMass(A1+A2-A3,Z1+Z2-Z3); 1004 } 630 } 1005 631 1006 G4double getTableSpeciesMass(const Partic 632 G4double getTableSpeciesMass(const ParticleSpecies &p) { 1007 if(p.theType == Composite) 633 if(p.theType == Composite) 1008 return (*getTableMass)(p.theA, p.theZ << 634 return (*getTableMass)(p.theA, p.theZ); 1009 else 635 else 1010 return (*getTableParticleMass)(p.theT 636 return (*getTableParticleMass)(p.theType); 1011 } 637 } 1012 638 1013 G4int getMassNumber(const ParticleType t) << 639 G4int getMassNumber(const ParticleType t) { 1014 << 1015 switch(t) { 640 switch(t) { 1016 case Proton: 641 case Proton: 1017 case Neutron: 642 case Neutron: 1018 case DeltaPlusPlus: 643 case DeltaPlusPlus: 1019 case DeltaPlus: 644 case DeltaPlus: 1020 case DeltaZero: 645 case DeltaZero: 1021 case DeltaMinus: 646 case DeltaMinus: 1022 case SigmaPlus: << 1023 case SigmaZero: << 1024 case SigmaMinus: << 1025 case Lambda: << 1026 case XiZero: << 1027 case XiMinus: << 1028 return 1; 647 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; 648 break; 1040 case PiPlus: 649 case PiPlus: 1041 case PiMinus: 650 case PiMinus: 1042 case PiZero: 651 case PiZero: 1043 case KPlus: << 1044 case KZero: << 1045 case KZeroBar: << 1046 case KShort: << 1047 case KLong: << 1048 case KMinus: << 1049 case Eta: 652 case Eta: 1050 case Omega: 653 case Omega: 1051 case EtaPrime: 654 case EtaPrime: 1052 case Photon: 655 case Photon: 1053 return 0; 656 return 0; 1054 break; 657 break; 1055 default: 658 default: 1056 return 0; 659 return 0; 1057 break; 660 break; 1058 } 661 } 1059 } 662 } 1060 663 1061 G4int getChargeNumber(const ParticleType 664 G4int getChargeNumber(const ParticleType t) { 1062 switch(t) { 665 switch(t) { 1063 case DeltaPlusPlus: 666 case DeltaPlusPlus: 1064 return 2; 667 return 2; 1065 break; 668 break; 1066 case Proton: 669 case Proton: 1067 case DeltaPlus: 670 case DeltaPlus: 1068 case PiPlus: 671 case PiPlus: 1069 case SigmaPlus: << 1070 case KPlus: << 1071 case antiSigmaMinus: << 1072 case antiXiMinus: << 1073 return 1; 672 return 1; 1074 break; << 673 break; 1075 case Neutron: 674 case Neutron: 1076 case DeltaZero: 675 case DeltaZero: 1077 case PiZero: 676 case PiZero: 1078 case SigmaZero: << 1079 case Lambda: << 1080 case KZero: << 1081 case KZeroBar: << 1082 case KShort: << 1083 case KLong: << 1084 case Eta: 677 case Eta: 1085 case Omega: 678 case Omega: 1086 case EtaPrime: 679 case EtaPrime: 1087 case Photon: 680 case Photon: 1088 case XiZero: << 1089 case antiNeutron: << 1090 case antiLambda: << 1091 case antiSigmaZero: << 1092 case antiXiZero: << 1093 return 0; 681 return 0; 1094 break; 682 break; 1095 case DeltaMinus: 683 case DeltaMinus: 1096 case PiMinus: 684 case PiMinus: 1097 case SigmaMinus: << 1098 case KMinus: << 1099 case antiProton: << 1100 case XiMinus: << 1101 case antiSigmaPlus: << 1102 return -1; 685 return -1; 1103 break; << 1104 default: << 1105 return 0; << 1106 break; 686 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: 687 default: 1160 return 0; 688 return 0; 1161 break; 689 break; 1162 } 690 } 1163 } 691 } 1164 692 1165 G4double getNuclearRadius(const ParticleT 693 G4double getNuclearRadius(const ParticleType t, const G4int A, const G4int Z) { 1166 // assert(A>=0); 694 // assert(A>=0); 1167 if(A > 19 || (A < 6 && A >= 2)) { << 695 if(A >= 19 || (A < 6 && A >= 2)) { 1168 // For large (Woods-Saxon or Modified 696 // For large (Woods-Saxon or Modified Harmonic Oscillator) or small 1169 // (Gaussian) nuclei, the radius para 697 // (Gaussian) nuclei, the radius parameter is just the nuclear radius 1170 return getRadiusParameter(t,A,Z); 698 return getRadiusParameter(t,A,Z); 1171 } else if(A < clusterTableASize && Z>=0 699 } else if(A < clusterTableASize && Z>=0 && Z < clusterTableZSize && A >= 6) { 1172 const G4double thisRMS = positionRMS[ 700 const G4double thisRMS = positionRMS[Z][A]; 1173 if(thisRMS>0.0) 701 if(thisRMS>0.0) 1174 return thisRMS; 702 return thisRMS; 1175 else { 703 else { 1176 INCL_DEBUG("getNuclearRadius: Radiu 704 INCL_DEBUG("getNuclearRadius: Radius for nucleus A = " << A << " Z = " << Z << " is not available" << '\n' 1177 << "returning radius for 705 << "returning radius for C12"); 1178 return positionRMS[6][12]; 706 return positionRMS[6][12]; 1179 } 707 } 1180 } else if(A <= 19) { << 708 } else if(A < 19) { 1181 const G4double theRadiusParameter = g 709 const G4double theRadiusParameter = getRadiusParameter(t, A, Z); 1182 const G4double theDiffusenessParamete 710 const G4double theDiffusenessParameter = getSurfaceDiffuseness(t, A, Z); 1183 // The formula yields the nuclear RMS 711 // The formula yields the nuclear RMS radius based on the parameters of 1184 // the nuclear-density function 712 // the nuclear-density function 1185 return 1.225*theDiffusenessParameter* 713 return 1.225*theDiffusenessParameter* 1186 std::sqrt((2.+5.*theRadiusParameter 714 std::sqrt((2.+5.*theRadiusParameter)/(2.+3.*theRadiusParameter)); 1187 } else { 715 } else { 1188 INCL_ERROR("getNuclearRadius: No radi 716 INCL_ERROR("getNuclearRadius: No radius for nucleus A = " << A << " Z = " << Z << '\n'); 1189 return 0.0; 717 return 0.0; 1190 } 718 } 1191 } 719 } 1192 720 1193 G4double getLargestNuclearRadius(const G4 721 G4double getLargestNuclearRadius(const G4int A, const G4int Z) { 1194 return Math::max(getNuclearRadius(Proto 722 return Math::max(getNuclearRadius(Proton, A, Z), getNuclearRadius(Neutron, A, Z)); 1195 } 723 } 1196 724 1197 G4double getRadiusParameter(const Particl 725 G4double getRadiusParameter(const ParticleType t, const G4int A, const G4int Z) { 1198 // assert(A>0); 726 // assert(A>0); 1199 if(A > 19) { << 727 if(A >= 28) { 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 728 // phenomenological radius fit 1206 G4double r0 = (2.745e-4 * A + 1.063) 729 G4double r0 = (2.745e-4 * A + 1.063) * std::pow(A, 1.0/3.0); 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) 730 if(t==Neutron) 1214 r0 += neutronSkin; 731 r0 += neutronSkin; 1215 return r0; 732 return r0; 1216 } else if(A < 6 && A >= 2) { 733 } else if(A < 6 && A >= 2) { 1217 if(Z<clusterTableZSize && Z>=0) { 734 if(Z<clusterTableZSize && Z>=0) { 1218 const G4double thisRMS = positionRM 735 const G4double thisRMS = positionRMS[Z][A]; 1219 if(thisRMS>0.0) 736 if(thisRMS>0.0) 1220 return thisRMS; 737 return thisRMS; 1221 else { 738 else { 1222 INCL_DEBUG("getRadiusParameter: R 739 INCL_DEBUG("getRadiusParameter: Radius for nucleus A = " << A << " Z = " << Z << " is not available" << '\n' 1223 << "returning radius f 740 << "returning radius for C12"); 1224 return positionRMS[6][12]; 741 return positionRMS[6][12]; 1225 } 742 } 1226 } else { 743 } else { 1227 INCL_DEBUG("getRadiusParameter: Rad 744 INCL_DEBUG("getRadiusParameter: Radius for nucleus A = " << A << " Z = " << Z << " is not available" << '\n' 1228 << "returning radius for 745 << "returning radius for C12"); 1229 return positionRMS[6][12]; 746 return positionRMS[6][12]; 1230 } 747 } 1231 } else if(A <= 19 && A >= 6) { << 748 } else if(A < 28 && 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]; 749 return mediumRadius[A-1]; 1242 // return 1.581*mediumDiffusenes 750 // return 1.581*mediumDiffuseness[A-1]*(2.+5.*mediumRadius[A-1])/(2.+3.*mediumRadius[A-1]); 1243 } else { 751 } else { 1244 INCL_ERROR("getRadiusParameter: No ra 752 INCL_ERROR("getRadiusParameter: No radius for nucleus A = " << A << " Z = " << Z << '\n'); 1245 return 0.0; 753 return 0.0; 1246 } 754 } 1247 } 755 } 1248 756 1249 G4double getMaximumNuclearRadius(const Pa 757 G4double getMaximumNuclearRadius(const ParticleType t, const G4int A, const G4int Z) { 1250 const G4double XFOISA = 8.0; 758 const G4double XFOISA = 8.0; 1251 if(A > 19) { << 759 if(A >= 19) { 1252 return getNuclearRadius(t,A,Z) + XFOI 760 return getNuclearRadius(t,A,Z) + XFOISA * getSurfaceDiffuseness(t,A,Z); 1253 } else if(A <= 19 && A >= 6) { << 761 } else if(A < 19 && A >= 6) { 1254 return 5.5 + 0.3 * (G4double(A) - 6.0 762 return 5.5 + 0.3 * (G4double(A) - 6.0)/12.0; 1255 } else if(A >= 2) { 763 } else if(A >= 2) { 1256 return getNuclearRadius(t, A, Z) + 4. 764 return getNuclearRadius(t, A, Z) + 4.5; 1257 } else { 765 } else { 1258 INCL_ERROR("getMaximumNuclearRadius : 766 INCL_ERROR("getMaximumNuclearRadius : No maximum radius for nucleus A = " << A << " Z = " << Z << '\n'); 1259 return 0.0; 767 return 0.0; 1260 } 768 } 1261 } 769 } 1262 770 1263 G4double getSurfaceDiffuseness(const Part 771 G4double getSurfaceDiffuseness(const ParticleType t, const G4int A, const G4int Z) { 1264 if(A > 19) { << 772 if(A >= 28) { 1265 // phenomenological fit << 1266 G4double a = 1.63e-4 * A + 0.510; 773 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) 774 if(t==Neutron) 1279 a += neutronHalo; 775 a += neutronHalo; 1280 return a; 776 return a; 1281 } else if(A <= 19 && A >= 6) { << 777 } else if(A < 28 && A >= 19) { 1282 // HFB calculations << 778 return mediumDiffuseness[A-1]; 1283 if(getRPCorrelationCoefficient(t)<1.) << 779 } else if(A < 19 && A >= 6) { 1284 G4double ahfb = HFB::getRadiusParam << 1285 if(ahfb>0.)return ahfb; << 1286 } << 1287 return mediumDiffuseness[A-1]; 780 return mediumDiffuseness[A-1]; 1288 } else if(A < 6 && A >= 2) { 781 } else if(A < 6 && A >= 2) { 1289 INCL_ERROR("getSurfaceDiffuseness: wa 782 INCL_ERROR("getSurfaceDiffuseness: was called for A = " << A << " Z = " << Z << '\n'); 1290 return 0.0; 783 return 0.0; 1291 } else { 784 } else { 1292 INCL_ERROR("getSurfaceDiffuseness: No 785 INCL_ERROR("getSurfaceDiffuseness: No diffuseness for nucleus A = " << A << " Z = " << Z << '\n'); 1293 return 0.0; 786 return 0.0; 1294 } 787 } 1295 } 788 } 1296 789 1297 G4double getMomentumRMS(const G4int A, co 790 G4double getMomentumRMS(const G4int A, const G4int Z) { 1298 // assert(Z>=0 && A>=0 && Z<=A); 791 // assert(Z>=0 && A>=0 && Z<=A); 1299 return getFermiMomentum(A,Z) * Math::sq 792 return getFermiMomentum(A,Z) * Math::sqrtThreeFifths; 1300 } 793 } 1301 794 1302 G4double getSeparationEnergyINCL(const Pa 795 G4double getSeparationEnergyINCL(const ParticleType t, const G4int /*A*/, const G4int /*Z*/) { 1303 if(t==Proton) 796 if(t==Proton) 1304 return theINCLProtonSeparationEnergy; 797 return theINCLProtonSeparationEnergy; 1305 else if(t==Neutron) 798 else if(t==Neutron) 1306 return theINCLNeutronSeparationEnergy 799 return theINCLNeutronSeparationEnergy; 1307 else if(t==Lambda) << 1308 return theINCLLambdaSeparationEnergy; << 1309 else if(t==antiProton) << 1310 return theINCLantiProtonSeparationEne << 1311 else { 800 else { 1312 INCL_ERROR("ParticleTable::getSeparat 801 INCL_ERROR("ParticleTable::getSeparationEnergyINCL : Unknown particle type." << '\n'); 1313 return 0.0; 802 return 0.0; 1314 } 803 } 1315 } 804 } 1316 805 1317 G4double getSeparationEnergyReal(const Pa 806 G4double getSeparationEnergyReal(const ParticleType t, const G4int A, const G4int Z) { 1318 // Real separation energies for all nuc 807 // Real separation energies for all nuclei 1319 if(t==Proton) 808 if(t==Proton) 1320 return (*getTableParticleMass)(Proton << 809 return (*getTableParticleMass)(Proton) + (*getTableMass)(A-1,Z-1) - (*getTableMass)(A,Z); 1321 else if(t==Neutron) 810 else if(t==Neutron) 1322 return (*getTableParticleMass)(Neutro << 811 return (*getTableParticleMass)(Neutron) + (*getTableMass)(A-1,Z) - (*getTableMass)(A,Z); 1323 else if(t==Lambda) << 1324 return (*getTableParticleMass)(Lambda << 1325 else { 812 else { 1326 INCL_ERROR("ParticleTable::getSeparat 813 INCL_ERROR("ParticleTable::getSeparationEnergyReal : Unknown particle type." << '\n'); 1327 return 0.0; 814 return 0.0; 1328 } 815 } 1329 } 816 } 1330 817 1331 G4double getSeparationEnergyRealForLight( 818 G4double getSeparationEnergyRealForLight(const ParticleType t, const G4int A, const G4int Z) { 1332 // Real separation energies for light n 819 // Real separation energies for light nuclei, fixed values for heavy nuclei 1333 if(Z<clusterTableZSize && A<clusterTabl 820 if(Z<clusterTableZSize && A<clusterTableASize) 1334 return getSeparationEnergyReal(t, A, 821 return getSeparationEnergyReal(t, A, Z); 1335 else 822 else 1336 return getSeparationEnergyINCL(t, A, 823 return getSeparationEnergyINCL(t, A, Z); 1337 } 824 } 1338 825 1339 G4double getProtonSeparationEnergy() { re 826 G4double getProtonSeparationEnergy() { return protonSeparationEnergy; } 1340 827 1341 G4double getNeutronSeparationEnergy() { r 828 G4double getNeutronSeparationEnergy() { return neutronSeparationEnergy; } 1342 829 1343 G4double getLambdaSeparationEnergy() { re << 830 void setProtonSeparationEnergy(const G4double s) { protonSeparationEnergy = s; } 1344 << 1345 void setProtonSeparationEnergy(const G4do << 1346 831 1347 void setNeutronSeparationEnergy(const G4d << 832 void setNeutronSeparationEnergy(const G4double s) { neutronSeparationEnergy = s; } 1348 << 1349 void setLambdaSeparationEnergy(const G4do << 1350 833 1351 std::string getElementName(const G4int Z) 834 std::string getElementName(const G4int Z) { 1352 if(Z<1) { 835 if(Z<1) { 1353 INCL_WARN("getElementName called with 836 INCL_WARN("getElementName called with Z<1" << '\n'); 1354 return elementTable[0]; 837 return elementTable[0]; 1355 } else if(Z<elementTableSize) 838 } else if(Z<elementTableSize) 1356 return elementTable[Z]; 839 return elementTable[Z]; 1357 else 840 else 1358 return getIUPACElementName(Z); 841 return getIUPACElementName(Z); 1359 } 842 } 1360 843 1361 std::string getIUPACElementName(const G4i 844 std::string getIUPACElementName(const G4int Z) { 1362 std::stringstream elementStream; 845 std::stringstream elementStream; 1363 elementStream << Z; 846 elementStream << Z; 1364 std::string elementName = elementStream 847 std::string elementName = elementStream.str(); 1365 std::transform(elementName.begin(), ele 848 std::transform(elementName.begin(), elementName.end(), elementName.begin(), intToIUPAC); 1366 elementName[0] = (char)std::toupper(ele << 849 elementName[0] = std::toupper(elementName.at(0)); 1367 return elementName; 850 return elementName; 1368 } 851 } 1369 852 1370 G4int parseElement(std::string pS) { 853 G4int parseElement(std::string pS) { 1371 // Normalize the element name 854 // Normalize the element name 1372 std::transform(pS.begin(), pS.end(), pS 855 std::transform(pS.begin(), pS.end(), pS.begin(), ::tolower); 1373 pS[0] = (char)std::toupper(pS[0]); << 856 pS[0] = ::toupper(pS[0]); 1374 857 1375 const std::string *iter = std::find(ele 858 const std::string *iter = std::find(elementTable, elementTable+elementTableSize, pS); 1376 if(iter != elementTable+elementTableSiz 859 if(iter != elementTable+elementTableSize) 1377 return G4int(iter - elementTable); << 860 return iter - elementTable; 1378 else 861 else 1379 return ParticleTable::parseIUPACEleme 862 return ParticleTable::parseIUPACElement(pS); 1380 } 863 } 1381 864 1382 G4int parseIUPACElement(std::string const << 865 G4int parseIUPACElement(std::string const &s) { 1383 // Normalise to lower case 866 // Normalise to lower case 1384 std::string elementName(sel); << 867 std::string elementName(s); 1385 std::transform(elementName.begin(), ele 868 std::transform(elementName.begin(), elementName.end(), elementName.begin(), ::tolower); 1386 // Return 0 if the element name contain 869 // Return 0 if the element name contains anything but IUPAC digits 1387 if(elementName.find_first_not_of(elemen 870 if(elementName.find_first_not_of(elementIUPACDigits)!=std::string::npos) 1388 return 0; 871 return 0; 1389 std::transform(elementName.begin(), ele 872 std::transform(elementName.begin(), elementName.end(), elementName.begin(), iupacToInt); 1390 std::stringstream elementStream(element 873 std::stringstream elementStream(elementName); 1391 G4int Z; 874 G4int Z; 1392 elementStream >> Z; 875 elementStream >> Z; 1393 return Z; 876 return Z; 1394 } 877 } 1395 878 1396 IsotopicDistribution const &getNaturalIso 879 IsotopicDistribution const &getNaturalIsotopicDistribution(const G4int Z) { 1397 return getNaturalIsotopicDistributions( 880 return getNaturalIsotopicDistributions()->getIsotopicDistribution(Z); 1398 } 881 } 1399 882 1400 G4int drawRandomNaturalIsotope(const G4in 883 G4int drawRandomNaturalIsotope(const G4int Z) { 1401 return getNaturalIsotopicDistributions( 884 return getNaturalIsotopicDistributions()->drawRandomIsotope(Z); 1402 } 885 } 1403 886 1404 G4double getFermiMomentumConstant(const G 887 G4double getFermiMomentumConstant(const G4int /*A*/, const G4int /*Z*/) { 1405 return constantFermiMomentum; 888 return constantFermiMomentum; 1406 } 889 } 1407 890 1408 G4double getFermiMomentumConstantLight(co 891 G4double getFermiMomentumConstantLight(const G4int A, const G4int Z) { 1409 // assert(Z>0 && A>0 && Z<=A); 892 // assert(Z>0 && A>0 && Z<=A); 1410 if(Z<clusterTableZSize && A<clusterTabl 893 if(Z<clusterTableZSize && A<clusterTableASize) { 1411 const G4double rms = momentumRMS[Z][A 894 const G4double rms = momentumRMS[Z][A]; 1412 return ((rms>0.) ? rms : momentumRMS[ 895 return ((rms>0.) ? rms : momentumRMS[6][12]) * Math::sqrtFiveThirds; 1413 } else 896 } else 1414 return getFermiMomentumConstant(A,Z); 897 return getFermiMomentumConstant(A,Z); 1415 } 898 } 1416 899 1417 G4double getFermiMomentumMassDependent(co 900 G4double getFermiMomentumMassDependent(const G4int A, const G4int /*Z*/) { 1418 // assert(A>0); 901 // assert(A>0); 1419 static const G4double alphaParam = 259. 902 static const G4double alphaParam = 259.416; // MeV/c 1420 static const G4double betaParam = 152. 903 static const G4double betaParam = 152.824; // MeV/c 1421 static const G4double gammaParam = 9.51 904 static const G4double gammaParam = 9.5157E-2; 1422 return alphaParam - betaParam*std::exp( 905 return alphaParam - betaParam*std::exp(-gammaParam*((G4double)A)); 1423 } 906 } 1424 907 1425 G4double getRPCorrelationCoefficient(cons 908 G4double getRPCorrelationCoefficient(const ParticleType t) { 1426 // assert(t==Proton || t==Neutron || t==Lambd << 909 // assert(t==Proton || t==Neutron); 1427 return rpCorrelationCoefficient[t]; 910 return rpCorrelationCoefficient[t]; 1428 } 911 } 1429 912 1430 G4double getNeutronSkin() { return neutro 913 G4double getNeutronSkin() { return neutronSkin; } 1431 914 1432 G4double getNeutronHalo() { return neutro 915 G4double getNeutronHalo() { return neutronHalo; } 1433 916 1434 G4ThreadLocal G4double minDeltaMass = 0.; 917 G4ThreadLocal G4double minDeltaMass = 0.; 1435 G4ThreadLocal G4double minDeltaMass2 = 0. 918 G4ThreadLocal G4double minDeltaMass2 = 0.; 1436 G4ThreadLocal G4double minDeltaMassRndm = 919 G4ThreadLocal G4double minDeltaMassRndm = 0.; 1437 G4ThreadLocal NuclearMassFn getTableMass 920 G4ThreadLocal NuclearMassFn getTableMass = NULL; 1438 G4ThreadLocal ParticleMassFn getTablePart 921 G4ThreadLocal ParticleMassFn getTableParticleMass = NULL; 1439 G4ThreadLocal SeparationEnergyFn getSepar 922 G4ThreadLocal SeparationEnergyFn getSeparationEnergy = NULL; 1440 G4ThreadLocal FermiMomentumFn getFermiMom 923 G4ThreadLocal FermiMomentumFn getFermiMomentum = NULL; 1441 924 1442 ParticleType getPionType(const G4int isos 925 ParticleType getPionType(const G4int isosp) { 1443 // assert(isosp == -2 || isosp == 0 || isosp 926 // assert(isosp == -2 || isosp == 0 || isosp == 2); 1444 if (isosp == -2) { 927 if (isosp == -2) { 1445 return PiMinus; 928 return PiMinus; 1446 } 929 } 1447 else if (isosp == 0) { 930 else if (isosp == 0) { 1448 return PiZero; 931 return PiZero; 1449 } 932 } 1450 else { 933 else { 1451 return PiPlus; 934 return PiPlus; 1452 } 935 } 1453 } 936 } 1454 937 1455 ParticleType getNucleonType(const G4int i 938 ParticleType getNucleonType(const G4int isosp) { 1456 // assert(isosp == -1 || isosp == 1); 939 // assert(isosp == -1 || isosp == 1); 1457 if (isosp == -1) { 940 if (isosp == -1) { 1458 return Neutron; 941 return Neutron; 1459 } 942 } 1460 else { 943 else { 1461 return Proton; 944 return Proton; 1462 } 945 } 1463 } 946 } 1464 947 1465 ParticleType getDeltaType(const G4int iso 948 ParticleType getDeltaType(const G4int isosp) { 1466 // assert(isosp == -3 || isosp == -1 || isosp 949 // assert(isosp == -3 || isosp == -1 || isosp == 1 || isosp == 3); 1467 if (isosp == -3) { 950 if (isosp == -3) { 1468 return DeltaMinus; 951 return DeltaMinus; 1469 } 952 } 1470 else if (isosp == -1) { 953 else if (isosp == -1) { 1471 return DeltaZero; 954 return DeltaZero; 1472 } 955 } 1473 else if (isosp == 1) { 956 else if (isosp == 1) { 1474 return DeltaPlus; 957 return DeltaPlus; 1475 } 958 } 1476 else { 959 else { 1477 return DeltaPlusPlus; 960 return DeltaPlusPlus; 1478 } 961 } 1479 } 962 } 1480 963 1481 ParticleType getSigmaType(const G4int iso << 964 G4double getWidth(const ParticleType pt) { 1482 // assert(isosp == -2 || isosp == 0 || isosp << 965 // assert(pt == PiPlus || pt == PiMinus || pt == PiZero || pt == Eta || pt == Omega || pt == EtaPrime); 1483 if (isosp == -2) { << 966 if(pt == PiPlus) { 1484 return SigmaMinus; << 967 return piPlusWidth; 1485 } << 968 } else if(pt == PiMinus) { 1486 else if (isosp == 0) { << 969 return piMinusWidth; 1487 return SigmaZero; << 970 } else if(pt == PiZero) { 1488 } << 971 return piZeroWidth; 1489 else { << 972 } else if(pt == Eta) { 1490 return SigmaPlus; << 973 return etaWidth; 1491 } << 974 } else if(pt == Omega) { 1492 } << 975 return omegaWidth; 1493 << 976 } else if(pt == EtaPrime) { 1494 ParticleType getXiType(const G4int isosp) << 977 return etaPrimeWidth; 1495 // assert(isosp == -1 || isosp == 1); << 978 } else { 1496 if (isosp == -1) { << 979 INCL_ERROR("getWidth : Unknown particle type." << '\n'); 1497 return XiMinus; << 980 return 0.0; 1498 } << 981 } 1499 else { << 982 } 1500 return XiZero; << 983 1501 } << 984 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 << 1557 G4double getWidth(const ParticleType pt) << 1558 // assert(pt == PiPlus || pt == PiMinus || pt << 1559 if(pt == PiPlus) { << 1560 return piPlusWidth; << 1561 } else if(pt == PiMinus) { << 1562 return piMinusWidth; << 1563 } else if(pt == PiZero) { << 1564 return piZeroWidth; << 1565 } else if(pt == Eta) { << 1566 return etaWidth; << 1567 } else if(pt == Omega) { << 1568 return omegaWidth; << 1569 } else if(pt == EtaPrime) { << 1570 return etaPrimeWidth; << 1571 } else if(pt == SigmaPlus) { << 1572 return SigmaPlusWidth; << 1573 } else if(pt == SigmaZero) { << 1574 return SigmaZeroWidth; << 1575 } else if(pt == SigmaMinus) { << 1576 return SigmaMinusWidth; << 1577 } else if(pt == KPlus) { << 1578 return KPlusWidth; << 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 985 } // namespace ParticleTable 1610 } // namespace G4INCL 986 } // namespace G4INCL 1611 987 1612 988