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