| 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 // 26 //
27 // ------------------------------------------- 27 // -------------------------------------------------------------------
28 // 28 //
29 // GEANT4 Class 29 // GEANT4 Class
30 // 30 //
31 // File name: G4EmCorrections 31 // File name: G4EmCorrections
32 // 32 //
33 // Author: Vladimir Ivanchenko 33 // Author: Vladimir Ivanchenko
34 // 34 //
35 // Creation date: 13.01.2005 35 // Creation date: 13.01.2005
36 // 36 //
37 // Modifications: 37 // Modifications:
38 // 05.05.2005 V.Ivanchenko Fix misprint in Mot 38 // 05.05.2005 V.Ivanchenko Fix misprint in Mott term
39 // 26.11.2005 V.Ivanchenko Fix effective charg 39 // 26.11.2005 V.Ivanchenko Fix effective charge for heavy ions using original paper
40 // 28.04.2006 V.Ivanchenko General cleanup, ad 40 // 28.04.2006 V.Ivanchenko General cleanup, add finite size corrections
41 // 13.05.2006 V.Ivanchenko Add corrections for 41 // 13.05.2006 V.Ivanchenko Add corrections for ion stopping
42 // 08.05.2007 V.Ivanchenko Use G4IonTable for 42 // 08.05.2007 V.Ivanchenko Use G4IonTable for ion mass instead of NistTable to avoid
43 // division by zero 43 // division by zero
44 // 29.02.2008 V.Ivanchenko use expantions for 44 // 29.02.2008 V.Ivanchenko use expantions for log and power function
45 // 21.04.2008 Updated computations for ions (V 45 // 21.04.2008 Updated computations for ions (V.Ivanchenko)
46 // 20.05.2008 Removed Finite Size correction ( 46 // 20.05.2008 Removed Finite Size correction (V.Ivanchenko)
47 // 19.04.2012 Fix reproducibility problem (A.R 47 // 19.04.2012 Fix reproducibility problem (A.Ribon)
48 // 48 //
49 // 49 //
50 // Class Description: 50 // Class Description:
51 // 51 //
52 // This class provides calculation of EM corre 52 // This class provides calculation of EM corrections to ionisation
53 // 53 //
54 54
55 // ------------------------------------------- 55 // -------------------------------------------------------------------
56 // 56 //
57 57
58 #include "G4EmCorrections.hh" 58 #include "G4EmCorrections.hh"
>> 59 #include "Randomize.hh"
59 #include "G4PhysicalConstants.hh" 60 #include "G4PhysicalConstants.hh"
60 #include "G4SystemOfUnits.hh" 61 #include "G4SystemOfUnits.hh"
61 #include "G4ParticleTable.hh" 62 #include "G4ParticleTable.hh"
62 #include "G4IonTable.hh" 63 #include "G4IonTable.hh"
63 #include "G4VEmModel.hh" 64 #include "G4VEmModel.hh"
64 #include "G4Proton.hh" 65 #include "G4Proton.hh"
65 #include "G4GenericIon.hh" 66 #include "G4GenericIon.hh"
66 #include "G4PhysicsLogVector.hh" 67 #include "G4PhysicsLogVector.hh"
67 #include "G4ProductionCutsTable.hh" 68 #include "G4ProductionCutsTable.hh"
68 #include "G4MaterialCutsCouple.hh" 69 #include "G4MaterialCutsCouple.hh"
69 #include "G4AtomicShells.hh" 70 #include "G4AtomicShells.hh"
70 #include "G4Log.hh" 71 #include "G4Log.hh"
71 #include "G4Exp.hh" 72 #include "G4Exp.hh"
72 #include "G4Pow.hh" 73 #include "G4Pow.hh"
>> 74 #include "G4Threading.hh"
73 75
74 //....oooOO0OOooo........oooOO0OOooo........oo 76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75 77
76 namespace << 78 const G4double inveplus = 1.0/CLHEP::eplus;
77 { << 79 const G4double alpha2 = CLHEP::fine_structure_const*CLHEP::fine_structure_const;
78 constexpr G4double inveplus = 1.0/CLHEP::epl << 80
79 constexpr G4double alpha2 = CLHEP::fine_stru <<
80 } <<
81 const G4double G4EmCorrections::ZD[11] = 81 const G4double G4EmCorrections::ZD[11] =
82 {0., 0., 0., 1.72, 2.09, 2.48, 2.82, 3.16, 82 {0., 0., 0., 1.72, 2.09, 2.48, 2.82, 3.16, 3.53, 3.84, 4.15};
83 const G4double G4EmCorrections::UK[20] = {1.99 83 const G4double G4EmCorrections::UK[20] = {1.9999, 2.0134, 2.0258, 2.0478, 2.0662,
84 2.0817, 2.0945, 2.0 84 2.0817, 2.0945, 2.0999, 2.1049, 2.1132,
85 2.1197, 2.1246, 2.1 85 2.1197, 2.1246, 2.1280, 2.1292, 2.1301,
86 2.1310, 2.1310, 2.1 86 2.1310, 2.1310, 2.1300, 2.1283, 2.1271};
87 const G4double G4EmCorrections::VK[20] = {8.34 87 const G4double G4EmCorrections::VK[20] = {8.3410, 8.3373, 8.3340, 8.3287, 8.3247,
88 8.3219, 8.3201, 8.3 88 8.3219, 8.3201, 8.3194, 8.3191, 8.3188,
89 8.3191, 8.3199, 8.3 89 8.3191, 8.3199, 8.3211, 8.3218, 8.3226,
90 8.3244, 8.3264, 8.3 90 8.3244, 8.3264, 8.3285, 8.3308, 8.3320};
91 G4double G4EmCorrections::ZK[] = {0.0}; 91 G4double G4EmCorrections::ZK[] = {0.0};
92 const G4double G4EmCorrections::Eta[29] = {0.0 92 const G4double G4EmCorrections::Eta[29] = {0.005, 0.007, 0.01, 0.015, 0.02,
93 0.03, 0.04, 0.05 93 0.03, 0.04, 0.05, 0.06, 0.08,
94 0.1, 0.15, 0.2, 94 0.1, 0.15, 0.2, 0.3, 0.4,
95 0.5, 0.6, 0.7, 95 0.5, 0.6, 0.7, 0.8, 1.0,
96 1.2, 1.4, 1.5, 96 1.2, 1.4, 1.5, 1.7, 2.0, 3.0, 5.0, 7.0, 10.};
97 G4double G4EmCorrections::CK[20][29]; 97 G4double G4EmCorrections::CK[20][29];
98 G4double G4EmCorrections::CL[26][28]; 98 G4double G4EmCorrections::CL[26][28];
99 const G4double G4EmCorrections::UL[] = {0.1215 99 const G4double G4EmCorrections::UL[] = {0.1215, 0.5265, 0.8411, 1.0878, 1.2828,
100 1.4379, 1.5032, 1.5 100 1.4379, 1.5032, 1.5617, 1.6608, 1.7401,
101 1.8036, 1.8543, 1.8 101 1.8036, 1.8543, 1.8756, 1.8945, 1.9262,
102 1.9508, 1.9696, 1.9 102 1.9508, 1.9696, 1.9836, 1.9890, 1.9935,
103 2.0001, 2.0039, 2.0 103 2.0001, 2.0039, 2.0053, 2.0049, 2.0040, 2.0028};
104 G4double G4EmCorrections::VL[] = {0.0}; 104 G4double G4EmCorrections::VL[] = {0.0};
105 G4double G4EmCorrections::sWmaxBarkas = 10.0; <<
106 105
107 G4PhysicsFreeVector* G4EmCorrections::sBarkasC 106 G4PhysicsFreeVector* G4EmCorrections::sBarkasCorr = nullptr;
108 G4PhysicsFreeVector* G4EmCorrections::sThetaK 107 G4PhysicsFreeVector* G4EmCorrections::sThetaK = nullptr;
109 G4PhysicsFreeVector* G4EmCorrections::sThetaL 108 G4PhysicsFreeVector* G4EmCorrections::sThetaL = nullptr;
110 109
>> 110 #ifdef G4MULTITHREADED
>> 111 G4Mutex G4EmCorrections::theCorrMutex = G4MUTEX_INITIALIZER;
>> 112 #endif
>> 113
111 G4EmCorrections::G4EmCorrections(G4int verb) 114 G4EmCorrections::G4EmCorrections(G4int verb)
112 : verbose(verb) <<
113 { 115 {
114 eth = 2.0*CLHEP::MeV; << 116 verbose = verb;
>> 117 eth = 2.0*CLHEP::MeV;
115 eCorrMin = 25.*CLHEP::keV; 118 eCorrMin = 25.*CLHEP::keV;
116 eCorrMax = 1.*CLHEP::GeV; 119 eCorrMax = 1.*CLHEP::GeV;
117 120
118 ionTable = G4ParticleTable::GetParticleTable 121 ionTable = G4ParticleTable::GetParticleTable()->GetIonTable();
119 g4calc = G4Pow::GetInstance(); 122 g4calc = G4Pow::GetInstance();
120 123
>> 124 // G.S. Khandelwal Nucl. Phys. A116(1968)97 - 111.
>> 125 // "Shell corrections for K- and L- electrons
>> 126 nK = 20;
>> 127 nL = 26;
>> 128 nEtaK = 29;
>> 129 nEtaL = 28;
>> 130
121 // fill vectors 131 // fill vectors
122 if (nullptr == sBarkasCorr) { << 132 if(sBarkasCorr == nullptr) {
123 Initialise(); << 133 #ifdef G4MULTITHREADED
124 isInitializer = true; << 134 G4MUTEXLOCK(&theCorrMutex);
>> 135 if (sBarkasCorr == nullptr) {
>> 136 #endif
>> 137 Initialise();
>> 138 isMaster = true;
>> 139 #ifdef G4MULTITHREADED
>> 140 }
>> 141 G4MUTEXUNLOCK(&theCorrMutex);
>> 142 #endif
125 } 143 }
126 } 144 }
127 145
128 //....oooOO0OOooo........oooOO0OOooo........oo 146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
129 147
130 G4EmCorrections::~G4EmCorrections() 148 G4EmCorrections::~G4EmCorrections()
131 { 149 {
132 for (G4int i=0; i<nIons; ++i) { delete stopD << 150 for(G4int i=0; i<nIons; ++i) {delete stopData[i];}
133 if (isInitializer) { << 151 if(isMaster) {
134 delete sBarkasCorr; 152 delete sBarkasCorr;
135 delete sThetaK; 153 delete sThetaK;
136 delete sThetaL; 154 delete sThetaL;
137 sBarkasCorr = sThetaK = sThetaL = nullptr; 155 sBarkasCorr = sThetaK = sThetaL = nullptr;
138 } 156 }
139 } 157 }
140 158
141 void G4EmCorrections::SetupKinematics(const G4 159 void G4EmCorrections::SetupKinematics(const G4ParticleDefinition* p,
142 const G4Material* mat, 160 const G4Material* mat,
143 const G4double kineticEnergy) << 161 G4double kineticEnergy)
144 { 162 {
145 if(kineticEnergy != kinEnergy || p != partic 163 if(kineticEnergy != kinEnergy || p != particle) {
146 particle = p; 164 particle = p;
147 kinEnergy = kineticEnergy; 165 kinEnergy = kineticEnergy;
148 mass = p->GetPDGMass(); 166 mass = p->GetPDGMass();
149 tau = kineticEnergy / mass; 167 tau = kineticEnergy / mass;
150 gamma = 1.0 + tau; 168 gamma = 1.0 + tau;
151 bg2 = tau * (tau+2.0); 169 bg2 = tau * (tau+2.0);
152 beta2 = bg2/(gamma*gamma); 170 beta2 = bg2/(gamma*gamma);
153 beta = std::sqrt(beta2); 171 beta = std::sqrt(beta2);
154 ba2 = beta2/alpha2; 172 ba2 = beta2/alpha2;
155 const G4double ratio = CLHEP::electron_mas << 173 G4double ratio = CLHEP::electron_mass_c2/mass;
156 tmax = 2.0*CLHEP::electron_mass_c2*bg2 174 tmax = 2.0*CLHEP::electron_mass_c2*bg2
157 /(1. + 2.0*gamma*ratio + ratio*ratio); 175 /(1. + 2.0*gamma*ratio + ratio*ratio);
158 charge = p->GetPDGCharge()*inveplus; 176 charge = p->GetPDGCharge()*inveplus;
159 if(charge > 1.5) { charge = effCharge.Effe 177 if(charge > 1.5) { charge = effCharge.EffectiveCharge(p,mat,kinEnergy); }
160 q2 = charge*charge; 178 q2 = charge*charge;
161 } 179 }
162 if(mat != material) { 180 if(mat != material) {
163 material = mat; 181 material = mat;
164 theElementVector = material->GetElementVec 182 theElementVector = material->GetElementVector();
165 atomDensity = material->GetAtomicNumDensi 183 atomDensity = material->GetAtomicNumDensityVector();
166 numberOfElements = (G4int)material->GetNum << 184 numberOfElements = material->GetNumberOfElements();
167 } 185 }
168 } 186 }
169 187
170 //....oooOO0OOooo........oooOO0OOooo........oo 188 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
171 189
172 G4double G4EmCorrections::HighOrderCorrections 190 G4double G4EmCorrections::HighOrderCorrections(const G4ParticleDefinition* p,
173 191 const G4Material* mat,
174 << 192 G4double e, G4double)
175 { 193 {
176 // . Z^3 Barkas effect in the stopping power 194 // . Z^3 Barkas effect in the stopping power of matter for charged particles
177 // J.C Ashley and R.H.Ritchie 195 // J.C Ashley and R.H.Ritchie
178 // Physical review B Vol.5 No.7 1 April 19 196 // Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
179 // and ICRU49 report 197 // and ICRU49 report
180 // valid for kineticEnergy < 0.5 MeV 198 // valid for kineticEnergy < 0.5 MeV
181 // Other corrections from S.P.Ahlen Rev. M 199 // Other corrections from S.P.Ahlen Rev. Mod. Phys., Vol 52, No1, 1980
182 200
183 SetupKinematics(p, mat, e); 201 SetupKinematics(p, mat, e);
184 if(tau <= 0.0) { return 0.0; } 202 if(tau <= 0.0) { return 0.0; }
185 203
186 const G4double Barkas = BarkasCorrection(p, << 204 G4double Barkas = BarkasCorrection (p, mat, e);
187 const G4double Bloch = BlochCorrection(p, m << 205 G4double Bloch = BlochCorrection (p, mat, e);
188 const G4double Mott = MottCorrection(p, mat, << 206 G4double Mott = MottCorrection (p, mat, e);
189 207
190 G4double sum = 2.0*(Barkas + Bloch) + Mott; << 208 G4double sum = (2.0*(Barkas + Bloch) + Mott);
191 209
192 if(verbose > 1) { 210 if(verbose > 1) {
193 G4cout << "EmCorrections: E(MeV)= " << e/M 211 G4cout << "EmCorrections: E(MeV)= " << e/MeV << " Barkas= " << Barkas
194 << " Bloch= " << Bloch << " Mott= " 212 << " Bloch= " << Bloch << " Mott= " << Mott
195 << " Sum= " << sum << " q2= " << q2 213 << " Sum= " << sum << " q2= " << q2 << G4endl;
196 G4cout << " ShellCorrection: " << ShellCor 214 G4cout << " ShellCorrection: " << ShellCorrection(p, mat, e)
197 << " Kshell= " << KShellCorrection( 215 << " Kshell= " << KShellCorrection(p, mat, e)
198 << " Lshell= " << LShellCorrection( 216 << " Lshell= " << LShellCorrection(p, mat, e)
199 << " " << mat->GetName() << G4end 217 << " " << mat->GetName() << G4endl;
200 } 218 }
201 sum *= material->GetElectronDensity()*q2*CLH << 219 sum *= material->GetElectronDensity() * q2 * twopi_mc2_rcl2 /beta2;
202 return sum; 220 return sum;
203 } 221 }
204 222
205 //....oooOO0OOooo........oooOO0OOooo........oo 223 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
206 224
207 G4double G4EmCorrections::IonBarkasCorrection( 225 G4double G4EmCorrections::IonBarkasCorrection(const G4ParticleDefinition* p,
208 226 const G4Material* mat,
209 << 227 G4double e)
210 { 228 {
211 SetupKinematics(p, mat, e); << 229 return 2.0*BarkasCorrection(p, mat, e)*
212 return 2.0*BarkasCorrection(p, mat, e, true) << 230 material->GetElectronDensity() * q2 * twopi_mc2_rcl2 /beta2;
213 material->GetElectronDensity() * q2 * CLHE <<
214 } 231 }
215 232
216 //....oooOO0OOooo........oooOO0OOooo........oo 233 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
217 234
218 G4double G4EmCorrections::ComputeIonCorrection 235 G4double G4EmCorrections::ComputeIonCorrections(const G4ParticleDefinition* p,
219 236 const G4Material* mat,
220 << 237 G4double e)
221 { 238 {
222 // . Z^3 Barkas effect in the stopping power 239 // . Z^3 Barkas effect in the stopping power of matter for charged particles
223 // J.C Ashley and R.H.Ritchie 240 // J.C Ashley and R.H.Ritchie
224 // Physical review B Vol.5 No.7 1 April 19 241 // Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
225 // and ICRU49 report 242 // and ICRU49 report
226 // valid for kineticEnergy < 0.5 MeV 243 // valid for kineticEnergy < 0.5 MeV
227 // Other corrections from S.P.Ahlen Rev. M 244 // Other corrections from S.P.Ahlen Rev. Mod. Phys., Vol 52, No1, 1980
228 SetupKinematics(p, mat, e); 245 SetupKinematics(p, mat, e);
229 if(tau <= 0.0) { return 0.0; } 246 if(tau <= 0.0) { return 0.0; }
230 247
231 const G4double Barkas = BarkasCorrection (p, << 248 G4double Barkas = BarkasCorrection (p, mat, e);
232 const G4double Bloch = BlochCorrection (p, << 249 G4double Bloch = BlochCorrection (p, mat, e);
233 const G4double Mott = MottCorrection (p, mat << 250 G4double Mott = MottCorrection (p, mat, e);
234 251
235 G4double sum = 2.0*(Barkas*(charge - 1.0)/ch 252 G4double sum = 2.0*(Barkas*(charge - 1.0)/charge + Bloch) + Mott;
236 253
237 if(verbose > 1) { 254 if(verbose > 1) {
238 G4cout << "EmCorrections: E(MeV)= " << e/M 255 G4cout << "EmCorrections: E(MeV)= " << e/MeV << " Barkas= " << Barkas
239 << " Bloch= " << Bloch << " Mott= " 256 << " Bloch= " << Bloch << " Mott= " << Mott
240 << " Sum= " << sum << G4endl; 257 << " Sum= " << sum << G4endl;
241 } 258 }
242 sum *= material->GetElectronDensity() * q2 * << 259 sum *= material->GetElectronDensity() * q2 * twopi_mc2_rcl2 /beta2;
243 260
244 if(verbose > 1) { G4cout << " Sum= " << sum 261 if(verbose > 1) { G4cout << " Sum= " << sum << G4endl; }
245 return sum; 262 return sum;
246 } 263 }
247 264
248 //....oooOO0OOooo........oooOO0OOooo........oo 265 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
249 266
250 G4double G4EmCorrections::IonHighOrderCorrecti 267 G4double G4EmCorrections::IonHighOrderCorrections(const G4ParticleDefinition* p,
251 268 const G4MaterialCutsCouple* couple,
252 << 269 G4double e)
253 { 270 {
254 // . Z^3 Barkas effect in the stopping power 271 // . Z^3 Barkas effect in the stopping power of matter for charged particles
255 // J.C Ashley and R.H.Ritchie 272 // J.C Ashley and R.H.Ritchie
256 // Physical review B Vol.5 No.7 1 April 19 273 // Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
257 // and ICRU49 report 274 // and ICRU49 report
258 // valid for kineticEnergy < 0.5 MeV 275 // valid for kineticEnergy < 0.5 MeV
259 // Other corrections from S.P.Ahlen Rev. M 276 // Other corrections from S.P.Ahlen Rev. Mod. Phys., Vol 52, No1, 1980
260 277
261 G4double sum = 0.0; 278 G4double sum = 0.0;
262 279
263 if (nullptr != ionHEModel) { << 280 if(ionHEModel) {
264 G4int Z = G4lrint(p->GetPDGCharge()*invepl 281 G4int Z = G4lrint(p->GetPDGCharge()*inveplus);
265 Z = std::max(std::min(Z, 99), 1); << 282 if(Z >= 100) Z = 99;
>> 283 else if(Z < 1) Z = 1;
266 284
267 const G4double ethscaled = eth*p->GetPDGMa << 285 G4double ethscaled = eth*p->GetPDGMass()/proton_mass_c2;
268 const G4int ionPDG = p->GetPDGEncoding(); << 286 G4int ionPDG = p->GetPDGEncoding();
269 auto iter = thcorr.find(ionPDG); << 287 if(thcorr.find(ionPDG)==thcorr.end()) { // Not found: fill the map
270 if (iter == thcorr.end()) { // Not found: <<
271 std::vector<G4double> v; 288 std::vector<G4double> v;
272 for(std::size_t i=0; i<ncouples; ++i){ << 289 for(size_t i=0; i<ncouples; ++i){
273 v.push_back(ethscaled*ComputeIonCorrec 290 v.push_back(ethscaled*ComputeIonCorrections(p,currmat[i],ethscaled));
274 } 291 }
275 thcorr.insert(std::pair< G4int, std::vec 292 thcorr.insert(std::pair< G4int, std::vector<G4double> >(ionPDG,v));
276 } 293 }
277 G4double rest = 0.0; << 294
278 iter = thcorr.find(ionPDG); << 295 //G4cout << " map size=" << thcorr.size() << G4endl;
279 if (iter != thcorr.end()) { rest = (iter-> << 296 //for(std::map< G4int, std::vector<G4double> >::iterator
>> 297 // it = thcorr.begin(); it != thcorr.end(); ++it){
>> 298 // G4cout << "\t map element: first (key)=" << it->first
>> 299 // << "\t second (vector): vec size=" << (it->second).size() << G4endl;
>> 300 // for(size_t i=0; i<(it->second).size(); ++i){
>> 301 // G4cout << "\t \t vec element: [" << i << "]=" << (it->second)[i]
>> 302 //<< G4endl; } }
>> 303
>> 304 G4double rest = (thcorr.find(ionPDG)->second)[couple->GetIndex()];
280 305
281 sum = ComputeIonCorrections(p,couple->GetM 306 sum = ComputeIonCorrections(p,couple->GetMaterial(),e) - rest/e;
282 307
283 if(verbose > 1) { 308 if(verbose > 1) {
284 G4cout << " Sum= " << sum << " dSum= " < 309 G4cout << " Sum= " << sum << " dSum= " << rest/e << G4endl;
285 } 310 }
286 } 311 }
287 return sum; 312 return sum;
288 } 313 }
289 314
290 //....oooOO0OOooo........oooOO0OOooo........oo 315 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
291 316
292 G4double G4EmCorrections::Bethe(const G4Partic 317 G4double G4EmCorrections::Bethe(const G4ParticleDefinition* p,
293 const G4Materi 318 const G4Material* mat,
294 const G4double << 319 G4double e)
295 { 320 {
296 SetupKinematics(p, mat, e); 321 SetupKinematics(p, mat, e);
297 const G4double eexc = material->GetIonisati 322 const G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
298 const G4double eexc2 = eexc*eexc; 323 const G4double eexc2 = eexc*eexc;
299 return 0.5*G4Log(2.0*electron_mass_c2*bg2*tm 324 return 0.5*G4Log(2.0*electron_mass_c2*bg2*tmax/eexc2)-beta2;
300 } 325 }
301 326
302 //....oooOO0OOooo........oooOO0OOooo........oo 327 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
303 328
304 G4double G4EmCorrections::SpinCorrection(const 329 G4double G4EmCorrections::SpinCorrection(const G4ParticleDefinition* p,
305 const 330 const G4Material* mat,
306 const << 331 G4double e)
307 { 332 {
308 SetupKinematics(p, mat, e); 333 SetupKinematics(p, mat, e);
309 const G4double dedx = 0.5*tmax/(kinEnergy + 334 const G4double dedx = 0.5*tmax/(kinEnergy + mass);
310 return 0.5*dedx*dedx; 335 return 0.5*dedx*dedx;
311 } 336 }
312 337
313 //....oooOO0OOooo........oooOO0OOooo........oo 338 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
314 339
315 G4double G4EmCorrections:: KShellCorrection(co 340 G4double G4EmCorrections:: KShellCorrection(const G4ParticleDefinition* p,
316 co 341 const G4Material* mat,
317 co << 342 G4double e)
318 { 343 {
319 SetupKinematics(p, mat, e); 344 SetupKinematics(p, mat, e);
320 G4double term = 0.0; 345 G4double term = 0.0;
321 for (G4int i = 0; i<numberOfElements; ++i) { 346 for (G4int i = 0; i<numberOfElements; ++i) {
322 347
323 G4double Z = (*theElementVector)[i]->GetZ( 348 G4double Z = (*theElementVector)[i]->GetZ();
324 G4int iz = (*theElementVector)[i]->GetZa 349 G4int iz = (*theElementVector)[i]->GetZasInt();
325 G4double f = 1.0; 350 G4double f = 1.0;
326 G4double Z2= (Z-0.3)*(Z-0.3); 351 G4double Z2= (Z-0.3)*(Z-0.3);
327 if(1 == iz) { 352 if(1 == iz) {
328 f = 0.5; 353 f = 0.5;
329 Z2 = 1.0; 354 Z2 = 1.0;
330 } 355 }
331 const G4double eta = ba2/Z2; << 356 G4double eta = ba2/Z2;
332 const G4double tet = (11 < iz) ? sThetaK-> << 357 G4double tet = Z2*(1. + Z2*0.25*alpha2);
>> 358 if(11 < iz) { tet = sThetaK->Value(Z); }
333 term += f*atomDensity[i]*KShell(tet,eta)/Z 359 term += f*atomDensity[i]*KShell(tet,eta)/Z;
334 } 360 }
335 361
336 term /= material->GetTotNbOfAtomsPerVolume() 362 term /= material->GetTotNbOfAtomsPerVolume();
337 363
338 return term; 364 return term;
339 } 365 }
340 366
341 //....oooOO0OOooo........oooOO0OOooo........oo 367 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
342 368
343 G4double G4EmCorrections:: LShellCorrection(co 369 G4double G4EmCorrections:: LShellCorrection(const G4ParticleDefinition* p,
344 co 370 const G4Material* mat,
345 co << 371 G4double e)
346 { 372 {
347 SetupKinematics(p, mat, e); 373 SetupKinematics(p, mat, e);
348 G4double term = 0.0; 374 G4double term = 0.0;
349 for (G4int i = 0; i<numberOfElements; ++i) { 375 for (G4int i = 0; i<numberOfElements; ++i) {
350 376
351 const G4double Z = (*theElementVector)[i]- << 377 G4double Z = (*theElementVector)[i]->GetZ();
352 const G4int iz = (*theElementVector)[i]->G << 378 G4int iz = (*theElementVector)[i]->GetZasInt();
353 if(2 < iz) { 379 if(2 < iz) {
354 const G4double Zeff = (iz < 10) ? Z - ZD << 380 G4double Zeff = Z - ZD[10];
355 const G4double Z2= Zeff*Zeff; << 381 if(iz < 10) { Zeff = Z - ZD[iz]; }
356 const G4double eta = ba2/Z2; << 382 G4double Z2= Zeff*Zeff;
>> 383 G4double f = 0.125;
>> 384 G4double eta = ba2/Z2;
357 G4double tet = sThetaL->Value(Z); 385 G4double tet = sThetaL->Value(Z);
358 G4int nmax = std::min(4, G4AtomicShells: << 386 G4int nmax = std::min(4,G4AtomicShells::GetNumberOfShells(iz));
359 for (G4int j=1; j<nmax; ++j) { << 387 for(G4int j=1; j<nmax; ++j) {
360 G4int ne = G4AtomicShells::GetNumberOf 388 G4int ne = G4AtomicShells::GetNumberOfElectrons(iz,j);
361 if (15 >= iz) { << 389 if(15 >= iz) {
362 tet = (3 > j) ? 0.25*Z2*(1.0 + 5*Z2* << 390 if(3 > j) { tet = 0.25*Z2*(1.0 + 5*Z2*alpha2/16.); }
363 0.25*Z2*(1.0 + Z2*alpha2/16.); << 391 else { tet = 0.25*Z2*(1.0 + Z2*alpha2/16.); }
364 } 392 }
365 //G4cout << " LShell: j= " << j << " n 393 //G4cout << " LShell: j= " << j << " ne= " << ne << " e(eV)= " << e/eV
366 // << " ThetaL= " << tet << G4en 394 // << " ThetaL= " << tet << G4endl;
367 term += 0.125*ne*atomDensity[i]*LShell << 395 term += f*ne*atomDensity[i]*LShell(tet,eta)/Z;
368 } 396 }
369 } 397 }
370 } 398 }
371 399
372 term /= material->GetTotNbOfAtomsPerVolume() 400 term /= material->GetTotNbOfAtomsPerVolume();
373 401
374 return term; 402 return term;
375 } 403 }
376 404
377 //....oooOO0OOooo........oooOO0OOooo........oo 405 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
378 406
379 G4double G4EmCorrections::KShell(const G4doubl << 407 G4double G4EmCorrections::KShell(G4double tet, G4double eta)
380 { 408 {
381 G4double corr = 0.0; 409 G4double corr = 0.0;
382 410
383 static const G4double TheK[20] = 411 static const G4double TheK[20] =
384 {0.64, 0.65, 0.66, 0.68, 0.70, 0.72, 0.74, 412 {0.64, 0.65, 0.66, 0.68, 0.70, 0.72, 0.74, 0.75, 0.76, 0.78,
385 0.80, 0.82, 0.84, 0.85, 0.86, 0.88, 0.90, 413 0.80, 0.82, 0.84, 0.85, 0.86, 0.88, 0.90, 0.92, 0.94, 0.95};
386 414
387 415
388 G4double x = tet; 416 G4double x = tet;
389 G4int itet = 0; 417 G4int itet = 0;
390 G4int ieta = 0; 418 G4int ieta = 0;
391 if(tet < TheK[0]) { 419 if(tet < TheK[0]) {
392 x = TheK[0]; 420 x = TheK[0];
393 } else if(tet > TheK[nK-1]) { 421 } else if(tet > TheK[nK-1]) {
394 x = TheK[nK-1]; 422 x = TheK[nK-1];
395 itet = nK-2; 423 itet = nK-2;
396 } else { 424 } else {
397 itet = Index(x, TheK, nK); 425 itet = Index(x, TheK, nK);
398 } 426 }
399 // asymptotic case << 427 // assimptotic case
400 if(eta >= Eta[nEtaK-1]) { 428 if(eta >= Eta[nEtaK-1]) {
401 corr = 429 corr =
402 (Value(x, TheK[itet], TheK[itet+1], UK[i 430 (Value(x, TheK[itet], TheK[itet+1], UK[itet], UK[itet+1]) +
403 Value(x, TheK[itet], TheK[itet+1], VK[i 431 Value(x, TheK[itet], TheK[itet+1], VK[itet], VK[itet+1])/eta +
404 Value(x, TheK[itet], TheK[itet+1], ZK[i 432 Value(x, TheK[itet], TheK[itet+1], ZK[itet], ZK[itet+1])/(eta*eta))/eta;
405 } else { 433 } else {
406 G4double y = eta; 434 G4double y = eta;
407 if(eta < Eta[0]) { 435 if(eta < Eta[0]) {
408 y = Eta[0]; << 436 y = Eta[0];
409 } else { 437 } else {
410 ieta = Index(y, Eta, nEtaK); 438 ieta = Index(y, Eta, nEtaK);
411 } 439 }
412 corr = Value2(x, y, TheK[itet], TheK[itet+ 440 corr = Value2(x, y, TheK[itet], TheK[itet+1], Eta[ieta], Eta[ieta+1],
413 CK[itet][ieta], CK[itet+1][i 441 CK[itet][ieta], CK[itet+1][ieta],
414 CK[itet][ieta+1], CK[itet+1] 442 CK[itet][ieta+1], CK[itet+1][ieta+1]);
415 //G4cout << " x= " <<x<<" y= "<<y<<" tet 443 //G4cout << " x= " <<x<<" y= "<<y<<" tet= " <<TheK[itet]
416 // <<" "<< TheK[itet+1]<<" eta= 444 // <<" "<< TheK[itet+1]<<" eta= "<< Eta[ieta]<<" "<< Eta[ieta+1]
417 // <<" CK= " << CK[itet][ieta]<< 445 // <<" CK= " << CK[itet][ieta]<<" "<< CK[itet+1][ieta]
418 // <<" "<< CK[itet][ieta+1]<<" " 446 // <<" "<< CK[itet][ieta+1]<<" "<< CK[itet+1][ieta+1]<<G4endl;
419 } 447 }
420 //G4cout << "Kshell: tet= " << tet << " eta 448 //G4cout << "Kshell: tet= " << tet << " eta= " << eta << " C= " << corr
421 // << " itet= " << itet << " ieta= 449 // << " itet= " << itet << " ieta= " << ieta <<G4endl;
422 return corr; 450 return corr;
423 } 451 }
424 452
425 //....oooOO0OOooo........oooOO0OOooo........oo 453 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
426 454
427 G4double G4EmCorrections::LShell(const G4doubl << 455 G4double G4EmCorrections::LShell(G4double tet, G4double eta)
428 { 456 {
429 G4double corr = 0.0; 457 G4double corr = 0.0;
430 458
431 static const G4double TheL[26] = 459 static const G4double TheL[26] =
432 {0.24, 0.26, 0.28, 0.30, 0.32, 0.34, 0.3 460 {0.24, 0.26, 0.28, 0.30, 0.32, 0.34, 0.35, 0.36, 0.38, 0.40,
433 0.42, 0.44, 0.45, 0.46, 0.48, 0.50, 0.5 461 0.42, 0.44, 0.45, 0.46, 0.48, 0.50, 0.52, 0.54, 0.55, 0.56,
434 0.58, 0.60, 0.62, 0.64, 0.65, 0.66}; 462 0.58, 0.60, 0.62, 0.64, 0.65, 0.66};
435 463
436 G4double x = tet; 464 G4double x = tet;
437 G4int itet = 0; 465 G4int itet = 0;
438 G4int ieta = 0; 466 G4int ieta = 0;
439 if(tet < TheL[0]) { 467 if(tet < TheL[0]) {
440 x = TheL[0]; 468 x = TheL[0];
441 } else if(tet > TheL[nL-1]) { 469 } else if(tet > TheL[nL-1]) {
442 x = TheL[nL-1]; 470 x = TheL[nL-1];
443 itet = nL-2; 471 itet = nL-2;
444 } else { 472 } else {
445 itet = Index(x, TheL, nL); 473 itet = Index(x, TheL, nL);
446 } 474 }
447 475
448 // asymptotic case << 476 // assimptotic case
449 if(eta >= Eta[nEtaL-1]) { 477 if(eta >= Eta[nEtaL-1]) {
450 corr = (Value(x, TheL[itet], TheL[itet+1], 478 corr = (Value(x, TheL[itet], TheL[itet+1], UL[itet], UL[itet+1])
451 + Value(x, TheL[itet], TheL[itet 479 + Value(x, TheL[itet], TheL[itet+1], VL[itet], VL[itet+1])/eta
452 )/eta; 480 )/eta;
453 } else { 481 } else {
454 G4double y = eta; 482 G4double y = eta;
455 if(eta < Eta[0]) { 483 if(eta < Eta[0]) {
456 y = Eta[0]; 484 y = Eta[0];
457 } else { 485 } else {
458 ieta = Index(y, Eta, nEtaL); 486 ieta = Index(y, Eta, nEtaL);
459 } 487 }
460 corr = Value2(x, y, TheL[itet], TheL[itet+ 488 corr = Value2(x, y, TheL[itet], TheL[itet+1], Eta[ieta], Eta[ieta+1],
461 CL[itet][ieta], CL[itet+1][i 489 CL[itet][ieta], CL[itet+1][ieta],
462 CL[itet][ieta+1], CL[itet+1] 490 CL[itet][ieta+1], CL[itet+1][ieta+1]);
463 //G4cout << " x= " <<x<<" y= "<<y<<" tet 491 //G4cout << " x= " <<x<<" y= "<<y<<" tet= " <<TheL[itet]
464 // <<" "<< TheL[itet+1]<<" eta= 492 // <<" "<< TheL[itet+1]<<" eta= "<< Eta[ieta]<<" "<< Eta[ieta+1]
465 // <<" CL= " << CL[itet][ieta]<< 493 // <<" CL= " << CL[itet][ieta]<<" "<< CL[itet+1][ieta]
466 // <<" "<< CL[itet][ieta+1]<<" " 494 // <<" "<< CL[itet][ieta+1]<<" "<< CL[itet+1][ieta+1]<<G4endl;
467 } 495 }
468 //G4cout<<"Lshell: tet= "<<tet<<" eta= "<<e 496 //G4cout<<"Lshell: tet= "<<tet<<" eta= "<<eta<<" itet= "<<itet
469 // <<" ieta= "<<ieta<<" Corr= "<<cor 497 // <<" ieta= "<<ieta<<" Corr= "<<corr<<G4endl;
470 return corr; 498 return corr;
471 } 499 }
472 500
473 //....oooOO0OOooo........oooOO0OOooo........oo 501 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
474 502
475 G4double G4EmCorrections::ShellCorrectionSTD(c 503 G4double G4EmCorrections::ShellCorrectionSTD(const G4ParticleDefinition* p,
476 c 504 const G4Material* mat,
477 c << 505 G4double e)
478 { 506 {
479 SetupKinematics(p, mat, e); 507 SetupKinematics(p, mat, e);
480 G4double taulim= 8.0*MeV/mass; 508 G4double taulim= 8.0*MeV/mass;
481 G4double bg2lim= taulim * (taulim+2.0); 509 G4double bg2lim= taulim * (taulim+2.0);
482 510
483 G4double* shellCorrectionVector = 511 G4double* shellCorrectionVector =
484 material->GetIonisation()->GetShel 512 material->GetIonisation()->GetShellCorrectionVector();
485 G4double sh = 0.0; 513 G4double sh = 0.0;
486 G4double x = 1.0; 514 G4double x = 1.0;
487 G4double taul = material->GetIonisation()-> 515 G4double taul = material->GetIonisation()->GetTaul();
488 516
489 if ( bg2 >= bg2lim ) { 517 if ( bg2 >= bg2lim ) {
490 for (G4int k=0; k<3; ++k) { 518 for (G4int k=0; k<3; ++k) {
491 x *= bg2 ; 519 x *= bg2 ;
492 sh += shellCorrectionVector[k]/x; 520 sh += shellCorrectionVector[k]/x;
493 } 521 }
494 522
495 } else { 523 } else {
496 for (G4int k=0; k<3; ++k) { 524 for (G4int k=0; k<3; ++k) {
497 x *= bg2lim ; 525 x *= bg2lim ;
498 sh += shellCorrectionVector[k]/x; 526 sh += shellCorrectionVector[k]/x;
499 } 527 }
500 sh *= G4Log(tau/taul)/G4Log(taulim/taul); 528 sh *= G4Log(tau/taul)/G4Log(taulim/taul);
501 } 529 }
502 sh *= 0.5; 530 sh *= 0.5;
503 return sh; 531 return sh;
504 } 532 }
505 533
506 534
507 //....oooOO0OOooo........oooOO0OOooo........oo 535 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
508 536
509 G4double G4EmCorrections::ShellCorrection(cons 537 G4double G4EmCorrections::ShellCorrection(const G4ParticleDefinition* p,
510 cons 538 const G4Material* mat,
511 cons << 539 G4double ekin)
512 { 540 {
513 SetupKinematics(p, mat, ekin); 541 SetupKinematics(p, mat, ekin);
>> 542
514 G4double term = 0.0; 543 G4double term = 0.0;
515 //G4cout << "### G4EmCorrections::ShellCorre 544 //G4cout << "### G4EmCorrections::ShellCorrection " << mat->GetName()
516 // << " " << ekin/MeV << " MeV " << << 545 // << " " << ekin/MeV << " MeV " << G4endl;
517 for (G4int i = 0; i<numberOfElements; ++i) { 546 for (G4int i = 0; i<numberOfElements; ++i) {
518 547
519 G4double res = 0.0; 548 G4double res = 0.0;
520 G4double res0 = 0.0; 549 G4double res0 = 0.0;
521 const G4double Z = (*theElementVector)[i]- << 550 G4double Z = (*theElementVector)[i]->GetZ();
522 const G4int iz = (*theElementVector)[i]->G << 551 G4int iz = (*theElementVector)[i]->GetZasInt();
523 G4double Z2= (Z-0.3)*(Z-0.3); 552 G4double Z2= (Z-0.3)*(Z-0.3);
524 G4double f = 1.0; 553 G4double f = 1.0;
525 if(1 == iz) { 554 if(1 == iz) {
526 f = 0.5; 555 f = 0.5;
527 Z2 = 1.0; 556 Z2 = 1.0;
528 } 557 }
529 G4double eta = ba2/Z2; 558 G4double eta = ba2/Z2;
530 G4double tet = (11 < iz) ? sThetaK->Value( << 559 G4double tet = Z2*(1. + Z2*0.25*alpha2);
>> 560 if(11 < iz) { tet = sThetaK->Value(Z); }
531 res0 = f*KShell(tet,eta); 561 res0 = f*KShell(tet,eta);
532 res += res0; 562 res += res0;
533 //G4cout << " Z= " << iz << " Shell 0" << 563 //G4cout << " Z= " << iz << " Shell 0" << " tet= " << tet
534 // << " eta= " << eta << " resK= " 564 // << " eta= " << eta << " resK= " << res0 << G4endl;
535 <<
536 if(2 < iz) { 565 if(2 < iz) {
537 const G4double Zeff = (iz < 10) ? Z - ZD << 566 G4double Zeff = Z - ZD[10];
>> 567 if(iz < 10) { Zeff = Z - ZD[iz]; }
538 Z2= Zeff*Zeff; 568 Z2= Zeff*Zeff;
539 eta = ba2/Z2; 569 eta = ba2/Z2;
540 tet = sThetaL->Value(Z); <<
541 f = 0.125; 570 f = 0.125;
542 const G4int ntot = G4AtomicShells::GetNu << 571 tet = sThetaL->Value(Z);
543 const G4int nmax = std::min(4, ntot); << 572 G4int ntot = G4AtomicShells::GetNumberOfShells(iz);
>> 573 G4int nmax = std::min(4, ntot);
544 G4double norm = 0.0; 574 G4double norm = 0.0;
545 G4double eshell = 0.0; 575 G4double eshell = 0.0;
546 for(G4int j=1; j<nmax; ++j) { 576 for(G4int j=1; j<nmax; ++j) {
547 G4int ne = G4AtomicShells::GetNumberOf 577 G4int ne = G4AtomicShells::GetNumberOfElectrons(iz,j);
548 if(15 >= iz) { 578 if(15 >= iz) {
549 tet = (3 > j) ? 0.25*Z2*(1.0 + 5*Z2* << 579 if(3 > j) { tet = 0.25*Z2*(1.0 + 5*Z2*alpha2/16.); }
550 0.25*Z2*(1.0 + Z2*alpha2/16.); << 580 else { tet = 0.25*Z2*(1.0 + Z2*alpha2/16.); }
551 } 581 }
552 norm += ne; 582 norm += ne;
553 eshell += tet*ne; 583 eshell += tet*ne;
554 res0 = f*ne*LShell(tet,eta); 584 res0 = f*ne*LShell(tet,eta);
555 res += res0; 585 res += res0;
556 //G4cout << " Zeff= " << Zeff << " She << 586 //G4cout << " Z= " << iz << " Shell " << j << " Ne= " << ne
557 // << " tet= " << tet << " eta= 587 // << " tet= " << tet << " eta= " << eta
558 // << " resL= " << res0 << G4en 588 // << " resL= " << res0 << G4endl;
559 } 589 }
560 if(ntot > nmax) { 590 if(ntot > nmax) {
561 if (norm > 0.0) { norm = 1.0/norm; } << 591 eshell /= norm;
562 eshell *= norm; <<
563 592
564 static const G4double HM[53] = { 593 static const G4double HM[53] = {
565 12.0, 12.0, 12.0, 12.0, 11.9, 11.7, 11.5, 594 12.0, 12.0, 12.0, 12.0, 11.9, 11.7, 11.5, 11.2, 10.8, 10.4,
566 10.0, 9.51, 8.97, 8.52, 8.03, 7.46, 6.95, 595 10.0, 9.51, 8.97, 8.52, 8.03, 7.46, 6.95, 6.53, 6.18, 5.87,
567 5.61, 5.39, 5.19, 5.01, 4.86, 4.72, 4.62, 596 5.61, 5.39, 5.19, 5.01, 4.86, 4.72, 4.62, 4.53, 4.44, 4.38,
568 4.32, 4.26, 4.20, 4.15, 4.1, 4.04, 4.00, 597 4.32, 4.26, 4.20, 4.15, 4.1, 4.04, 4.00, 3.95, 3.93, 3.91,
569 3.90, 3.89, 3.89, 3.88, 3.88, 3.88, 3.88, 598 3.90, 3.89, 3.89, 3.88, 3.88, 3.88, 3.88, 3.88, 3.89, 3.89,
570 3.90, 3.92, 3.93 }; 599 3.90, 3.92, 3.93 };
571 static const G4double HN[31] = { 600 static const G4double HN[31] = {
572 75.5, 61.9, 52.2, 45.1, 39.6, 35.4, 31.9, 601 75.5, 61.9, 52.2, 45.1, 39.6, 35.4, 31.9, 29.1, 27.2, 25.8,
573 24.5, 23.6, 22.7, 22.0, 21.4, 20.9, 20.5, 602 24.5, 23.6, 22.7, 22.0, 21.4, 20.9, 20.5, 20.2, 19.9, 19.7,
574 19.5, 19.3, 19.2, 19.1, 18.4, 18.8, 18.7, 603 19.5, 19.3, 19.2, 19.1, 18.4, 18.8, 18.7, 18.6, 18.5, 18.4,
575 18.2}; 604 18.2};
576 605
577 // Add M-shell 606 // Add M-shell
578 if(28 > iz) { 607 if(28 > iz) {
579 res += f*(iz - 10)*LShell(eshell,HM[ 608 res += f*(iz - 10)*LShell(eshell,HM[iz-11]*eta);
580 } else if(63 > iz) { 609 } else if(63 > iz) {
581 res += f*18*LShell(eshell,HM[iz-11]* 610 res += f*18*LShell(eshell,HM[iz-11]*eta);
582 } else { 611 } else {
583 res += f*18*LShell(eshell,HM[52]*eta 612 res += f*18*LShell(eshell,HM[52]*eta);
584 } 613 }
585 // Add N-shell 614 // Add N-shell
586 if(32 < iz) { 615 if(32 < iz) {
587 if(60 > iz) { 616 if(60 > iz) {
588 res += f*(iz - 28)*LShell(eshell,H 617 res += f*(iz - 28)*LShell(eshell,HN[iz-33]*eta);
589 } else if(63 > iz) { 618 } else if(63 > iz) {
590 res += 4*LShell(eshell,HN[iz-33]*e 619 res += 4*LShell(eshell,HN[iz-33]*eta);
591 } else { 620 } else {
592 res += 4*LShell(eshell,HN[30]*eta) 621 res += 4*LShell(eshell,HN[30]*eta);
593 } 622 }
594 // Add O-P-shells 623 // Add O-P-shells
595 if(60 < iz) { 624 if(60 < iz) {
596 res += f*(iz - 60)*LShell(eshell,1 625 res += f*(iz - 60)*LShell(eshell,150*eta);
597 } 626 }
598 } 627 }
599 } 628 }
600 } 629 }
601 term += res*atomDensity[i]/Z; 630 term += res*atomDensity[i]/Z;
602 } 631 }
603 632
604 term /= material->GetTotNbOfAtomsPerVolume() 633 term /= material->GetTotNbOfAtomsPerVolume();
605 //G4cout << "##Shell Correction=" << term << << 634 //G4cout << "# Shell Correction= " << term << G4endl;
606 return term; 635 return term;
607 } 636 }
608 637
609 //....oooOO0OOooo........oooOO0OOooo........oo 638 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
610 639
611 G4double G4EmCorrections::DensityCorrection(co 640 G4double G4EmCorrections::DensityCorrection(const G4ParticleDefinition* p,
612 co 641 const G4Material* mat,
613 co << 642 G4double e)
614 { 643 {
615 SetupKinematics(p, mat, e); 644 SetupKinematics(p, mat, e);
616 645
617 G4double cden = material->GetIonisation()-> 646 G4double cden = material->GetIonisation()->GetCdensity();
618 G4double mden = material->GetIonisation()-> 647 G4double mden = material->GetIonisation()->GetMdensity();
619 G4double aden = material->GetIonisation()-> 648 G4double aden = material->GetIonisation()->GetAdensity();
620 G4double x0den = material->GetIonisation()-> 649 G4double x0den = material->GetIonisation()->GetX0density();
621 G4double x1den = material->GetIonisation()-> 650 G4double x1den = material->GetIonisation()->GetX1density();
622 651
623 G4double dedx = 0.0; 652 G4double dedx = 0.0;
624 653
625 // density correction 654 // density correction
626 static const G4double twoln10 = 2.0*G4Log(10 655 static const G4double twoln10 = 2.0*G4Log(10.0);
627 G4double x = G4Log(bg2)/twoln10; 656 G4double x = G4Log(bg2)/twoln10;
628 if ( x >= x0den ) { 657 if ( x >= x0den ) {
629 dedx = twoln10*x - cden ; 658 dedx = twoln10*x - cden ;
630 if ( x < x1den ) dedx += aden*G4Exp(G4Log( 659 if ( x < x1den ) dedx += aden*G4Exp(G4Log(x1den-x)*mden) ;
631 } 660 }
632 661
633 return dedx; 662 return dedx;
634 } 663 }
635 664
636 //....oooOO0OOooo........oooOO0OOooo........oo 665 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
637 666
638 G4double G4EmCorrections::BarkasCorrection(con 667 G4double G4EmCorrections::BarkasCorrection(const G4ParticleDefinition* p,
639 con 668 const G4Material* mat,
640 con << 669 G4double e)
641 con <<
642 { 670 {
643 // . Z^3 Barkas effect in the stopping power 671 // . Z^3 Barkas effect in the stopping power of matter for charged particles
644 // J.C Ashley and R.H.Ritchie 672 // J.C Ashley and R.H.Ritchie
645 // Physical review B Vol.5 No.7 1 April 19 673 // Physical review B Vol.5 No.7 1 April 1972 pp. 2393-2397
646 // valid for kineticEnergy > 0.5 MeV 674 // valid for kineticEnergy > 0.5 MeV
647 675
648 if (!isInitialized) { SetupKinematics(p, mat << 676 SetupKinematics(p, mat, e);
649 G4double BarkasTerm = 0.0; 677 G4double BarkasTerm = 0.0;
650 678
651 for (G4int i = 0; i<numberOfElements; ++i) { 679 for (G4int i = 0; i<numberOfElements; ++i) {
652 680
653 const G4int iz = (*theElementVector)[i]->G << 681 G4double Z = (*theElementVector)[i]->GetZ();
>> 682 G4int iz = (*theElementVector)[i]->GetZasInt();
654 if(iz == 47) { 683 if(iz == 47) {
655 BarkasTerm += atomDensity[i]*0.006812*G4 684 BarkasTerm += atomDensity[i]*0.006812*G4Exp(-G4Log(beta)*0.9);
656 } else if(iz >= 64) { 685 } else if(iz >= 64) {
657 BarkasTerm += atomDensity[i]*0.002833*G4 686 BarkasTerm += atomDensity[i]*0.002833*G4Exp(-G4Log(beta)*1.2);
658 } else { 687 } else {
659 688
660 const G4double Z = (*theElementVector)[i << 689 G4double X = ba2 / Z;
661 const G4double X = ba2 / Z; <<
662 G4double b = 1.3; 690 G4double b = 1.3;
663 if(1 == iz) { b = (material->GetName() = << 691 if(1 == iz) {
>> 692 if(material->GetName() == "G4_lH2") { b = 0.6; }
>> 693 else { b = 1.8; }
>> 694 }
664 else if(2 == iz) { b = 0.6; } 695 else if(2 == iz) { b = 0.6; }
665 else if(10 >= iz) { b = 1.8; } 696 else if(10 >= iz) { b = 1.8; }
666 else if(17 >= iz) { b = 1.4; } 697 else if(17 >= iz) { b = 1.4; }
667 else if(18 == iz) { b = 1.8; } 698 else if(18 == iz) { b = 1.8; }
668 else if(25 >= iz) { b = 1.4; } 699 else if(25 >= iz) { b = 1.4; }
669 else if(50 >= iz) { b = 1.35;} 700 else if(50 >= iz) { b = 1.35;}
670 701
671 const G4double W = b/std::sqrt(X); << 702 G4double W = b/std::sqrt(X);
672 703
673 G4double val = sBarkasCorr->Value(W, idx << 704 G4double val = sBarkasCorr->Value(W);
674 if (W > sWmaxBarkas) { val *= (sWmaxBark << 705 if(W > sBarkasCorr->Energy(46)) {
>> 706 val *= sBarkasCorr->Energy(46)/W;
>> 707 }
675 // G4cout << "i= " << i << " b= " << 708 // G4cout << "i= " << i << " b= " << b << " W= " << W
676 // << " Z= " << Z << " X= " << X << " va 709 // << " Z= " << Z << " X= " << X << " val= " << val<< G4endl;
677 BarkasTerm += val*atomDensity[i] / (std: 710 BarkasTerm += val*atomDensity[i] / (std::sqrt(Z*X)*X);
678 } 711 }
679 } 712 }
680 713
681 BarkasTerm *= 1.29*charge/material->GetTotNb 714 BarkasTerm *= 1.29*charge/material->GetTotNbOfAtomsPerVolume();
682 715
683 return BarkasTerm; 716 return BarkasTerm;
684 } 717 }
685 718
686 //....oooOO0OOooo........oooOO0OOooo........oo 719 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
687 720
688 G4double G4EmCorrections::BlochCorrection(cons 721 G4double G4EmCorrections::BlochCorrection(const G4ParticleDefinition* p,
689 cons 722 const G4Material* mat,
690 cons << 723 G4double e)
691 cons <<
692 { 724 {
693 if (!isInitialized) { SetupKinematics(p, mat << 725 SetupKinematics(p, mat, e);
694 726
695 G4double y2 = q2/ba2; 727 G4double y2 = q2/ba2;
696 728
697 G4double term = 1.0/(1.0 + y2); 729 G4double term = 1.0/(1.0 + y2);
698 G4double del; 730 G4double del;
699 G4double j = 1.0; 731 G4double j = 1.0;
700 do { 732 do {
701 j += 1.0; 733 j += 1.0;
702 del = 1.0/(j* (j*j + y2)); 734 del = 1.0/(j* (j*j + y2));
703 term += del; 735 term += del;
704 // Loop checking, 03-Aug-2015, Vladimir Iv 736 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
705 } while (del > 0.01*term); 737 } while (del > 0.01*term);
706 738
707 return -y2*term; << 739 G4double res = -y2*term;
>> 740 return res;
708 } 741 }
709 742
710 //....oooOO0OOooo........oooOO0OOooo........oo 743 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
711 744
712 G4double G4EmCorrections::MottCorrection(const 745 G4double G4EmCorrections::MottCorrection(const G4ParticleDefinition* p,
713 const 746 const G4Material* mat,
714 const << 747 G4double e)
715 const <<
716 { 748 {
717 if (!isInitialized) { SetupKinematics(p, mat << 749 SetupKinematics(p, mat, e);
718 return CLHEP::pi*CLHEP::fine_structure_const << 750 G4double mterm = CLHEP::pi*fine_structure_const*beta*charge;
>> 751 return mterm;
719 } 752 }
720 753
721 //....oooOO0OOooo........oooOO0OOooo........oo 754 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
722 755
723 G4double 756 G4double
724 G4EmCorrections::EffectiveChargeCorrection(con 757 G4EmCorrections::EffectiveChargeCorrection(const G4ParticleDefinition* p,
725 con 758 const G4Material* mat,
726 con << 759 G4double ekin)
727 { 760 {
728 G4double factor = 1.0; 761 G4double factor = 1.0;
729 if(p->GetPDGCharge() <= 2.5*CLHEP::eplus || 762 if(p->GetPDGCharge() <= 2.5*CLHEP::eplus || nIons <= 0) { return factor; }
730 763
731 if(verbose > 1) { 764 if(verbose > 1) {
732 G4cout << "EffectiveChargeCorrection: " << 765 G4cout << "EffectiveChargeCorrection: " << p->GetParticleName()
733 << " in " << mat->GetName() 766 << " in " << mat->GetName()
734 << " ekin(MeV)= " << ekin/MeV << G4 767 << " ekin(MeV)= " << ekin/MeV << G4endl;
735 } 768 }
736 769
737 if(p != curParticle || mat != curMaterial) { 770 if(p != curParticle || mat != curMaterial) {
738 curParticle = p; 771 curParticle = p;
739 curMaterial = mat; 772 curMaterial = mat;
740 curVector = nullptr; 773 curVector = nullptr;
741 currentZ = p->GetAtomicNumber(); 774 currentZ = p->GetAtomicNumber();
742 if(verbose > 1) { 775 if(verbose > 1) {
743 G4cout << "G4EmCorrections::EffectiveCha 776 G4cout << "G4EmCorrections::EffectiveChargeCorrection: Zion= "
744 << currentZ << " Aion= " << p->Ge 777 << currentZ << " Aion= " << p->GetPDGMass()/amu_c2 << G4endl;
745 } 778 }
746 massFactor = CLHEP::proton_mass_c2/p->GetP 779 massFactor = CLHEP::proton_mass_c2/p->GetPDGMass();
747 idx = -1; 780 idx = -1;
748 781
749 for(G4int i=0; i<nIons; ++i) { 782 for(G4int i=0; i<nIons; ++i) {
750 if(materialList[i] == mat && currentZ == 783 if(materialList[i] == mat && currentZ == Zion[i]) {
751 idx = i; 784 idx = i;
752 break; 785 break;
753 } 786 }
754 } 787 }
755 //G4cout << " idx= " << idx << " dz= " << 788 //G4cout << " idx= " << idx << " dz= " << G4endl;
756 if(idx >= 0) { 789 if(idx >= 0) {
757 if(nullptr == ionList[idx]) { BuildCorre 790 if(nullptr == ionList[idx]) { BuildCorrectionVector(); }
758 curVector = stopData[idx]; 791 curVector = stopData[idx];
759 } else { 792 } else {
760 return factor; 793 return factor;
761 } 794 }
762 } 795 }
763 if(nullptr != curVector) { 796 if(nullptr != curVector) {
764 factor = curVector->Value(ekin*massFactor) 797 factor = curVector->Value(ekin*massFactor);
765 if(verbose > 1) { 798 if(verbose > 1) {
766 G4cout << "E= " << ekin << " factor= " < 799 G4cout << "E= " << ekin << " factor= " << factor << " massfactor= "
767 << massFactor << G4endl; 800 << massFactor << G4endl;
768 } 801 }
769 } 802 }
770 return factor; 803 return factor;
771 } 804 }
772 805
773 //....oooOO0OOooo........oooOO0OOooo........oo 806 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
774 807
775 void G4EmCorrections::AddStoppingData(const G4 << 808 void G4EmCorrections::AddStoppingData(G4int Z, G4int A,
776 const G4 809 const G4String& mname,
777 G4Physic 810 G4PhysicsVector* dVector)
778 { 811 {
779 G4int i = 0; 812 G4int i = 0;
780 for(; i<nIons; ++i) { 813 for(; i<nIons; ++i) {
781 if(Z == Zion[i] && A == Aion[i] && mname = 814 if(Z == Zion[i] && A == Aion[i] && mname == materialName[i]) break;
782 } 815 }
783 if(i == nIons) { 816 if(i == nIons) {
784 Zion.push_back(Z); 817 Zion.push_back(Z);
785 Aion.push_back(A); 818 Aion.push_back(A);
786 materialName.push_back(mname); 819 materialName.push_back(mname);
787 materialList.push_back(nullptr); 820 materialList.push_back(nullptr);
788 ionList.push_back(nullptr); 821 ionList.push_back(nullptr);
789 stopData.push_back(dVector); 822 stopData.push_back(dVector);
790 nIons++; 823 nIons++;
791 if(verbose > 1) { 824 if(verbose > 1) {
792 G4cout << "AddStoppingData Z= " << Z << 825 G4cout << "AddStoppingData Z= " << Z << " A= " << A << " " << mname
793 << " idx= " << i << G4endl; 826 << " idx= " << i << G4endl;
794 } 827 }
795 } 828 }
796 } 829 }
797 830
798 //....oooOO0OOooo........oooOO0OOooo........oo 831 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
799 832
800 void G4EmCorrections::BuildCorrectionVector() 833 void G4EmCorrections::BuildCorrectionVector()
801 { 834 {
802 if(nullptr == ionLEModel || nullptr == ionHE 835 if(nullptr == ionLEModel || nullptr == ionHEModel) {
803 return; 836 return;
804 } 837 }
805 838
806 const G4ParticleDefinition* ion = curParticl 839 const G4ParticleDefinition* ion = curParticle;
807 const G4ParticleDefinition* gion = G4Generic 840 const G4ParticleDefinition* gion = G4GenericIon::GenericIon();
808 G4int Z = Zion[idx]; 841 G4int Z = Zion[idx];
809 G4double A = Aion[idx]; 842 G4double A = Aion[idx];
810 G4PhysicsVector* v = stopData[idx]; 843 G4PhysicsVector* v = stopData[idx];
811 844
812 if(verbose > 1) { 845 if(verbose > 1) {
813 G4cout << "BuildCorrectionVector: Stopping 846 G4cout << "BuildCorrectionVector: Stopping for "
814 << curParticle->GetParticleName() < 847 << curParticle->GetParticleName() << " in "
815 << materialName[idx] << " Ion Z= " 848 << materialName[idx] << " Ion Z= " << Z << " A= " << A
816 << " massFactor= " << massFactor << 849 << " massFactor= " << massFactor << G4endl;
817 G4cout << " Nbins=" << nbinCorr << " Em 850 G4cout << " Nbins=" << nbinCorr << " Emin(MeV)=" << eCorrMin
818 << " Emax(MeV)=" << eCorrMax << " ion: " 851 << " Emax(MeV)=" << eCorrMax << " ion: "
819 << ion->GetParticleName() << G4endl 852 << ion->GetParticleName() << G4endl;
820 } 853 }
821 854
822 auto vv = new G4PhysicsLogVector(eCorrMin,eC << 855 G4PhysicsLogVector* vv =
>> 856 new G4PhysicsLogVector(eCorrMin,eCorrMax,nbinCorr,false);
823 const G4double eth0 = v->Energy(0); 857 const G4double eth0 = v->Energy(0);
824 const G4double escal = eth/massFactor; 858 const G4double escal = eth/massFactor;
825 G4double qe = 859 G4double qe =
826 effCharge.EffectiveChargeSquareRatio(curPa 860 effCharge.EffectiveChargeSquareRatio(curParticle, curMaterial, escal);
827 const G4double dedxt = 861 const G4double dedxt =
828 ionLEModel->ComputeDEDXPerVolume(curMateri 862 ionLEModel->ComputeDEDXPerVolume(curMaterial, gion, eth, eth)*qe;
829 const G4double dedx1t = 863 const G4double dedx1t =
830 ionHEModel->ComputeDEDXPerVolume(curMateri 864 ionHEModel->ComputeDEDXPerVolume(curMaterial, gion, eth, eth)*qe
831 + ComputeIonCorrections(curParticle, curMa 865 + ComputeIonCorrections(curParticle, curMaterial, escal);
832 const G4double rest = escal*(dedxt - dedx1t) 866 const G4double rest = escal*(dedxt - dedx1t);
833 if(verbose > 1) { 867 if(verbose > 1) {
834 G4cout << "Escal(MeV)= " << escal << " qe= 868 G4cout << "Escal(MeV)= " << escal << " qe=" << qe
835 << " dedxt= " << dedxt << " dedx1t= 869 << " dedxt= " << dedxt << " dedx1t= " << dedx1t << G4endl;
836 } 870 }
837 for(G4int i=0; i<=nbinCorr; ++i) { 871 for(G4int i=0; i<=nbinCorr; ++i) {
838 // energy in the local table (GenericIon) 872 // energy in the local table (GenericIon)
839 G4double e = vv->Energy(i); 873 G4double e = vv->Energy(i);
840 // energy of the real ion 874 // energy of the real ion
841 G4double eion = e/massFactor; 875 G4double eion = e/massFactor;
842 // energy in the imput stopping data vecto 876 // energy in the imput stopping data vector
843 G4double e1 = eion/A; 877 G4double e1 = eion/A;
844 G4double dedx = (e1 < eth0) 878 G4double dedx = (e1 < eth0)
845 ? v->Value(eth0)*std::sqrt(e1/eth0) : v- 879 ? v->Value(eth0)*std::sqrt(e1/eth0) : v->Value(e1);
846 qe = effCharge.EffectiveChargeSquareRatio( 880 qe = effCharge.EffectiveChargeSquareRatio(curParticle, curMaterial, eion);
847 G4double dedx1 = (e <= eth) 881 G4double dedx1 = (e <= eth)
848 ? ionLEModel->ComputeDEDXPerVolume(curMa 882 ? ionLEModel->ComputeDEDXPerVolume(curMaterial, gion, e, e)*qe
849 : ionHEModel->ComputeDEDXPerVolume(curMa 883 : ionHEModel->ComputeDEDXPerVolume(curMaterial, gion, e, e)*qe +
850 ComputeIonCorrections(curParticle, curMa 884 ComputeIonCorrections(curParticle, curMaterial, eion) + rest/eion;
851 vv->PutValue(i, dedx/dedx1); 885 vv->PutValue(i, dedx/dedx1);
852 if(verbose > 1) { 886 if(verbose > 1) {
853 G4cout << "E(MeV)=" << e/CLHEP::MeV << " 887 G4cout << "E(MeV)=" << e/CLHEP::MeV << " Eion=" << eion/CLHEP::MeV
854 << " e1=" << e1 << " dedxRatio= " << de 888 << " e1=" << e1 << " dedxRatio= " << dedx/dedx1
855 << " dedx=" << dedx << " dedx1=" 889 << " dedx=" << dedx << " dedx1=" << dedx1
856 << " qe=" << qe << " rest/eion=" 890 << " qe=" << qe << " rest/eion=" << rest/eion << G4endl;
857 } 891 }
858 } 892 }
859 delete v; 893 delete v;
860 ionList[idx] = ion; 894 ionList[idx] = ion;
861 stopData[idx] = vv; 895 stopData[idx] = vv;
862 if(verbose > 1) { G4cout << "End data set " 896 if(verbose > 1) { G4cout << "End data set " << G4endl; }
863 } 897 }
864 898
865 //....oooOO0OOooo........oooOO0OOooo........oo 899 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
866 900
867 void G4EmCorrections::InitialiseForNewRun() 901 void G4EmCorrections::InitialiseForNewRun()
868 { 902 {
869 G4ProductionCutsTable* tb = G4ProductionCuts 903 G4ProductionCutsTable* tb = G4ProductionCutsTable::GetProductionCutsTable();
870 ncouples = tb->GetTableSize(); 904 ncouples = tb->GetTableSize();
871 if(currmat.size() != ncouples) { 905 if(currmat.size() != ncouples) {
872 currmat.resize(ncouples); 906 currmat.resize(ncouples);
873 for(auto it = thcorr.begin(); it != thcorr << 907 for(std::map< G4int, std::vector<G4double> >::iterator it =
>> 908 thcorr.begin(); it != thcorr.end(); ++it){
874 (it->second).clear(); 909 (it->second).clear();
875 } 910 }
876 thcorr.clear(); 911 thcorr.clear();
877 for(std::size_t i=0; i<ncouples; ++i) { << 912 for(size_t i=0; i<ncouples; ++i) {
878 currmat[i] = tb->GetMaterialCutsCouple(( << 913 currmat[i] = tb->GetMaterialCutsCouple(i)->GetMaterial();
879 G4String nam = currmat[i]->GetName(); 914 G4String nam = currmat[i]->GetName();
880 for(G4int j=0; j<nIons; ++j) { 915 for(G4int j=0; j<nIons; ++j) {
881 if(nam == materialName[j]) { materialL 916 if(nam == materialName[j]) { materialList[j] = currmat[i]; }
882 } 917 }
883 } 918 }
884 } 919 }
885 } 920 }
886 921
887 //....oooOO0OOooo........oooOO0OOooo........oo 922 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
888 923
889 void G4EmCorrections::Initialise() 924 void G4EmCorrections::Initialise()
890 { 925 {
891 // Z^3 Barkas effect in the stopping power o 926 // Z^3 Barkas effect in the stopping power of matter for charged particles
892 // J.C Ashley and R.H.Ritchie 927 // J.C Ashley and R.H.Ritchie
893 // Physical review B Vol.5 No.7 1 April 1972 928 // Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
894 // G.S. Khandelwal Nucl. Phys. A116(1968)97 <<
895 // "Shell corrections for K- and L- electron <<
896 <<
897 G4int i, j; 929 G4int i, j;
898 static const G4double fTable[47][2] = { 930 static const G4double fTable[47][2] = {
899 { 0.02, 21.5}, 931 { 0.02, 21.5},
900 { 0.03, 20.0}, 932 { 0.03, 20.0},
901 { 0.04, 18.0}, 933 { 0.04, 18.0},
902 { 0.05, 15.6}, 934 { 0.05, 15.6},
903 { 0.06, 15.0}, 935 { 0.06, 15.0},
904 { 0.07, 14.0}, 936 { 0.07, 14.0},
905 { 0.08, 13.5}, 937 { 0.08, 13.5},
906 { 0.09, 13.}, 938 { 0.09, 13.},
907 { 0.1, 12.2}, 939 { 0.1, 12.2},
908 { 0.2, 9.25}, 940 { 0.2, 9.25},
909 { 0.3, 7.0}, 941 { 0.3, 7.0},
910 { 0.4, 6.0}, 942 { 0.4, 6.0},
911 { 0.5, 4.5}, 943 { 0.5, 4.5},
912 { 0.6, 3.5}, 944 { 0.6, 3.5},
913 { 0.7, 3.0}, 945 { 0.7, 3.0},
914 { 0.8, 2.5}, 946 { 0.8, 2.5},
915 { 0.9, 2.0}, 947 { 0.9, 2.0},
916 { 1.0, 1.7}, 948 { 1.0, 1.7},
917 { 1.2, 1.2}, 949 { 1.2, 1.2},
918 { 1.3, 1.0}, 950 { 1.3, 1.0},
919 { 1.4, 0.86}, 951 { 1.4, 0.86},
920 { 1.5, 0.7}, 952 { 1.5, 0.7},
921 { 1.6, 0.61}, 953 { 1.6, 0.61},
922 { 1.7, 0.52}, 954 { 1.7, 0.52},
923 { 1.8, 0.5}, 955 { 1.8, 0.5},
924 { 1.9, 0.43}, 956 { 1.9, 0.43},
925 { 2.0, 0.42}, 957 { 2.0, 0.42},
926 { 2.1, 0.3}, 958 { 2.1, 0.3},
927 { 2.4, 0.2}, 959 { 2.4, 0.2},
928 { 3.0, 0.13}, 960 { 3.0, 0.13},
929 { 3.08, 0.1}, 961 { 3.08, 0.1},
930 { 3.1, 0.09}, 962 { 3.1, 0.09},
931 { 3.3, 0.08}, 963 { 3.3, 0.08},
932 { 3.5, 0.07}, 964 { 3.5, 0.07},
933 { 3.8, 0.06}, 965 { 3.8, 0.06},
934 { 4.0, 0.051}, 966 { 4.0, 0.051},
935 { 4.1, 0.04}, 967 { 4.1, 0.04},
936 { 4.8, 0.03}, 968 { 4.8, 0.03},
937 { 5.0, 0.024}, 969 { 5.0, 0.024},
938 { 5.1, 0.02}, 970 { 5.1, 0.02},
939 { 6.0, 0.013}, 971 { 6.0, 0.013},
940 { 6.5, 0.01}, 972 { 6.5, 0.01},
941 { 7.0, 0.009}, 973 { 7.0, 0.009},
942 { 7.1, 0.008}, 974 { 7.1, 0.008},
943 { 8.0, 0.006}, 975 { 8.0, 0.006},
944 { 9.0, 0.0032}, 976 { 9.0, 0.0032},
945 { 10.0, 0.0025} }; 977 { 10.0, 0.0025} };
946 978
947 sBarkasCorr = new G4PhysicsFreeVector(47, fa 979 sBarkasCorr = new G4PhysicsFreeVector(47, false);
948 for(i=0; i<47; ++i) { sBarkasCorr->PutValues 980 for(i=0; i<47; ++i) { sBarkasCorr->PutValues(i, fTable[i][0], fTable[i][1]); }
949 981
950 const G4double SK[20] = {1.9477, 1.9232, 1.8 << 982 static const G4double SK[20] = {1.9477, 1.9232, 1.8996, 1.8550, 1.8137,
951 1.7754, 1.7396, 1.7 983 1.7754, 1.7396, 1.7223, 1.7063, 1.6752,
952 1.6461, 1.6189, 1.5 984 1.6461, 1.6189, 1.5933, 1.5811, 1.5693,
953 1.5467, 1.5254, 1.5 985 1.5467, 1.5254, 1.5053, 1.4863, 1.4772};
954 const G4double TK[20] = {2.5222, 2.5125, 2.5 << 986 static const G4double TK[20] = {2.5222, 2.5125, 2.5026, 2.4821, 2.4608,
955 2.4388, 2.4163, 2.4 987 2.4388, 2.4163, 2.4044, 2.3933, 2.3701,
956 2.3466, 2.3229, 2.2 988 2.3466, 2.3229, 2.2992, 2.2872, 2.2753,
957 2.2515, 2.2277, 2.2 989 2.2515, 2.2277, 2.2040, 2.1804, 2.1686};
958 990
959 const G4double SL[26] = {15.3343, 13.9389, 1 << 991 static const G4double SL[26] = {15.3343, 13.9389, 12.7909, 11.8343, 11.0283,
960 10.3424, 10.0371, 992 10.3424, 10.0371, 9.7537, 9.2443, 8.8005,
961 8.4114, 8.0683, 993 8.4114, 8.0683, 7.9117, 7.7641, 7.4931,
962 7.2506, 7.0327, 994 7.2506, 7.0327, 6.8362, 6.7452, 6.6584,
963 6.4969, 6.3498, 995 6.4969, 6.3498, 6.2154, 6.0923, 6.0345, 5.9792};
964 const G4double TL[26] = {35.0669, 33.4344, 3 << 996 static const G4double TL[26] = {35.0669, 33.4344, 32.0073, 30.7466, 29.6226,
965 28.6128, 28.1449, 2 997 28.6128, 28.1449, 27.6991, 26.8674, 26.1061,
966 25.4058, 24.7587, 2 998 25.4058, 24.7587, 24.4531, 24.1583, 23.5992,
967 23.0771, 22.5880, 2 999 23.0771, 22.5880, 22.1285, 21.9090, 21.6958,
968 21.2872, 20.9006, 2 1000 21.2872, 20.9006, 20.5341, 20.1859, 20.0183, 19.8546};
969 1001
970 const G4double bk1[29][11] = { << 1002 static const G4double bk1[29][11] = {
971 {0.005, 1.34782E-8, 1.46132E-8, 1.72179E-8, 1003 {0.005, 1.34782E-8, 1.46132E-8, 1.72179E-8, 2.03521E-8, 2.41370E-8, 2.87247E-8, 3.13778E-8, 3.43072E-8, 4.11274E-8, 4.94946E-8},
972 {0.007, 6.87555E-8, 7.44373E-8, 8.74397E-8, 1004 {0.007, 6.87555E-8, 7.44373E-8, 8.74397E-8, 1.03022E-7, 1.21760E-7, 1.44370E-7, 1.57398E-7, 1.71747E-7, 2.05023E-7, 2.45620E-7},
973 {0.01, 3.78413E-7, 4.08831E-7, 4.78154E-7, 5 1005 {0.01, 3.78413E-7, 4.08831E-7, 4.78154E-7, 5.60760E-7, 6.59478E-7, 7.77847E-7, 8.45709E-7, 9.20187E-7, 1.09192E-6, 1.29981E-6},
974 {0.015, 2.53200E-6, 2.72664E-6, 3.16738E-6, 1006 {0.015, 2.53200E-6, 2.72664E-6, 3.16738E-6, 3.68791E-6, 4.30423E-6, 5.03578E-6, 5.45200E-6, 5.90633E-6, 6.94501E-6, 8.18757E-6},
975 {0.02, 9.43891E-6, 1.01339E-5, 1.16984E-5, 1 1007 {0.02, 9.43891E-6, 1.01339E-5, 1.16984E-5, 1.35313E-5, 1.56829E-5, 1.82138E-5, 1.96439E-5, 2.11973E-5, 2.47216E-5, 2.88935E-5},
976 {0.03, 5.67088E-5, 6.05576E-5, 6.91311E-5, 7 1008 {0.03, 5.67088E-5, 6.05576E-5, 6.91311E-5, 7.90324E-5, 9.04832E-5, 1.03744E-4, 1.11147E-4, 1.19122E-4, 1.36980E-4, 1.57744E-4},
977 {0.04, 1.91576E-4, 2.03626E-4, 2.30230E-4, 2 1009 {0.04, 1.91576E-4, 2.03626E-4, 2.30230E-4, 2.60584E-4, 2.95248E-4, 3.34870E-4, 3.56771E-4, 3.80200E-4, 4.32104E-4, 4.91584E-4},
978 {0.05, 4.74226E-4, 5.02006E-4, 5.62872E-4, 6 1010 {0.05, 4.74226E-4, 5.02006E-4, 5.62872E-4, 6.31597E-4, 7.09244E-4, 7.97023E-4, 8.45134E-4, 8.96410E-4, 1.00867E-3, 1.13590E-3},
979 {0.06, 9.67285E-4, 1.02029E-3, 1.13566E-3, 1 1011 {0.06, 9.67285E-4, 1.02029E-3, 1.13566E-3, 1.26476E-3, 1.46928E-3, 1.57113E-3, 1.65921E-3, 1.75244E-3, 1.95562E-3, 2.18336E-3},
980 {0.08, 2.81537E-3, 2.95200E-3, 3.24599E-3, 3 1012 {0.08, 2.81537E-3, 2.95200E-3, 3.24599E-3, 3.57027E-3, 3.92793E-3, 4.32246E-3, 4.53473E-3, 4.75768E-3, 5.23785E-3, 5.76765E-3},
981 {0.1, 6.14216E-3, 6.40864E-3, 6.97750E-3, 7. 1013 {0.1, 6.14216E-3, 6.40864E-3, 6.97750E-3, 7.59781E-3, 8.27424E-3, 9.01187E-3, 9.40534E-3, 9.81623E-3, 1.06934E-2, 1.16498E-2},
982 {0.15, 2.23096E-2, 2.30790E-2, 2.46978E-2, 2 1014 {0.15, 2.23096E-2, 2.30790E-2, 2.46978E-2, 2.64300E-2, 2.82832E-2, 3.02661E-2, 3.13090E-2, 3.23878E-2, 3.46580E-2, 3.70873E-2},
983 {0.2, 5.04022E-2, 5.18492E-2, 5.48682E-2, 5. 1015 {0.2, 5.04022E-2, 5.18492E-2, 5.48682E-2, 5.80617E-2, 6.14403E-2, 6.50152E-2, 6.68798E-2, 6.87981E-2, 7.28020E-2, 7.70407E-2},
984 {0.3, 1.38018E-1, 1.41026E-1, 1.47244E-1, 1. 1016 {0.3, 1.38018E-1, 1.41026E-1, 1.47244E-1, 1.53743E-1, 1.60536E-1, 1.67641E-1, 1.71315E-1, 1.75074E-1, 1.82852E-1, 1.90997E-1},
985 {0.4, 2.56001E-1, 2.60576E-1, 2.69992E-1, 2. 1017 {0.4, 2.56001E-1, 2.60576E-1, 2.69992E-1, 2.79776E-1, 2.89946E-1, 3.00525E-1, 3.05974E-1, 3.11533E-1, 3.22994E-1, 3.34935E-1},
986 {0.5, 3.92181E-1, 3.98213E-1, 4.10597E-1, 4. 1018 {0.5, 3.92181E-1, 3.98213E-1, 4.10597E-1, 4.23427E-1, 4.36726E-1, 4.50519E-1, 4.57610E-1, 4.64834E-1, 4.79700E-1, 4.95148E-1},
987 {0.6, 5.37913E-1, 5.45268E-1, 5.60350E-1, 5. 1019 {0.6, 5.37913E-1, 5.45268E-1, 5.60350E-1, 5.75948E-1, 5.92092E-1, 6.08811E-1, 6.17396E-1, 6.26136E-1, 6.44104E-1, 6.62752E-1},
988 {0.7, 6.87470E-1, 6.96021E-1, 7.13543E-1, 7. 1020 {0.7, 6.87470E-1, 6.96021E-1, 7.13543E-1, 7.31650E-1, 7.50373E-1, 7.69748E-1, 7.79591E-1, 7.89811E-1, 8.10602E-1, 8.32167E-1},
989 {0.8, 8.37159E-1, 8.46790E-1, 8.66525E-1, 8. 1021 {0.8, 8.37159E-1, 8.46790E-1, 8.66525E-1, 8.86910E-1, 9.07979E-1, 9.29772E-1, 9.40953E-1, 9.52331E-1, 9.75701E-1, 9.99934E-1},
990 {1.0, 1.12850, 1.14002, 1.16362, 1.18799, 1. 1022 {1.0, 1.12850, 1.14002, 1.16362, 1.18799, 1.21317, 1.23921, 1.25257, 1.26616, 1.29408, 1.32303},
991 {1.2, 1.40232, 1.41545, 1.44232, 1.47007, 1. 1023 {1.2, 1.40232, 1.41545, 1.44232, 1.47007, 1.49875, 1.52842, 1.54364, 1.55913, 1.59095, 1.62396},
992 {1.4, 1.65584, 1.67034, 1.70004, 1.73072, 1. 1024 {1.4, 1.65584, 1.67034, 1.70004, 1.73072, 1.76244, 1.79526, 1.81210, 1.82925, 1.86448, 1.90104},
993 {1.5, 1.77496, 1.79009, 1.82107, 1.85307, 1. 1025 {1.5, 1.77496, 1.79009, 1.82107, 1.85307, 1.88617, 1.92042, 1.93800, 1.95590, 1.99269, 2.03087},
994 {1.7, 1.99863, 2.01490, 2.04822, 2.08265, 2. 1026 {1.7, 1.99863, 2.01490, 2.04822, 2.08265, 2.11827, 2.15555, 2.17409, 2.19337, 2.23302, 2.27419},
995 {2.0, 2.29325, 2.31100, 2.34738, 2.38501, 2. 1027 {2.0, 2.29325, 2.31100, 2.34738, 2.38501, 2.42395, 2.46429, 2.48401, 2.50612, 2.54955, 2.59468},
996 {3.0, 3.08887, 3.11036, 3.15445, 3.20013, 3. 1028 {3.0, 3.08887, 3.11036, 3.15445, 3.20013, 3.24748, 3.29664, 3.32192, 3.34770, 3.40081, 3.45611},
997 {5.0, 4.07599, 4.10219, 4.15606, 4.21199, 4. 1029 {5.0, 4.07599, 4.10219, 4.15606, 4.21199, 4.27010, 4.33056, 4.36172, 4.39353, 4.45918, 4.52772},
998 {7.0, 4.69647, 4.72577, 4.78608, 4.84877, 4. 1030 {7.0, 4.69647, 4.72577, 4.78608, 4.84877, 4.91402, 4.98200, 5.01707, 5.05290, 5.12695, 5.20436},
999 {10.0, 5.32590, 5.35848, 5.42560, 5.49547, 5 1031 {10.0, 5.32590, 5.35848, 5.42560, 5.49547, 5.56830, 5.64429, 5.68353, 5.72366, 5.80666, 5.89359}
1000 }; 1032 };
1001 1033
1002 const G4double bk2[29][11] = { << 1034 static const G4double bk2[29][11] = {
1003 {0.005, 5.98040E-8, 7.25636E-8, 8.00602E-8, 1035 {0.005, 5.98040E-8, 7.25636E-8, 8.00602E-8, 8.84294E-8, 1.08253E-7, 1.33148E-7, 1.64573E-7, 2.04459E-7, 2.28346E-7, 2.55370E-7},
1004 {0.007, 2.95345E-7, 3.56497E-7, 3.92247E-7, 1036 {0.007, 2.95345E-7, 3.56497E-7, 3.92247E-7, 4.32017E-7, 5.25688E-7, 6.42391E-7, 7.88464E-7, 9.72171E-7, 1.08140E-6, 1.20435E-6},
1005 {0.01, 1.55232E-6, 1.86011E-6, 2.03881E-6, 1037 {0.01, 1.55232E-6, 1.86011E-6, 2.03881E-6, 2.23662E-6, 2.69889E-6, 3.26860E-6, 3.26860E-6, 4.84882E-6, 5.36428E-6, 5.94048E-6},
1006 {0.015, 9.67802E-6, 1.14707E-5, 1.25008E-5, 1038 {0.015, 9.67802E-6, 1.14707E-5, 1.25008E-5, 1.36329E-5, 1.62480E-5, 1.94200E-5, 2.32783E-5, 2.79850E-5, 3.07181E-5, 3.37432E-5},
1007 {0.02, 3.38425E-5, 3.97259E-5, 4.30763E-5, 1039 {0.02, 3.38425E-5, 3.97259E-5, 4.30763E-5, 4.67351E-5, 5.51033E-5, 6.51154E-5, 7.71154E-5, 9.15431E-5, 9.98212E-5, 1.08909E-4},
1008 {0.03, 1.81920E-4, 2.10106E-4, 2.25918E-4, 1040 {0.03, 1.81920E-4, 2.10106E-4, 2.25918E-4, 2.43007E-4, 2.81460E-4, 3.26458E-4, 3.79175E-4, 4.41006E-4, 4.75845E-4, 5.13606E-4},
1009 {0.04, 5.59802E-4, 6.38103E-4, 6.81511E-4, 1041 {0.04, 5.59802E-4, 6.38103E-4, 6.81511E-4, 7.28042E-4, 8.31425E-4, 9.50341E-4, 1.08721E-3, 1.24485E-3, 1.33245E-3, 1.42650E-3},
1010 {0.05, 1.28002E-3, 1.44336E-3, 1.53305E-3, 1042 {0.05, 1.28002E-3, 1.44336E-3, 1.53305E-3, 1.62855E-3, 1.83861E-3, 2.07396E-3, 2.34750E-3, 2.65469E-3, 2.82358E-3, 3.00358E-3},
1011 {0.06, 2.42872E-3, 2.72510E-3, 2.88111E-3, 1043 {0.06, 2.42872E-3, 2.72510E-3, 2.88111E-3, 3.04636E-3, 3.40681E-3, 3.81132E-3, 4.26536E-3, 4.77507E-3, 5.05294E-3, 5.34739E-3},
1012 {0.08, 6.35222E-3, 6.99730E-3, 7.34446E-3, 1044 {0.08, 6.35222E-3, 6.99730E-3, 7.34446E-3, 7.70916E-3, 8.49478E-3, 9.36187E-3, 1.03189E-2, 1.13754E-2, 1.19441E-2, 1.25417E-2},
1013 {0.1, 1.26929E-2, 1.38803E-2, 1.44371E-2, 1 1045 {0.1, 1.26929E-2, 1.38803E-2, 1.44371E-2, 1.50707E-2, 1.64235E-2, 1.78989E-2, 1.95083E-2, 2.12640E-2, 2.22009E-2, 2.31795E-2},
1014 {0.15, 3.96872E-2, 4.24699E-2, 4.39340E-2, 1046 {0.15, 3.96872E-2, 4.24699E-2, 4.39340E-2, 4.54488E-2, 4.86383E-2, 5.20542E-2, 5.57135E-2, 5.96350E-2, 6.17003E-2, 6.38389E-2},
1015 {0.2, 8.15290E-2, 8.62830E-2, 8.87650E-2, 9 1047 {0.2, 8.15290E-2, 8.62830E-2, 8.87650E-2, 9.13200E-2, 9.66589E-2, 1.02320E-1, 1.08326E-1, 1.14701E-1, 1.18035E-1, 1.21472E-1},
1016 {0.3, 1.99528E-1, 2.08471E-1, 2.13103E-1, 2 1048 {0.3, 1.99528E-1, 2.08471E-1, 2.13103E-1, 2.17848E-1, 2.27689E-1, 2.38022E-1, 2.48882E-1, 2.60304E-1, 2.66239E-1, 2.72329E-1},
1017 {0.4, 3.47383E-1, 3.60369E-1, 3.67073E-1, 3 1049 {0.4, 3.47383E-1, 3.60369E-1, 3.67073E-1, 3.73925E-1, 3.88089E-1, 4.02900E-1, 4.18404E-1, 4.34647E-1, 4.43063E-1, 4.51685E-1},
1018 {0.5, 5.11214E-1, 5.27935E-1, 5.36554E-1, 5 1050 {0.5, 5.11214E-1, 5.27935E-1, 5.36554E-1, 5.45354E-1, 5.63515E-1, 5.82470E-1, 6.02273E-1, 6.22986E-1, 6.33705E-1, 6.44677E-1},
1019 {0.6, 6.82122E-1, 7.02260E-1, 7.12631E-1, 7 1051 {0.6, 6.82122E-1, 7.02260E-1, 7.12631E-1, 7.23214E-1, 7.45041E-1, 7.67800E-1, 7.91559E-1, 8.16391E-1, 8.29235E-1, 8.42380E-1},
1020 {0.7, 8.54544E-1, 8.77814E-1, 8.89791E-1, 9 1052 {0.7, 8.54544E-1, 8.77814E-1, 8.89791E-1, 9.02008E-1, 9.27198E-1, 9.53454E-1, 9.80856E-1, 1.00949, 1.02430, 1.03945},
1021 {0.8, 1.02508, 1.05121, 1.06466, 1.07838, 1 1053 {0.8, 1.02508, 1.05121, 1.06466, 1.07838, 1.10667, 1.13615, 1.16692, 1.19907, 1.21570, 1.23272},
1022 {1.0, 1.35307, 1.38429, 1.40036, 1.41676, 1 1054 {1.0, 1.35307, 1.38429, 1.40036, 1.41676, 1.45057, 1.48582, 1.52263, 1.56111, 1.58102, 1.60140},
1023 {1.2, 1.65823, 1.69385, 1.71220, 1.73092, 1 1055 {1.2, 1.65823, 1.69385, 1.71220, 1.73092, 1.76954, 1.80983, 1.85192, 1.89596, 1.91876, 1.94211},
1024 {1.4, 1.93902, 1.97852, 1.99887, 2.01964, 2 1056 {1.4, 1.93902, 1.97852, 1.99887, 2.01964, 2.06251, 2.10727, 2.15406, 2.20304, 2.22842, 2.25442},
1025 {1.5, 2.07055, 2.11182, 2.13309, 2.15480, 2 1057 {1.5, 2.07055, 2.11182, 2.13309, 2.15480, 2.19963, 2.24644, 2.29539, 2.34666, 2.37323, 2.40045},
1026 {1.7, 2.31700, 2.36154, 2.38451, 2.40798, 2 1058 {1.7, 2.31700, 2.36154, 2.38451, 2.40798, 2.45641, 2.50703, 2.56000, 2.61552, 2.64430, 2.67381},
1027 {2.0, 2.64162, 2.69053, 2.71576, 2.74154, 2 1059 {2.0, 2.64162, 2.69053, 2.71576, 2.74154, 2.79481, 2.85052, 2.90887, 2.97009, 3.00185, 3.03442},
1028 {3.0, 3.51376, 3.57394, 3.60503, 3.63684, 3 1060 {3.0, 3.51376, 3.57394, 3.60503, 3.63684, 3.70268, 3.77170, 3.84418, 3.92040, 3.96003, 4.00073},
1029 {5.0, 4.59935, 4.67433, 4.71316, 4.75293, 4 1061 {5.0, 4.59935, 4.67433, 4.71316, 4.75293, 4.83543, 4.92217, 5.01353, 5.10992, 5.16014, 5.21181},
1030 {7.0, 5.28542, 5.37040, 5.41445, 5.45962, 5 1062 {7.0, 5.28542, 5.37040, 5.41445, 5.45962, 5.55344, 5.65226, 5.79496, 5.90517, 5.96269, 6.02191},
1031 {10.0, 5.98474, 6.08046, 6.13015, 6.18112, 1063 {10.0, 5.98474, 6.08046, 6.13015, 6.18112, 6.28715, 6.39903, 6.51728, 6.64249, 6.70792, 6.77535}
1032 }; 1064 };
1033 1065
1034 const G4double bls1[28][10] = { << 1066 static const G4double bls1[28][10] = {
1035 {0.005, 2.4111E-4, 2.5612E-4, 2.7202E-4, 3. 1067 {0.005, 2.4111E-4, 2.5612E-4, 2.7202E-4, 3.0658E-4, 3.4511E-4, 3.8795E-4, 4.3542E-4, 4.6100E-4, 4.8786E-4},
1036 {0.007, 6.3947E-4, 6.7058E-4, 7.0295E-4, 7. 1068 {0.007, 6.3947E-4, 6.7058E-4, 7.0295E-4, 7.7167E-4, 8.4592E-4, 9.2605E-4, 1.0125E-3, 1.0583E-3, 1.1058E-3},
1037 {0.01, 1.5469E-3, 1.6036E-3, 1.6622E-3, 1.7 1069 {0.01, 1.5469E-3, 1.6036E-3, 1.6622E-3, 1.7856E-3, 1.9181E-3, 2.1615E-3, 2.3178E-3, 2.4019E-3, 2.4904E-3},
1038 {0.015, 3.7221E-3, 3.8404E-3, 3.9650E-3, 4. 1070 {0.015, 3.7221E-3, 3.8404E-3, 3.9650E-3, 4.2354E-3, 4.5396E-3, 4.8853E-3, 5.2820E-3, 5.5031E-3, 5.7414E-3},
1039 {0.02, 6.9449E-3, 7.1910E-3, 7.4535E-3, 8.0 1071 {0.02, 6.9449E-3, 7.1910E-3, 7.4535E-3, 8.0336E-3, 8.6984E-3, 9.4638E-3, 1.0348E-2, 1.0841E-2, 1.1372E-2},
1040 {0.03, 1.7411E-2, 1.8180E-2, 1.8997E-2, 2.0 1072 {0.03, 1.7411E-2, 1.8180E-2, 1.8997E-2, 2.0784E-2, 2.2797E-2, 2.5066E-2, 2.7622E-2, 2.9020E-2, 3.0503E-2},
1041 {0.04, 3.8474E-2, 4.0056E-2, 4.1718E-2, 4.5 1073 {0.04, 3.8474E-2, 4.0056E-2, 4.1718E-2, 4.5300E-2, 4.9254E-2, 5.3619E-2, 5.8436E-2, 6.1028E-2, 6.3752E-2},
1042 {0.05, 6.7371E-2, 6.9928E-2, 7.2596E-2, 7.8 1074 {0.05, 6.7371E-2, 6.9928E-2, 7.2596E-2, 7.8282E-2, 8.4470E-2, 9.1206E-2, 9.8538E-2, 1.0244E-1, 1.0652E-1},
1043 {0.06, 1.0418E-1, 1.0778E-1, 1.1152E-1, 1.1 1075 {0.06, 1.0418E-1, 1.0778E-1, 1.1152E-1, 1.1943E-1, 1.2796E-1, 1.3715E-1, 1.4706E-1, 1.5231E-1, 1.5776E-1},
1044 {0.08, 1.9647E-1, 2.0217E-1, 2.0805E-1, 2.2 1076 {0.08, 1.9647E-1, 2.0217E-1, 2.0805E-1, 2.2038E-1, 2.3351E-1, 2.4751E-1, 2.6244E-1, 2.7027E-1, 2.7837E-1},
1045 {0.1, 3.0594E-1, 3.1361E-1, 3.2150E-1, 3.37 1077 {0.1, 3.0594E-1, 3.1361E-1, 3.2150E-1, 3.3796E-1, 3.5537E-1, 3.7381E-1, 3.9336E-1, 4.0357E-1, 4.1410E-1},
1046 {0.15, 6.1411E-1, 6.2597E-1, 6.3811E-1, 6.6 1078 {0.15, 6.1411E-1, 6.2597E-1, 6.3811E-1, 6.6330E-1, 6.8974E-1, 7.1753E-1, 7.4678E-1, 7.6199E-1, 7.7761E-1},
1047 {0.2, 9.3100E-1, 9.5614E-1, 9.7162E-1, 1.00 1079 {0.2, 9.3100E-1, 9.5614E-1, 9.7162E-1, 1.0037, 1.0372, 1.0723, 1.1092, 1.1284, 1.1480},
1048 {0.3, 1.5172, 1.5372, 1.5576, 1.5998, 1.643 1080 {0.3, 1.5172, 1.5372, 1.5576, 1.5998, 1.6438, 1.6899, 1.7382, 1.7632, 1.7889},
1049 {0.4, 2.0173, 2.0408, 2.0647, 2.1142, 2.165 1081 {0.4, 2.0173, 2.0408, 2.0647, 2.1142, 2.1659, 2.2199, 2.2765, 2.3059, 2.3360},
1050 {0.5, 2.3932, 2.4193, 2.4460, 2.5011, 2.558 1082 {0.5, 2.3932, 2.4193, 2.4460, 2.5011, 2.5587, 2.6190, 2.6821, 2.7148, 2.7484},
1051 {0.6, 2.7091, 2.7374, 2.7663, 2.8260, 2.888 1083 {0.6, 2.7091, 2.7374, 2.7663, 2.8260, 2.8884, 2.9538, 3.0222, 3.0577, 3.0941},
1052 {0.7, 2.9742, 3.0044, 3.0352, 3.0988, 3.165 1084 {0.7, 2.9742, 3.0044, 3.0352, 3.0988, 3.1652, 3.2349, 3.3079, 3.3457, 3.3845},
1053 {0.8, 3.2222, 3.2539, 3.2863, 3.3532, 3.423 1085 {0.8, 3.2222, 3.2539, 3.2863, 3.3532, 3.4232, 3.4965, 3.5734, 3.6133, 3.6542},
1054 {1.0, 3.6690, 3.7033, 3.7384, 3.8108, 3.886 1086 {1.0, 3.6690, 3.7033, 3.7384, 3.8108, 3.8866, 3.9661, 4.0495, 4.0928, 4.1371},
1055 {1.2, 3.9819, 4.0183, 4.0555, 4.1324, 4.213 1087 {1.2, 3.9819, 4.0183, 4.0555, 4.1324, 4.2130, 4.2974, 4.3861, 4.4321, 4.4794},
1056 {1.4, 4.2745, 4.3127, 4.3517, 4.4324, 4.517 1088 {1.4, 4.2745, 4.3127, 4.3517, 4.4324, 4.5170, 4.6056, 4.6988, 4.7471, 4.7968},
1057 {1.5, 4.4047, 4.4436, 4.4834, 4.5658, 4.652 1089 {1.5, 4.4047, 4.4436, 4.4834, 4.5658, 4.6521, 4.7426, 4.8378, 4.8872, 4.9379},
1058 {1.7, 4.6383, 4.6787, 4.7200, 4.8054, 4.894 1090 {1.7, 4.6383, 4.6787, 4.7200, 4.8054, 4.8949, 4.9888, 5.0876, 5.1388, 5.1915},
1059 {2.0, 4.9369, 4.9791, 5.0223, 5.1116, 5.205 1091 {2.0, 4.9369, 4.9791, 5.0223, 5.1116, 5.2053, 5.3036, 5.4070, 5.4607, 5.5159},
1060 {3.0, 5.6514, 5.6981, 5.7459, 5.8450, 5.948 1092 {3.0, 5.6514, 5.6981, 5.7459, 5.8450, 5.9489, 6.0581, 6.1730, 6.2328, 6.2943},
1061 {5.0, 6.4665, 6.5189, 6.5724, 6.6835, 6.800 1093 {5.0, 6.4665, 6.5189, 6.5724, 6.6835, 6.8003, 6.9231, 7.0525, 7.1199, 7.1892},
1062 {7.0, 6.8634, 6.9194, 6.9767, 7.0957, 7.220 1094 {7.0, 6.8634, 6.9194, 6.9767, 7.0957, 7.2208, 7.3526, 7.4915, 7.5639, 7.6384}
1063 }; 1095 };
1064 1096
1065 const G4double bls2[28][10] = { << 1097 static const G4double bls2[28][10] = {
1066 {0.005, 5.4561E-4, 6.0905E-4, 6.7863E-4, 7. 1098 {0.005, 5.4561E-4, 6.0905E-4, 6.7863E-4, 7.5494E-4, 7.9587E-4, 8.3883E-4, 9.3160E-4, 1.0352E-3, 1.1529E-3},
1067 {0.007, 1.2068E-3, 1.3170E-3, 1.4377E-3, 1. 1099 {0.007, 1.2068E-3, 1.3170E-3, 1.4377E-3, 1.5719E-3, 1.6451E-3, 1.7231E-3, 1.8969E-3, 2.1009E-3, 2.3459E-3},
1068 {0.01, 2.6832E-3, 2.9017E-3, 3.1534E-3, 3.4 1100 {0.01, 2.6832E-3, 2.9017E-3, 3.1534E-3, 3.4479E-3, 3.6149E-3, 3.7976E-3, 4.2187E-3, 4.7320E-3, 5.3636E-3},
1069 {0.015, 6.2775E-3, 6.9077E-3, 7.6525E-3, 8. 1101 {0.015, 6.2775E-3, 6.9077E-3, 7.6525E-3, 8.5370E-2, 9.0407E-3, 9.5910E-3, 1.0850E-2, 1.2358E-2, 1.4165E-2},
1070 {0.02, 1.2561E-2, 1.3943E-2, 1.5553E-2, 1.7 1102 {0.02, 1.2561E-2, 1.3943E-2, 1.5553E-2, 1.7327E-2, 1.8478E-2, 1.9612E-2, 2.2160E-2, 2.5130E-2, 2.8594E-2},
1071 {0.03, 3.3750E-2, 3.7470E-2, 4.1528E-2, 4.6 1103 {0.03, 3.3750E-2, 3.7470E-2, 4.1528E-2, 4.6170E-2, 4.8708E-2, 5.1401E-2, 5.7297E-2, 6.3943E-2, 7.1441E-2},
1072 {0.04, 6.9619E-2, 7.6098E-2, 8.3249E-2, 9.1 1104 {0.04, 6.9619E-2, 7.6098E-2, 8.3249E-2, 9.1150E-2, 9.5406E-2, 9.9881E-2, 1.0954E-1, 1.2023E-1, 1.3208E-1},
1073 {0.05, 1.1522E-1, 1.2470E-1, 1.3504E-1, 1.4 1105 {0.05, 1.1522E-1, 1.2470E-1, 1.3504E-1, 1.4632E-1, 1.5234E-1, 1.5864E-1, 1.7211E-1, 1.8686E-1, 2.0304E-1},
1074 {0.06, 1.6931E-1, 1.8179E-1, 1.9530E-1, 2.0 1106 {0.06, 1.6931E-1, 1.8179E-1, 1.9530E-1, 2.0991E-1, 2.1767E-1, 2.2576E-1, 2.4295E-1, 2.6165E-1, 2.8201E-1},
1075 {0.08, 2.9540E-1, 3.1361E-1, 3.3312E-1, 3.5 1107 {0.08, 2.9540E-1, 3.1361E-1, 3.3312E-1, 3.5403E-1, 3.6506E-1, 3.7650E-1, 4.0067E-1, 4.2673E-1, 4.5488E-1},
1076 {0.1, 4.3613E-1, 4.5956E-1, 4.852E-1, 5.1 1108 {0.1, 4.3613E-1, 4.5956E-1, 4.852E-1, 5.1115E-1, 5.2514E-1, 5.3961E-1, 5.7008E-1, 6.0277E-1, 6.3793E-1},
1077 {0.15, 8.1014E-1, 8.4453E-1, 8.8093E-1, 9.1 1109 {0.15, 8.1014E-1, 8.4453E-1, 8.8093E-1, 9.1954E-1, 9.3973E-1, 9.6056E-1, 1.0043, 1.0509, 1.1008},
1078 {0.2, 1.1888, 1.2319, 1.2774, 1.3255, 1.350 1110 {0.2, 1.1888, 1.2319, 1.2774, 1.3255, 1.3506, 1.3765, 1.4308, 1.4886, 1.5504},
1079 {0.3, 1.8422, 1.8983, 1.9575, 2.0201, 2.052 1111 {0.3, 1.8422, 1.8983, 1.9575, 2.0201, 2.0528, 2.0864, 2.1569, 2.2319, 2.3120},
1080 {0.4, 2.3984, 2.4642, 2.5336, 2.6070, 2.645 1112 {0.4, 2.3984, 2.4642, 2.5336, 2.6070, 2.6452, 2.6847, 2.7674, 2.8554, 2.9494},
1081 {0.5, 2.8181, 2.8915, 2.9690, 3.0509, 3.093 1113 {0.5, 2.8181, 2.8915, 2.9690, 3.0509, 3.0937, 3.1378, 3.2301, 3.3285, 3.4337},
1082 {0.6, 3.1698, 3.2494, 3.3336, 3.4226, 3.469 1114 {0.6, 3.1698, 3.2494, 3.3336, 3.4226, 3.4691, 3.5171, 3.6175, 3.7246, 3.8391},
1083 {0.7, 3.4652, 3.5502, 3.6400, 3.7351, 3.784 1115 {0.7, 3.4652, 3.5502, 3.6400, 3.7351, 3.7848, 3.8360, 3.9433, 4.0578, 4.1804},
1084 {0.8, 3.7392, 3.8289, 3.9236, 4.0239, 4.076 1116 {0.8, 3.7392, 3.8289, 3.9236, 4.0239, 4.0764, 4.1304, 4.2438, 4.3648, 4.4944},
1085 {1.0, 4.2295, 4.3269, 4.4299, 4.5391, 4.596 1117 {1.0, 4.2295, 4.3269, 4.4299, 4.5391, 4.5962, 4.6551, 4.7786, 4.9106, 5.0520},
1086 {1.2, 4.5777, 4.6814, 4.7912, 4.9076, 4.968 1118 {1.2, 4.5777, 4.6814, 4.7912, 4.9076, 4.9685, 5.0314, 5.1633, 5.3043, 5.4555},
1087 {1.4, 4.9001, 5.0092, 5.1247, 5.2473, 5.311 1119 {1.4, 4.9001, 5.0092, 5.1247, 5.2473, 5.3114, 5.3776, 5.5166, 5.6653, 5.8249},
1088 {1.5, 5.0434, 5.1550, 5.2731, 5.3984, 5.464 1120 {1.5, 5.0434, 5.1550, 5.2731, 5.3984, 5.4640, 5.5317, 5.6739, 5.8260, 5.9893},
1089 {1.7, 5.3011, 5.4170, 5.5398, 5.6701, 5.737 1121 {1.7, 5.3011, 5.4170, 5.5398, 5.6701, 5.7373, 5.8088, 5.9568, 6.1152, 6.2853},
1090 {2.0, 5.6308, 5.7523, 5.8811, 6.0174, 6.089 1122 {2.0, 5.6308, 5.7523, 5.8811, 6.0174, 6.0896, 6.1636, 6.3192, 6.4857, 6.6647},
1091 {3.0, 6.4224, 6.5580, 6.7019, 6.8549, 6.935 1123 {3.0, 6.4224, 6.5580, 6.7019, 6.8549, 6.9351, 7.0180, 7.1925, 7.3795, 7.5808},
1092 {5.0, 7.3338, 7.4872, 7.6500, 7.8235, 7.914 1124 {5.0, 7.3338, 7.4872, 7.6500, 7.8235, 7.9146, 8.0087, 8.2071, 8.4200, 8.6496},
1093 {7.0, 7.7938, 7.9588, 8.1342, 8.3211, 8.419 1125 {7.0, 7.7938, 7.9588, 8.1342, 8.3211, 8.4193, 8.5209, 8.7350, 8.9651, 9.2133}
1094 }; 1126 };
1095 1127
1096 const G4double bls3[28][9] = { << 1128 static const G4double bls3[28][9] = {
1097 {0.005, 1.2895E-3, 1.3670E-3, 1.4524E-3, 1. 1129 {0.005, 1.2895E-3, 1.3670E-3, 1.4524E-3, 1.6524E-3, 1.9078E-3, 2.2414E-3, 2.6889E-3, 3.3006E-3},
1098 {0.007, 2.6467E-3, 2.8242E-3, 3.0238E-3, 3. 1130 {0.007, 2.6467E-3, 2.8242E-3, 3.0238E-3, 3.5045E-3, 4.1260E-3, 4.9376E-3, 6.0050E-3, 7.4152E-3},
1099 {0.01, 6.1472E-3, 6.6086E-3, 7.1246E-3, 8.3 1131 {0.01, 6.1472E-3, 6.6086E-3, 7.1246E-3, 8.3491E-3, 9.8871E-3, 1.1822E-2, 1.4261E-2, 1.7335E-2},
1100 {0.015, 1.63316E-2, 1.7572E-2, 1.8932E-2, 2 1132 {0.015, 1.63316E-2, 1.7572E-2, 1.8932E-2, 2.2053E-2, 2.5803E-2, 3.0308E-2, 3.5728E-2, 4.2258E-2},
1101 {0.02, 3.2634E-2, 3.4900E-2, 3.7348E-2, 4.2 1133 {0.02, 3.2634E-2, 3.4900E-2, 3.7348E-2, 4.2850E-2, 4.9278E-2, 5.6798E-2, 6.5610E-2, 7.5963E-2},
1102 {0.03, 7.9907E-2, 8.4544E-2, 8.9486E-2, 1.0 1134 {0.03, 7.9907E-2, 8.4544E-2, 8.9486E-2, 1.0032E-1, 1.1260E-1, 1.2656E-1, 1.4248E-1, 1.6071E-1},
1103 {0.04, 1.4523E-1, 1.5237E-1, 1.5985E-1, 1.7 1135 {0.04, 1.4523E-1, 1.5237E-1, 1.5985E-1, 1.7614E-1, 1.9434E-1, 2.1473E-1, 2.3766E-1, 2.6357E-1},
1104 {0.05, 2.2082E-1, 2.3036E-1, 2.4038E-1, 2.6 1136 {0.05, 2.2082E-1, 2.3036E-1, 2.4038E-1, 2.6199E-1, 2.8590E-1, 3.1248E-1, 3.4212E-1, 3.7536E-1},
1105 {0.06, 3.0423E-1, 3.1611E-1, 3.2854E-1, 3.5 1137 {0.06, 3.0423E-1, 3.1611E-1, 3.2854E-1, 3.5522E-1, 3.8459E-1, 4.1704E-1, 4.5306E-1, 4.9326E-1},
1106 {0.08, 4.8536E-1, 5.0156E-1, 5.1846E-1, 5.5 1138 {0.08, 4.8536E-1, 5.0156E-1, 5.1846E-1, 5.5453E-1, 5.9397E-1, 6.3728E-1, 6.8507E-1, 7.3810E-1},
1107 {0.1, 6.7586E-1, 6.9596E-1, 7.1688E-1, 7.61 1139 {0.1, 6.7586E-1, 6.9596E-1, 7.1688E-1, 7.6141E-1, 8.0992E-1, 8.6301E-1, 9.2142E-1, 9.8604E-1},
1108 {0.15, 1.1544, 1.1828, 1.2122, 1.2746, 1.34 1140 {0.15, 1.1544, 1.1828, 1.2122, 1.2746, 1.3424, 1.4163, 1.4974, 1.5868},
1109 {0.2, 1.6167, 1.6517, 1.6880, 1.7650, 1.848 1141 {0.2, 1.6167, 1.6517, 1.6880, 1.7650, 1.8484, 1.9394, 2.0390, 2.1489},
1110 {0.3, 2.3979, 2.4432, 2.4902, 2.5899, 2.698 1142 {0.3, 2.3979, 2.4432, 2.4902, 2.5899, 2.6980, 2.8159, 2.9451, 3.0876},
1111 {0.4, 3.0502, 3.1034, 3.1586, 3.2758, 3.403 1143 {0.4, 3.0502, 3.1034, 3.1586, 3.2758, 3.4030, 3.5416, 3.6938, 3.8620},
1112 {0.5, 3.5466, 3.6062, 3.6681, 3.7994, 3.942 1144 {0.5, 3.5466, 3.6062, 3.6681, 3.7994, 3.9421, 4.0978, 4.2688, 4.4580},
1113 {0.6, 3.9620, 4.0270, 4.0945, 4.2378, 4.393 1145 {0.6, 3.9620, 4.0270, 4.0945, 4.2378, 4.3935, 4.5636, 4.7506, 4.9576},
1114 {0.7, 4.3020, 4.3715, 4.4438, 4.5974, 4.764 1146 {0.7, 4.3020, 4.3715, 4.4438, 4.5974, 4.7644, 4.9470, 5.1478, 5.3703},
1115 {0.8, 4.6336, 4.7072, 4.7837, 4.9463, 5.123 1147 {0.8, 4.6336, 4.7072, 4.7837, 4.9463, 5.1233, 5.3169, 5.5300, 5.7661},
1116 {1.0, 5.2041, 5.2845, 5.3682, 5.5462, 5.740 1148 {1.0, 5.2041, 5.2845, 5.3682, 5.5462, 5.7400, 5.9523, 6.1863, 6.4458},
1117 {1.2, 5.6182, 5.7044, 5.7940, 5.9848, 6.192 1149 {1.2, 5.6182, 5.7044, 5.7940, 5.9848, 6.1927, 6.4206, 6.6719, 6.9510},
1118 {1.4, 5.9967, 6.0876, 6.1823, 6.3839, 6.603 1150 {1.4, 5.9967, 6.0876, 6.1823, 6.3839, 6.6038, 6.8451, 7.1113, 7.4071},
1119 {1.5, 6.1652, 6.2583, 6.3553, 6.5618, 6.787 1151 {1.5, 6.1652, 6.2583, 6.3553, 6.5618, 6.7871, 7.0344, 7.3073, 7.6107},
1120 {1.7, 6.4686, 6.5657, 6.6668, 6.8823, 7.117 1152 {1.7, 6.4686, 6.5657, 6.6668, 6.8823, 7.1175, 7.3757, 7.6609, 7.9782},
1121 {2.0, 6.8577, 6.9600, 7.0666, 7.2937, 7.541 1153 {2.0, 6.8577, 6.9600, 7.0666, 7.2937, 7.5418, 7.8143, 8.1156, 8.4510},
1122 {3.0, 7.7981, 7.9134, 8.0336, 8.2901, 8.570 1154 {3.0, 7.7981, 7.9134, 8.0336, 8.2901, 8.5708, 8.8796, 9.2214, 9.6027},
1123 {5.0, 8.8978, 9.0297, 9.1673, 9.4612, 9.783 1155 {5.0, 8.8978, 9.0297, 9.1673, 9.4612, 9.7834, 10.1384, 10.5323, 10.9722},
1124 {7.0, 9.4819, 9.6248, 9.7739, 10.0926, 10.4 1156 {7.0, 9.4819, 9.6248, 9.7739, 10.0926, 10.4423, 10.8282, 11.2565, 11.7356}
1125 }; 1157 };
1126 1158
1127 const G4double bll1[28][10] = { << 1159 static const G4double bll1[28][10] = {
1128 {0.005, 3.6324E-5, 4.0609E-5, 4.5430E-5, 5. 1160 {0.005, 3.6324E-5, 4.0609E-5, 4.5430E-5, 5.6969E-5, 7.1625E-5, 9.0279E-5, 1.1407E-4, 1.2834E-4, 1.4447E-4},
1129 {0.007, 1.8110E-4, 2.0001E-4, 2.2099E-4, 2. 1161 {0.007, 1.8110E-4, 2.0001E-4, 2.2099E-4, 2.7006E-4, 3.3049E-4, 4.0498E-4, 4.9688E-4, 5.5061E-4, 6.1032E-4},
1130 {0.01, 8.6524E-4, 9.4223E-4, 1.0262E-3, 1.2 1162 {0.01, 8.6524E-4, 9.4223E-4, 1.0262E-3, 1.2178E-3, 1.4459E-3, 1.7174E-3, 2.0405E-3, 2.2245E-3, 2.4252E-3},
1131 {0.015, 4.2293E-3, 4.5314E-3, 4.8551E-3, 5. 1163 {0.015, 4.2293E-3, 4.5314E-3, 4.8551E-3, 5.5731E-3, 6.3968E-3, 7.3414E-3, 8.4242E-3, 9.0236E-3, 9.6652E-3},
1132 {0.02, 1.1485E-2, 1.2172E-2, 1.2900E-2, 1.4 1164 {0.02, 1.1485E-2, 1.2172E-2, 1.2900E-2, 1.4486E-2, 1.6264E-2, 1.8256E-2, 2.0487E-2, 2.1702E-2, 2.2989E-2},
1133 {0.03, 3.9471E-2, 4.1270E-2, 4.3149E-2, 4.7 1165 {0.03, 3.9471E-2, 4.1270E-2, 4.3149E-2, 4.7163E-2, 5.1543E-2, 5.6423E-2, 6.1540E-2, 6.4326E-2, 6.7237E-2},
1134 {0.04, 8.4454E-2, 8.7599E-2, 9.0860E-2, 9.7 1166 {0.04, 8.4454E-2, 8.7599E-2, 9.0860E-2, 9.7747E-2, 1.0516E-1, 1.1313E-1, 1.2171E-1, 1.2625E-1, 1.3096E-1},
1135 {0.05, 1.4339E-1, 1.4795E-1, 1.5266E-1, 1.6 1167 {0.05, 1.4339E-1, 1.4795E-1, 1.5266E-1, 1.6253E-1, 1.7306E-1, 1.8430E-1, 1.9630E-1, 2.0261E-1, 2.0924E-1},
1136 {0.06, 2.1304E-1, 2.1899E-1, 2.2512E-1, 2.3 1168 {0.06, 2.1304E-1, 2.1899E-1, 2.2512E-1, 2.3794E-1, 2.5153E-1, 2.6596E-1, 2.8130E-1, 2.8934E-1, 2.9763E-1},
1137 {0.08, 3.7382E-1, 3.8241E-1, 3.9122E-1, 4.0 1169 {0.08, 3.7382E-1, 3.8241E-1, 3.9122E-1, 4.0955E-1, 4.2888E-1, 4.4928E-1, 4.7086E-1, 4.8212E-1, 4.9371E-1},
1138 {0.1, 5.5056E-1, 5.6151E-1, 5.7273E-1, 5.96 1170 {0.1, 5.5056E-1, 5.6151E-1, 5.7273E-1, 5.9601E-1, 6.2049E-1, 6.4627E-1, 6.7346E-1, 6.8762E-1, 7.0218E-1},
1139 {0.15, 1.0066, 1.0224, 1.0386, 1.0721, 1.10 1171 {0.15, 1.0066, 1.0224, 1.0386, 1.0721, 1.1073, 1.1443, 1.1832, 1.2035, 1.2243},
1140 {0.2, 1.4376, 1.4572, 1.4773, 1.5188, 1.562 1172 {0.2, 1.4376, 1.4572, 1.4773, 1.5188, 1.5624, 1.6083, 1.6566, 1.6817, 1.7076},
1141 {0.3, 2.1712, 2.1964, 2.2223, 2.2758, 2.332 1173 {0.3, 2.1712, 2.1964, 2.2223, 2.2758, 2.3322, 2.3915, 2.4542, 2.4869, 2.5205},
1142 {0.4, 2.7500, 2.7793, 2.8094, 2.8719, 2.937 1174 {0.4, 2.7500, 2.7793, 2.8094, 2.8719, 2.9377, 3.0072, 3.0807, 3.1192, 3.1587},
1143 {0.5, 3.2033, 3.2359, 3.2693, 3.3389, 3.412 1175 {0.5, 3.2033, 3.2359, 3.2693, 3.3389, 3.4122, 3.4898, 3.5721, 3.6151, 3.6595},
1144 {0.6, 3.6038, 3.6391, 3.6753, 3.7506, 3.830 1176 {0.6, 3.6038, 3.6391, 3.6753, 3.7506, 3.8303, 3.9146, 4.0042, 4.0511, 4.0995},
1145 {0.7, 3.9106, 3.9482, 3.9867, 4.0670, 4.152 1177 {0.7, 3.9106, 3.9482, 3.9867, 4.0670, 4.1520, 4.2421, 4.3380, 4.3882, 4.4401},
1146 {0.8, 4.1790, 4.2185, 4.2590, 4.3437, 4.433 1178 {0.8, 4.1790, 4.2185, 4.2590, 4.3437, 4.4333, 4.5285, 4.6298, 4.6830, 4.7380},
1147 {1.0, 4.6344, 4.6772, 4.7212, 4.8131, 4.910 1179 {1.0, 4.6344, 4.6772, 4.7212, 4.8131, 4.9106, 5.0144, 5.1250, 5.1831, 5.2432},
1148 {1.2, 4.9787, 5.0242, 5.0711, 5.1689, 5.272 1180 {1.2, 4.9787, 5.0242, 5.0711, 5.1689, 5.2729, 5.3837, 5.5050, 5.5642, 5.6287},
1149 {1.4, 5.2688, 5.3166, 5.3658, 5.4687, 5.578 1181 {1.4, 5.2688, 5.3166, 5.3658, 5.4687, 5.5782, 5.6950, 5.8198, 5.8855, 5.9554},
1150 {1.5, 5.3966, 5.4454, 5.4957, 5.6009, 5.712 1182 {1.5, 5.3966, 5.4454, 5.4957, 5.6009, 5.7128, 5.8323, 5.9601, 6.0274, 6.0972},
1151 {1.7, 5.6255, 5.6762, 5.7284, 5.8377, 5.954 1183 {1.7, 5.6255, 5.6762, 5.7284, 5.8377, 5.9541, 6.0785, 6.2116, 6.2818, 6.3546},
1152 {2.0, 5.9170, 5.9701, 6.0248, 6.1395, 6.261 1184 {2.0, 5.9170, 5.9701, 6.0248, 6.1395, 6.2618, 6.3925, 6.5327, 6.6066, 6.6833},
1153 {3.0, 6.6210, 6.6801, 6.7411, 6.8692, 7.006 1185 {3.0, 6.6210, 6.6801, 6.7411, 6.8692, 7.0062, 7.1529, 7.3107, 7.3941, 7.4807},
1154 {5.0, 7.4620, 7.5288, 7.5977, 7.7428, 7.898 1186 {5.0, 7.4620, 7.5288, 7.5977, 7.7428, 7.8982, 8.0653, 8.2454, 8.3409, 8.4402},
1155 {7.0, 7.7362, 7.8079, 7.8821, 8.0383, 8.206 1187 {7.0, 7.7362, 7.8079, 7.8821, 8.0383, 8.2061, 8.3866, 8.5816, 8.6850, 8.7927}
1156 }; 1188 };
1157 1189
1158 const G4double bll2[28][10] = { << 1190 static const G4double bll2[28][10] = {
1159 {0.005, 1.8339E-4, 2.3330E-4, 2.9738E-4, 3. 1191 {0.005, 1.8339E-4, 2.3330E-4, 2.9738E-4, 3.7977E-4, 4.2945E-4, 4.8582E-4, 6.2244E-4, 7.9858E-4, 1.0258E-3},
1160 {0.007, 7.5042E-4, 9.2355E-4, 1.1375E-3, 1. 1192 {0.007, 7.5042E-4, 9.2355E-4, 1.1375E-3, 1.4021E-3, 1.5570E-3, 1.7292E-3, 2.1335E-3, 2.6335E-3, 3.2515E-3},
1161 {0.01, 2.8829E-3, 3.4275E-3, 4.0758E-3, 4.8 1193 {0.01, 2.8829E-3, 3.4275E-3, 4.0758E-3, 4.8457E-3, 5.2839E-3, 5.7617E-3, 6.8504E-3, 8.1442E-3, 9.6816E-3},
1162 {0.015, 1.1087E-2, 1.2716E-2, 1.4581E-2, 1. 1194 {0.015, 1.1087E-2, 1.2716E-2, 1.4581E-2, 1.6717E-2, 1.7898E-2, 1.9163E-2, 2.1964E-2, 2.5173E-2, 2.8851E-2},
1163 {0.02, 2.5786E-2, 2.8922E-2, 3.2435E-2, 3.6 1195 {0.02, 2.5786E-2, 2.8922E-2, 3.2435E-2, 3.6371E-2, 3.8514E-2, 4.0784E-2, 4.5733E-2, 5.1288E-2, 5.7531E-2},
1164 {0.03, 7.3461E-2, 8.0264E-2, 8.7705E-2, 9.5 1196 {0.03, 7.3461E-2, 8.0264E-2, 8.7705E-2, 9.5852E-2, 1.0021E-1, 1.0478E-1, 1.1458E-1, 1.2535E-1, 1.3721E-1},
1165 {0.04, 1.4094E-1, 1.5172E-1, 1.6336E-1, 1.7 1197 {0.04, 1.4094E-1, 1.5172E-1, 1.6336E-1, 1.7596E-1, 1.8265E-1, 1.8962E-1, 2.0445E-1, 2.2058E-1, 2.3818E-1},
1166 {0.05, 2.2289E-1, 2.3762E-1, 2.5344E-1, 2.7 1198 {0.05, 2.2289E-1, 2.3762E-1, 2.5344E-1, 2.7045E-1, 2.7944E-1, 2.8877E-1, 3.0855E-1, 3.2995E-1, 3.5318E-1},
1167 {0.06, 3.1503E-1, 3.3361E-1, 3.5346E-1, 3.7 1199 {0.06, 3.1503E-1, 3.3361E-1, 3.5346E-1, 3.7473E-1, 3.8594E-1, 3.9756E-1, 4.2212E-1, 4.4861E-1, 4.7727E-1},
1168 {0.08, 5.1793E-1, 5.4368E-1, 5.7109E-1, 6.0 1200 {0.08, 5.1793E-1, 5.4368E-1, 5.7109E-1, 6.0032E-1, 6.1569E-1, 6.3159E-1, 6.6512E-1, 7.0119E-1, 7.4012E-1},
1169 {0.1, 7.3258E-1, 7.6481E-1, 7.9907E-1, 8.35 1201 {0.1, 7.3258E-1, 7.6481E-1, 7.9907E-1, 8.3556E-1, 8.5472E-1, 8.7454E-1, 9.1630E-1, 9.6119E-1, 1.0096},
1170 {0.15, 1.2677, 1.3137, 1.3626, 1.4147, 1.44 1202 {0.15, 1.2677, 1.3137, 1.3626, 1.4147, 1.4421, 1.4703, 1.5300, 1.5942, 1.6636},
1171 {0.2, 1.7615, 1.8188, 1.8797, 1.9446, 1.978 1203 {0.2, 1.7615, 1.8188, 1.8797, 1.9446, 1.9788, 2.0142, 2.0889, 2.1694, 2.2567},
1172 {0.3, 2.5909, 2.6658, 2.7457, 2.8312, 2.876 1204 {0.3, 2.5909, 2.6658, 2.7457, 2.8312, 2.8762, 2.9231, 3.0222, 3.1295, 3.2463},
1173 {0.4, 3.2417, 3.3302, 3.4249, 3.5266, 3.580 1205 {0.4, 3.2417, 3.3302, 3.4249, 3.5266, 3.5803, 3.6361, 3.7546, 3.8835, 4.0242},
1174 {0.5, 3.7527, 3.8523, 3.9591, 4.0741, 4.135 1206 {0.5, 3.7527, 3.8523, 3.9591, 4.0741, 4.1350, 4.1983, 4.3330, 4.4799, 4.6408},
1175 {0.6, 4.2013, 4.3103, 4.4274, 4.5537, 4.620 1207 {0.6, 4.2013, 4.3103, 4.4274, 4.5537, 4.6206, 4.6904, 4.8390, 5.0013, 5.1796},
1176 {0.7, 4.5493, 4.6664, 4.7925, 4.9286, 5.000 1208 {0.7, 4.5493, 4.6664, 4.7925, 4.9286, 5.0009, 5.0762, 5.2370, 5.4129, 5.6066},
1177 {0.8, 4.8537, 4.9780, 5.1119, 5.2568, 5.333 1209 {0.8, 4.8537, 4.9780, 5.1119, 5.2568, 5.3338, 5.4141, 5.5857, 5.7738, 5.9811},
1178 {1.0, 5.3701, 5.5066, 5.6540, 5.8138, 5.898 1210 {1.0, 5.3701, 5.5066, 5.6540, 5.8138, 5.8989, 5.9878, 6.1780, 6.3870, 6.6179},
1179 {1.2, 5.7648, 5.9114, 6.0701, 6.2424, 6.334 1211 {1.2, 5.7648, 5.9114, 6.0701, 6.2424, 6.3343, 6.4303, 6.6361, 6.8626, 7.1137},
1180 {1.4, 6.0976, 6.2530, 6.4214, 6.6044, 6.702 1212 {1.4, 6.0976, 6.2530, 6.4214, 6.6044, 6.7021, 6.8043, 7.0237, 7.2655, 7.5338},
1181 {1.5, 6.2447, 6.4041, 6.5768, 6.7647, 6.865 1213 {1.5, 6.2447, 6.4041, 6.5768, 6.7647, 6.8650, 6.9700, 7.1954, 7.4442, 7.7203},
1182 {1.7, 6.5087, 6.6752, 6.8558, 7.0526, 7.157 1214 {1.7, 6.5087, 6.6752, 6.8558, 7.0526, 7.1578, 7.2679, 7.5045, 7.7660, 8.0565},
1183 {2.0, 6.8458, 7.0218, 7.2129, 7.4213, 7.532 1215 {2.0, 6.8458, 7.0218, 7.2129, 7.4213, 7.5328, 7.6496, 7.9010, 8.1791, 8.4886},
1184 {3.0, 7.6647, 7.8644, 8.0819, 8.3189, 8.447 1216 {3.0, 7.6647, 7.8644, 8.0819, 8.3189, 8.4475, 8.5814, 8.8702, 9.1908, 9.5488},
1185 {5.0, 8.6515, 8.8816, 9.1330, 9.4090, 9.557 1217 {5.0, 8.6515, 8.8816, 9.1330, 9.4090, 9.5572, 9.7132, 10.0504, 10.4259, 10.8466},
1186 {7.0, 9.0221, 9.2724, 9.5464, 9.8477, 10.00 1218 {7.0, 9.0221, 9.2724, 9.5464, 9.8477, 10.0099, 10.1805, 10.5499, 10.9622, 11.4250}
1187 }; 1219 };
1188 1220
1189 const G4double bll3[28][9] = { << 1221 static const G4double bll3[28][9] = {
1190 {0.005, 1.3190E-3, 1.4961E-3, 1.6974E-3, 2. 1222 {0.005, 1.3190E-3, 1.4961E-3, 1.6974E-3, 2.1858E-3, 2.8163E-3, 3.6302E-3, 4.6814E-3, 6.0395E-3},
1191 {0.007, 4.0158E-3, 4.4623E-3, 4.9592E-3, 6. 1223 {0.007, 4.0158E-3, 4.4623E-3, 4.9592E-3, 6.1257E-3, 7.5675E-3, 9.3502E-3, 1.1556E-2, 1.4290E-2},
1192 {0.01, 1.1509E-2, 1.2548E-2, 1.3681E-2, 1.6 1224 {0.01, 1.1509E-2, 1.2548E-2, 1.3681E-2, 1.6263E-2, 1.9336E-2, 2.2999E-2, 2.7370E-2, 3.2603E-2},
1193 {0.015, 3.3070E-2, 3.5408E-2, 3.7914E-2, 4. 1225 {0.015, 3.3070E-2, 3.5408E-2, 3.7914E-2, 4.3483E-2, 4.9898E-2, 5.7304E-2, 6.5884E-2, 7.5861E-2},
1194 {0.02, 6.4555E-2, 6.8394E-2, 7.2472E-2, 8.1 1226 {0.02, 6.4555E-2, 6.8394E-2, 7.2472E-2, 8.1413E-2, 9.1539E-2, 1.0304E-1, 1.1617E-1, 1.3121E-1},
1195 {0.03, 1.5030E-1, 1.5101E-1, 1.5844E-1, 1.7 1227 {0.03, 1.5030E-1, 1.5101E-1, 1.5844E-1, 1.7451E-1, 1.9244E-1, 2.1244E-1, 2.3496E-1, 2.6044E-1},
1196 {0.04, 2.5743E-1, 2.6774E-1, 2.7855E-1, 3.0 1228 {0.04, 2.5743E-1, 2.6774E-1, 2.7855E-1, 3.0180E-1, 3.2751E-1, 3.5608E-1, 3.8803E-1, 4.2401E-1},
1197 {0.05, 3.7846E-1, 3.9195E-1, 4.0607E-1, 4.3 1229 {0.05, 3.7846E-1, 3.9195E-1, 4.0607E-1, 4.3635E-1, 4.6973E-1, 5.0672E-1, 5.4798E-1, 5.9436E-1},
1198 {0.06, 5.0839E-1, 5.2497E-1, 5.4230E-1, 5.7 1230 {0.06, 5.0839E-1, 5.2497E-1, 5.4230E-1, 5.7943E-1, 6.2028E-1, 6.6549E-1, 7.1589E-1, 7.7252E-1},
1199 {0.08, 7.8230E-1, 8.0474E-1, 8.2818E-1, 8.7 1231 {0.08, 7.8230E-1, 8.0474E-1, 8.2818E-1, 8.7836E-1, 9.3355E-1, 9.9462E-1, 1.0627, 1.1394},
1200 {0.1, 1.0621, 1.0900, 1.1192, 1.1816, 1.250 1232 {0.1, 1.0621, 1.0900, 1.1192, 1.1816, 1.2503, 1.3265, 1.4116, 1.5076},
1201 {0.15, 1.7389, 1.7790, 1.8210, 1.9112, 2.01 1233 {0.15, 1.7389, 1.7790, 1.8210, 1.9112, 2.0108, 2.1217, 2.2462, 2.3876},
1202 {0.2, 2.3516, 2.4024, 2.4556, 2.5701, 2.697 1234 {0.2, 2.3516, 2.4024, 2.4556, 2.5701, 2.6971, 2.8391, 2.9994, 3.1822},
1203 {0.3, 3.3741, 3.4427, 3.5148, 3.6706, 3.844 1235 {0.3, 3.3741, 3.4427, 3.5148, 3.6706, 3.8445, 4.0404, 4.2631, 4.5193},
1204 {0.4, 4.1788, 4.2620, 4.3496, 4.5398, 4.753 1236 {0.4, 4.1788, 4.2620, 4.3496, 4.5398, 4.7530, 4.9944, 5.2703, 5.5895},
1205 {0.5, 4.8180, 4.9137, 5.0146, 5.2341, 5.481 1237 {0.5, 4.8180, 4.9137, 5.0146, 5.2341, 5.4811, 5.7618, 6.0840, 6.4583},
1206 {0.6, 5.3765, 5.4830, 5.5954, 5.8406, 6.117 1238 {0.6, 5.3765, 5.4830, 5.5954, 5.8406, 6.1173, 6.4326, 6.7958, 7.2191},
1207 {0.7, 5.8208, 5.9369, 6.0596, 6.3275, 6.630 1239 {0.7, 5.8208, 5.9369, 6.0596, 6.3275, 6.6306, 6.9769, 7.3767, 7.8440},
1208 {0.8, 6.2109, 6.3355, 6.4674, 6.7558, 7.082 1240 {0.8, 6.2109, 6.3355, 6.4674, 6.7558, 7.0827, 7.4570, 7.8900, 8.3972},
1209 {1.0, 6.8747, 7.0142, 7.1621, 7.4861, 7.854 1241 {1.0, 6.8747, 7.0142, 7.1621, 7.4861, 7.8546, 8.2778, 8.7690, 9.3464},
1210 {1.2, 7.3933, 7.5454, 7.7068, 8.0612, 8.465 1242 {1.2, 7.3933, 7.5454, 7.7068, 8.0612, 8.4652, 8.9302, 9.4713, 10.1090},
1211 {1.4, 7.8331, 7.9967, 8.1694, 8.5502, 8.985 1243 {1.4, 7.8331, 7.9967, 8.1694, 8.5502, 8.9851, 9.4866, 10.0713, 10.7619},
1212 {1.5, 8.0286, 8.1967, 8.3753, 8.7681, 9.218 1244 {1.5, 8.0286, 8.1967, 8.3753, 8.7681, 9.2181, 9.7352, 10.3399, 11.0546},
1213 {1.7, 8.3813, 8.5585, 8.7469, 9.1618, 9.636 1245 {1.7, 8.3813, 8.5585, 8.7469, 9.1618, 9.6367, 10.1856, 10.8270, 11.5863},
1214 {2.0, 8.8352, 9.0245, 9.2260, 9.6701, 10.17 1246 {2.0, 8.8352, 9.0245, 9.2260, 9.6701, 10.1793, 10.7688, 11.4590, 12.2775},
1215 {3.0, 9.9511, 10.1714, 10.4062, 10.9254, 11 1247 {3.0, 9.9511, 10.1714, 10.4062, 10.9254, 11.5229, 12.2172, 13.0332, 14.0048},
1216 {5.0, 11.3211, 11.5818, 11.8601, 12.4771, 1 1248 {5.0, 11.3211, 11.5818, 11.8601, 12.4771, 13.1898, 14.0213, 15.0024, 16.1752},
1217 {7.0, 11.9480, 12.2357, 12.5432, 13.2260, 1 1249 {7.0, 11.9480, 12.2357, 12.5432, 13.2260, 14.0164, 14.9404, 16.0330, 17.3420}
1218 }; 1250 };
1219 1251
1220 G4double b, bs; 1252 G4double b, bs;
1221 for(i=0; i<nEtaK; ++i) { 1253 for(i=0; i<nEtaK; ++i) {
1222 1254
1223 G4double et = Eta[i]; 1255 G4double et = Eta[i];
1224 G4double loget = G4Log(et); 1256 G4double loget = G4Log(et);
1225 1257
1226 for(j=0; j<nK; ++j) { 1258 for(j=0; j<nK; ++j) {
1227 1259
1228 if(j < 10) { b = bk2[i][10-j]; } 1260 if(j < 10) { b = bk2[i][10-j]; }
1229 else { b = bk1[i][20-j]; } 1261 else { b = bk1[i][20-j]; }
1230 1262
1231 CK[j][i] = SK[j]*loget + TK[j] - b; 1263 CK[j][i] = SK[j]*loget + TK[j] - b;
1232 1264
1233 if(i == nEtaK-1) { 1265 if(i == nEtaK-1) {
1234 ZK[j] = et*(et*et*CK[j][i] - et*UK[j] 1266 ZK[j] = et*(et*et*CK[j][i] - et*UK[j] - VK[j]);
1235 //G4cout << "i= " << i << " j= " << j 1267 //G4cout << "i= " << i << " j= " << j
1236 // << " CK[j][i]= " << CK[j][i 1268 // << " CK[j][i]= " << CK[j][i]
1237 // << " ZK[j]= " << ZK[j] << " 1269 // << " ZK[j]= " << ZK[j] << " b= " << b << G4endl;
1238 } 1270 }
1239 } 1271 }
1240 //G4cout << G4endl; 1272 //G4cout << G4endl;
1241 if(i < nEtaL) { 1273 if(i < nEtaL) {
1242 //G4cout << " LShell:" <<G4endl; 1274 //G4cout << " LShell:" <<G4endl;
1243 for(j=0; j<nL; ++j) { 1275 for(j=0; j<nL; ++j) {
1244 1276
1245 if(j < 8) { 1277 if(j < 8) {
1246 bs = bls3[i][8-j]; 1278 bs = bls3[i][8-j];
1247 b = bll3[i][8-j]; 1279 b = bll3[i][8-j];
1248 } else if(j < 17) { 1280 } else if(j < 17) {
1249 bs = bls2[i][17-j]; 1281 bs = bls2[i][17-j];
1250 b = bll2[i][17-j]; 1282 b = bll2[i][17-j];
1251 } else { 1283 } else {
1252 bs = bls1[i][26-j]; 1284 bs = bls1[i][26-j];
1253 b = bll1[i][26-j]; 1285 b = bll1[i][26-j];
1254 } 1286 }
1255 G4double c = SL[j]*loget + TL[j]; 1287 G4double c = SL[j]*loget + TL[j];
1256 CL[j][i] = c - bs - 3.0*b; 1288 CL[j][i] = c - bs - 3.0*b;
1257 if(i == nEtaL-1) { 1289 if(i == nEtaL-1) {
1258 VL[j] = et*(et*CL[j][i] - UL[j]); 1290 VL[j] = et*(et*CL[j][i] - UL[j]);
1259 //G4cout << "i= " << i << " j= " << 1291 //G4cout << "i= " << i << " j= " << j
1260 // << " CL[j][i]= " << CL[ 1292 // << " CL[j][i]= " << CL[j][i]
1261 // << " VL[j]= " << VL[j] < 1293 // << " VL[j]= " << VL[j] << " b= " << b << " bs= " << bs
1262 // << " et= " << et << G4en 1294 // << " et= " << et << G4endl;
1263 //" UL= " << UL[j] << " TL= " << TL 1295 //" UL= " << UL[j] << " TL= " << TL[j] << " SL= " << SL[j] <<G4endl;
1264 } 1296 }
1265 } 1297 }
1266 //G4cout << G4endl; 1298 //G4cout << G4endl;
1267 } 1299 }
1268 } 1300 }
1269 1301
1270 const G4double xzk[34] = { 11.7711, << 1302 static const G4double xzk[34] = { 11.7711,
1271 13.3669, 15.5762, 17.1715, 18.7667, 20.85 1303 13.3669, 15.5762, 17.1715, 18.7667, 20.8523, 23.0606, 24.901, 26.9861, 29.4394, 31.77,
1272 34.3457, 37.4119, 40.3555, 42.3177, 44.77 1304 34.3457, 37.4119, 40.3555, 42.3177, 44.7705, 47.2234, 50.78, 53.8458, 56.4214, 58.3834,
1273 60.9586, 63.6567, 66.5998, 68.807, 71.872 1305 60.9586, 63.6567, 66.5998, 68.807, 71.8728, 74.5706, 77.3911, 81.8056, 85.7297, 89.8988,
1274 93.4549, 96.2753 1306 93.4549, 96.2753, 99.709};
1275 const G4double yzk[34] = { 0.70663, << 1307 static const G4double yzk[34] = { 0.70663,
1276 0.72033, 0.73651, 0.74647, 0.75518, 0.763 1308 0.72033, 0.73651, 0.74647, 0.75518, 0.76388, 0.77258, 0.78129, 0.78625, 0.7937, 0.79991,
1277 0.80611, 0.8123, 0.8185, 0.82097, 0.82467 1309 0.80611, 0.8123, 0.8185, 0.82097, 0.82467, 0.82838, 0.83457, 0.83702, 0.84198, 0.8432,
1278 0.84565, 0.84936, 0.85181, 0.85303, 0.855 1310 0.84565, 0.84936, 0.85181, 0.85303, 0.85548, 0.85794, 0.8604, 0.86283, 0.86527, 0.86646,
1279 0.86891, 0.87011 1311 0.86891, 0.87011, 0.87381};
1280 1312
1281 const G4double xzl[36] = { 15.5102, << 1313 static const G4double xzl[36] = { 15.5102,
1282 16.7347, 17.9592, 19.551, 21.0204, 22.612 1314 16.7347, 17.9592, 19.551, 21.0204, 22.6122, 24.9388, 27.3878, 29.5918, 31.3061, 32.898,
1283 34.4898, 36.2041, 38.4082, 40.3674, 42.57 1315 34.4898, 36.2041, 38.4082, 40.3674, 42.5714, 44.898, 47.4694, 49.9184, 52.7347, 55.9184,
1284 59.3469, 61.9184, 64.6122, 67.4286, 71.46 1316 59.3469, 61.9184, 64.6122, 67.4286, 71.4694, 75.2653, 78.3265, 81.2653, 85.551, 88.7347,
1285 91.551, 94.2449, 1317 91.551, 94.2449, 96.449, 98.4082, 99.7551};
1286 const G4double yzl[36] = { 0.29875, << 1318 static const G4double yzl[36] = { 0.29875,
1287 0.31746, 0.33368, 0.35239, 0.36985, 0.387 1319 0.31746, 0.33368, 0.35239, 0.36985, 0.38732, 0.41102, 0.43472, 0.45343, 0.4659, 0.47713,
1288 0.4896, 0.50083, 0.51331, 0.52328, 0.5307 1320 0.4896, 0.50083, 0.51331, 0.52328, 0.53077, 0.54075, 0.54823, 0.55572, 0.56445, 0.57193,
1289 0.58191, 0.5869, 0.59189, 0.60062, 0.6068 1321 0.58191, 0.5869, 0.59189, 0.60062, 0.60686, 0.61435, 0.61809, 0.62183, 0.62931, 0.6343,
1290 0.6368, 0.64054 1322 0.6368, 0.64054, 0.64304, 0.64428, 0.64678};
1291 1323
1292 sThetaK = new G4PhysicsFreeVector(34, false 1324 sThetaK = new G4PhysicsFreeVector(34, false);
1293 for(i=0; i<34; ++i) { sThetaK->PutValues(i, 1325 for(i=0; i<34; ++i) { sThetaK->PutValues(i, xzk[i], yzk[i]); }
1294 sThetaL = new G4PhysicsFreeVector(36, false 1326 sThetaL = new G4PhysicsFreeVector(36, false);
1295 for(i=0; i<36; ++i) { sThetaL->PutValues(i, 1327 for(i=0; i<36; ++i) { sThetaL->PutValues(i, xzl[i], yzl[i]); }
1296 } 1328 }
1297 1329
1298 //....oooOO0OOooo........oooOO0OOooo........o 1330 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1299 1331
1300 1332