| 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 // $Id: G4EmCorrections.cc,v 1.64 2010-12-02 12:19:40 vnivanch Exp $
>> 27 // GEANT4 tag $Name: not supported by cvs2svn $
26 // 28 //
27 // ------------------------------------------- 29 // -------------------------------------------------------------------
28 // 30 //
29 // GEANT4 Class 31 // GEANT4 Class
30 // 32 //
31 // File name: G4EmCorrections 33 // File name: G4EmCorrections
32 // 34 //
33 // Author: Vladimir Ivanchenko 35 // Author: Vladimir Ivanchenko
34 // 36 //
35 // Creation date: 13.01.2005 37 // Creation date: 13.01.2005
36 // 38 //
37 // Modifications: 39 // Modifications:
38 // 05.05.2005 V.Ivanchenko Fix misprint in Mot 40 // 05.05.2005 V.Ivanchenko Fix misprint in Mott term
39 // 26.11.2005 V.Ivanchenko Fix effective charg 41 // 26.11.2005 V.Ivanchenko Fix effective charge for heavy ions using original paper
40 // 28.04.2006 V.Ivanchenko General cleanup, ad 42 // 28.04.2006 V.Ivanchenko General cleanup, add finite size corrections
41 // 13.05.2006 V.Ivanchenko Add corrections for 43 // 13.05.2006 V.Ivanchenko Add corrections for ion stopping
42 // 08.05.2007 V.Ivanchenko Use G4IonTable for 44 // 08.05.2007 V.Ivanchenko Use G4IonTable for ion mass instead of NistTable to avoid
43 // division by zero 45 // division by zero
44 // 29.02.2008 V.Ivanchenko use expantions for 46 // 29.02.2008 V.Ivanchenko use expantions for log and power function
45 // 21.04.2008 Updated computations for ions (V 47 // 21.04.2008 Updated computations for ions (V.Ivanchenko)
46 // 20.05.2008 Removed Finite Size correction ( 48 // 20.05.2008 Removed Finite Size correction (V.Ivanchenko)
47 // 19.04.2012 Fix reproducibility problem (A.R <<
48 // 49 //
49 // 50 //
50 // Class Description: 51 // Class Description:
51 // 52 //
52 // This class provides calculation of EM corre 53 // This class provides calculation of EM corrections to ionisation
53 // 54 //
54 55
55 // ------------------------------------------- 56 // -------------------------------------------------------------------
56 // 57 //
57 58
58 #include "G4EmCorrections.hh" 59 #include "G4EmCorrections.hh"
59 #include "G4PhysicalConstants.hh" << 60 #include "Randomize.hh"
60 #include "G4SystemOfUnits.hh" <<
61 #include "G4ParticleTable.hh" 61 #include "G4ParticleTable.hh"
62 #include "G4IonTable.hh" 62 #include "G4IonTable.hh"
63 #include "G4VEmModel.hh" 63 #include "G4VEmModel.hh"
64 #include "G4Proton.hh" 64 #include "G4Proton.hh"
65 #include "G4GenericIon.hh" 65 #include "G4GenericIon.hh"
>> 66 #include "G4LPhysicsFreeVector.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 "G4LPhysicsFreeVector.hh"
71 #include "G4Exp.hh" <<
72 #include "G4Pow.hh" <<
73 72
74 //....oooOO0OOooo........oooOO0OOooo........oo 73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75 74
76 namespace << 75 G4EmCorrections::G4EmCorrections()
77 { 76 {
78 constexpr G4double inveplus = 1.0/CLHEP::epl << 77 particle = 0;
79 constexpr G4double alpha2 = CLHEP::fine_stru << 78 curParticle= 0;
80 } << 79 material = 0;
81 const G4double G4EmCorrections::ZD[11] = << 80 curMaterial= 0;
82 {0., 0., 0., 1.72, 2.09, 2.48, 2.82, 3.16, << 81 curVector = 0;
83 const G4double G4EmCorrections::UK[20] = {1.99 << 82 kinEnergy = 0.0;
84 2.0817, 2.0945, 2.0 << 83 ionLEModel = 0;
85 2.1197, 2.1246, 2.1 << 84 ionHEModel = 0;
86 2.1310, 2.1310, 2.1 << 85 nIons = 0;
87 const G4double G4EmCorrections::VK[20] = {8.34 << 86 verbose = 1;
88 8.3219, 8.3201, 8.3 << 87 ncouples = 0;
89 8.3191, 8.3199, 8.3 << 88 massFactor = 1.0;
90 8.3244, 8.3264, 8.3 << 89 eth = 2.0*MeV;
91 G4double G4EmCorrections::ZK[] = {0.0}; << 90 nbinCorr = 20;
92 const G4double G4EmCorrections::Eta[29] = {0.0 << 91 eCorrMin = 25.*keV;
93 0.03, 0.04, 0.05 << 92 eCorrMax = 250.*MeV;
94 0.1, 0.15, 0.2, << 93 nist = G4NistManager::Instance();
95 0.5, 0.6, 0.7, <<
96 1.2, 1.4, 1.5, <<
97 G4double G4EmCorrections::CK[20][29]; <<
98 G4double G4EmCorrections::CL[26][28]; <<
99 const G4double G4EmCorrections::UL[] = {0.1215 <<
100 1.4379, 1.5032, 1.5 <<
101 1.8036, 1.8543, 1.8 <<
102 1.9508, 1.9696, 1.9 <<
103 2.0001, 2.0039, 2.0 <<
104 G4double G4EmCorrections::VL[] = {0.0}; <<
105 G4double G4EmCorrections::sWmaxBarkas = 10.0; <<
106 <<
107 G4PhysicsFreeVector* G4EmCorrections::sBarkasC <<
108 G4PhysicsFreeVector* G4EmCorrections::sThetaK <<
109 G4PhysicsFreeVector* G4EmCorrections::sThetaL <<
110 <<
111 G4EmCorrections::G4EmCorrections(G4int verb) <<
112 : verbose(verb) <<
113 { <<
114 eth = 2.0*CLHEP::MeV; <<
115 eCorrMin = 25.*CLHEP::keV; <<
116 eCorrMax = 1.*CLHEP::GeV; <<
117 <<
118 ionTable = G4ParticleTable::GetParticleTable 94 ionTable = G4ParticleTable::GetParticleTable()->GetIonTable();
119 g4calc = G4Pow::GetInstance(); << 95 BarkasCorr = ThetaK = ThetaL = 0;
120 << 96 Initialise();
121 // fill vectors <<
122 if (nullptr == sBarkasCorr) { <<
123 Initialise(); <<
124 isInitializer = true; <<
125 } <<
126 } 97 }
127 98
128 //....oooOO0OOooo........oooOO0OOooo........oo 99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
129 100
130 G4EmCorrections::~G4EmCorrections() 101 G4EmCorrections::~G4EmCorrections()
131 { 102 {
132 for (G4int i=0; i<nIons; ++i) { delete stopD << 103 for(G4int i=0; i<nIons; ++i) {delete stopData[i];}
133 if (isInitializer) { << 104 delete BarkasCorr;
134 delete sBarkasCorr; << 105 delete ThetaK;
135 delete sThetaK; << 106 delete ThetaL;
136 delete sThetaL; <<
137 sBarkasCorr = sThetaK = sThetaL = nullptr; <<
138 } <<
139 } <<
140 <<
141 void G4EmCorrections::SetupKinematics(const G4 <<
142 const G4Material* mat, <<
143 const G4double kineticEnergy) <<
144 { <<
145 if(kineticEnergy != kinEnergy || p != partic <<
146 particle = p; <<
147 kinEnergy = kineticEnergy; <<
148 mass = p->GetPDGMass(); <<
149 tau = kineticEnergy / mass; <<
150 gamma = 1.0 + tau; <<
151 bg2 = tau * (tau+2.0); <<
152 beta2 = bg2/(gamma*gamma); <<
153 beta = std::sqrt(beta2); <<
154 ba2 = beta2/alpha2; <<
155 const G4double ratio = CLHEP::electron_mas <<
156 tmax = 2.0*CLHEP::electron_mass_c2*bg2 <<
157 /(1. + 2.0*gamma*ratio + ratio*ratio); <<
158 charge = p->GetPDGCharge()*inveplus; <<
159 if(charge > 1.5) { charge = effCharge.Effe <<
160 q2 = charge*charge; <<
161 } <<
162 if(mat != material) { <<
163 material = mat; <<
164 theElementVector = material->GetElementVec <<
165 atomDensity = material->GetAtomicNumDensi <<
166 numberOfElements = (G4int)material->GetNum <<
167 } <<
168 } 107 }
169 108
170 //....oooOO0OOooo........oooOO0OOooo........oo 109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
171 110
172 G4double G4EmCorrections::HighOrderCorrections 111 G4double G4EmCorrections::HighOrderCorrections(const G4ParticleDefinition* p,
173 112 const G4Material* mat,
174 << 113 G4double e, G4double)
175 { 114 {
176 // . Z^3 Barkas effect in the stopping power 115 // . Z^3 Barkas effect in the stopping power of matter for charged particles
177 // J.C Ashley and R.H.Ritchie 116 // J.C Ashley and R.H.Ritchie
178 // Physical review B Vol.5 No.7 1 April 19 117 // Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
179 // and ICRU49 report 118 // and ICRU49 report
180 // valid for kineticEnergy < 0.5 MeV 119 // valid for kineticEnergy < 0.5 MeV
181 // Other corrections from S.P.Ahlen Rev. M 120 // Other corrections from S.P.Ahlen Rev. Mod. Phys., Vol 52, No1, 1980
182 121
183 SetupKinematics(p, mat, e); 122 SetupKinematics(p, mat, e);
184 if(tau <= 0.0) { return 0.0; } 123 if(tau <= 0.0) { return 0.0; }
185 124
186 const G4double Barkas = BarkasCorrection(p, << 125 G4double Barkas = BarkasCorrection (p, mat, e);
187 const G4double Bloch = BlochCorrection(p, m << 126 G4double Bloch = BlochCorrection (p, mat, e);
188 const G4double Mott = MottCorrection(p, mat, << 127 G4double Mott = MottCorrection (p, mat, e);
189 128
190 G4double sum = 2.0*(Barkas + Bloch) + Mott; << 129 G4double sum = (2.0*(Barkas + Bloch) + Mott);
191 130
192 if(verbose > 1) { 131 if(verbose > 1) {
193 G4cout << "EmCorrections: E(MeV)= " << e/M 132 G4cout << "EmCorrections: E(MeV)= " << e/MeV << " Barkas= " << Barkas
194 << " Bloch= " << Bloch << " Mott= " << 133 << " Bloch= " << Bloch << " Mott= " << Mott
195 << " Sum= " << sum << " q2= " << q2 << 134 << " Sum= " << sum << " q2= " << q2 << G4endl;
196 G4cout << " ShellCorrection: " << ShellCor 135 G4cout << " ShellCorrection: " << ShellCorrection(p, mat, e)
197 << " Kshell= " << KShellCorrection( << 136 << " Kshell= " << KShellCorrection(p, mat, e)
198 << " Lshell= " << LShellCorrection( << 137 << " Lshell= " << LShellCorrection(p, mat, e)
199 << " " << mat->GetName() << G4end << 138 << " " << mat->GetName() << G4endl;
200 } 139 }
201 sum *= material->GetElectronDensity()*q2*CLH << 140 sum *= material->GetElectronDensity() * q2 * twopi_mc2_rcl2 /beta2;
202 return sum; 141 return sum;
203 } 142 }
204 143
205 //....oooOO0OOooo........oooOO0OOooo........oo 144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
206 145
207 G4double G4EmCorrections::IonBarkasCorrection( 146 G4double G4EmCorrections::IonBarkasCorrection(const G4ParticleDefinition* p,
208 << 147 const G4Material* mat,
209 << 148 G4double e)
210 { 149 {
>> 150 // . Z^3 Barkas effect in the stopping power of matter for charged particles
>> 151 // J.C Ashley and R.H.Ritchie
>> 152 // Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
>> 153 // and ICRU49 report
>> 154 // valid for kineticEnergy < 0.5 MeV
>> 155
211 SetupKinematics(p, mat, e); 156 SetupKinematics(p, mat, e);
212 return 2.0*BarkasCorrection(p, mat, e, true) << 157 G4double res = 0.0;
213 material->GetElectronDensity() * q2 * CLHE << 158 if(tau > 0.0)
>> 159 res = 2.0*BarkasCorrection(p, mat, e)*
>> 160 material->GetElectronDensity() * q2 * twopi_mc2_rcl2 /beta2;
>> 161 return res;
214 } 162 }
215 163
216 //....oooOO0OOooo........oooOO0OOooo........oo 164 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
217 165
218 G4double G4EmCorrections::ComputeIonCorrection 166 G4double G4EmCorrections::ComputeIonCorrections(const G4ParticleDefinition* p,
219 << 167 const G4Material* mat,
220 << 168 G4double e)
221 { 169 {
222 // . Z^3 Barkas effect in the stopping power 170 // . Z^3 Barkas effect in the stopping power of matter for charged particles
223 // J.C Ashley and R.H.Ritchie 171 // J.C Ashley and R.H.Ritchie
224 // Physical review B Vol.5 No.7 1 April 19 172 // Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
225 // and ICRU49 report 173 // and ICRU49 report
226 // valid for kineticEnergy < 0.5 MeV 174 // valid for kineticEnergy < 0.5 MeV
227 // Other corrections from S.P.Ahlen Rev. M 175 // Other corrections from S.P.Ahlen Rev. Mod. Phys., Vol 52, No1, 1980
228 SetupKinematics(p, mat, e); 176 SetupKinematics(p, mat, e);
229 if(tau <= 0.0) { return 0.0; } 177 if(tau <= 0.0) { return 0.0; }
230 178
231 const G4double Barkas = BarkasCorrection (p, << 179 G4double Barkas = BarkasCorrection (p, mat, e);
232 const G4double Bloch = BlochCorrection (p, << 180 G4double Bloch = BlochCorrection (p, mat, e);
233 const G4double Mott = MottCorrection (p, mat << 181 G4double Mott = MottCorrection (p, mat, e);
234 182
235 G4double sum = 2.0*(Barkas*(charge - 1.0)/ch 183 G4double sum = 2.0*(Barkas*(charge - 1.0)/charge + Bloch) + Mott;
236 184
237 if(verbose > 1) { 185 if(verbose > 1) {
238 G4cout << "EmCorrections: E(MeV)= " << e/M 186 G4cout << "EmCorrections: E(MeV)= " << e/MeV << " Barkas= " << Barkas
239 << " Bloch= " << Bloch << " Mott= " << 187 << " Bloch= " << Bloch << " Mott= " << Mott
240 << " Sum= " << sum << G4endl; << 188 << " Sum= " << sum << G4endl;
241 } 189 }
242 sum *= material->GetElectronDensity() * q2 * << 190 sum *= material->GetElectronDensity() * q2 * twopi_mc2_rcl2 /beta2;
243 191
244 if(verbose > 1) { G4cout << " Sum= " << sum 192 if(verbose > 1) { G4cout << " Sum= " << sum << G4endl; }
245 return sum; 193 return sum;
246 } 194 }
247 195
248 //....oooOO0OOooo........oooOO0OOooo........oo 196 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
249 197
250 G4double G4EmCorrections::IonHighOrderCorrecti 198 G4double G4EmCorrections::IonHighOrderCorrections(const G4ParticleDefinition* p,
251 << 199 const G4MaterialCutsCouple* couple,
252 << 200 G4double e)
253 { 201 {
254 // . Z^3 Barkas effect in the stopping power << 202 // . Z^3 Barkas effect in the stopping power of matter for charged particles
255 // J.C Ashley and R.H.Ritchie << 203 // J.C Ashley and R.H.Ritchie
256 // Physical review B Vol.5 No.7 1 April 19 << 204 // Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
257 // and ICRU49 report << 205 // and ICRU49 report
258 // valid for kineticEnergy < 0.5 MeV << 206 // valid for kineticEnergy < 0.5 MeV
259 // Other corrections from S.P.Ahlen Rev. M << 207 // Other corrections from S.P.Ahlen Rev. Mod. Phys., Vol 52, No1, 1980
260 208
261 G4double sum = 0.0; 209 G4double sum = 0.0;
262 210
263 if (nullptr != ionHEModel) { << 211 if(ionHEModel) {
264 G4int Z = G4lrint(p->GetPDGCharge()*invepl << 212 G4int Z = G4int(p->GetPDGCharge()/eplus + 0.5);
265 Z = std::max(std::min(Z, 99), 1); << 213 if(Z >= 100) Z = 99;
266 << 214 else if(Z < 1) Z = 1;
267 const G4double ethscaled = eth*p->GetPDGMa << 215
268 const G4int ionPDG = p->GetPDGEncoding(); << 216 // fill vector
269 auto iter = thcorr.find(ionPDG); << 217 if(thcorr[Z].size() == 0) {
270 if (iter == thcorr.end()) { // Not found: << 218 thcorr[Z].resize(ncouples);
271 std::vector<G4double> v; << 219 G4double ethscaled = eth*p->GetPDGMass()/proton_mass_c2;
272 for(std::size_t i=0; i<ncouples; ++i){ << 220
273 v.push_back(ethscaled*ComputeIonCorrec << 221 for(size_t i=0; i<ncouples; ++i) {
274 } << 222 (thcorr[Z])[i] = ethscaled*ComputeIonCorrections(p, currmat[i], ethscaled);
275 thcorr.insert(std::pair< G4int, std::vec << 223 //G4cout << i << ". ethscaled= " << ethscaled
276 } << 224 //<< " corr= " << (thcorr[Z])[i]/ethscaled << G4endl;
277 G4double rest = 0.0; << 225 }
278 iter = thcorr.find(ionPDG); << 226 }
279 if (iter != thcorr.end()) { rest = (iter-> << 227 G4double rest = (thcorr[Z])[couple->GetIndex()];
280 228
281 sum = ComputeIonCorrections(p,couple->GetM 229 sum = ComputeIonCorrections(p,couple->GetMaterial(),e) - rest/e;
282 230
283 if(verbose > 1) { << 231 if(verbose > 1) { G4cout << " Sum= " << sum << " dSum= " << rest/e << G4endl; }
284 G4cout << " Sum= " << sum << " dSum= " < <<
285 } <<
286 } 232 }
287 return sum; 233 return sum;
288 } 234 }
289 235
290 //....oooOO0OOooo........oooOO0OOooo........oo 236 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
291 237
292 G4double G4EmCorrections::Bethe(const G4Partic 238 G4double G4EmCorrections::Bethe(const G4ParticleDefinition* p,
293 const G4Materi << 239 const G4Material* mat,
294 const G4double << 240 G4double e)
295 { 241 {
296 SetupKinematics(p, mat, e); 242 SetupKinematics(p, mat, e);
297 const G4double eexc = material->GetIonisati << 243 G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
298 const G4double eexc2 = eexc*eexc; << 244 G4double eexc2 = eexc*eexc;
299 return 0.5*G4Log(2.0*electron_mass_c2*bg2*tm << 245 G4double dedx = 0.5*std::log(2.0*electron_mass_c2*bg2*tmax/eexc2)-beta2;
>> 246 return dedx;
300 } 247 }
301 248
302 //....oooOO0OOooo........oooOO0OOooo........oo 249 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
303 250
304 G4double G4EmCorrections::SpinCorrection(const 251 G4double G4EmCorrections::SpinCorrection(const G4ParticleDefinition* p,
305 const << 252 const G4Material* mat,
306 const << 253 G4double e)
307 { 254 {
308 SetupKinematics(p, mat, e); 255 SetupKinematics(p, mat, e);
309 const G4double dedx = 0.5*tmax/(kinEnergy + << 256 G4double dedx = 0.5*tmax/(kinEnergy + mass);
310 return 0.5*dedx*dedx; 257 return 0.5*dedx*dedx;
311 } 258 }
312 259
313 //....oooOO0OOooo........oooOO0OOooo........oo 260 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
314 261
315 G4double G4EmCorrections:: KShellCorrection(co 262 G4double G4EmCorrections:: KShellCorrection(const G4ParticleDefinition* p,
316 co << 263 const G4Material* mat,
317 co << 264 G4double e)
318 { 265 {
319 SetupKinematics(p, mat, e); 266 SetupKinematics(p, mat, e);
320 G4double term = 0.0; 267 G4double term = 0.0;
321 for (G4int i = 0; i<numberOfElements; ++i) { 268 for (G4int i = 0; i<numberOfElements; ++i) {
322 269
323 G4double Z = (*theElementVector)[i]->GetZ( 270 G4double Z = (*theElementVector)[i]->GetZ();
324 G4int iz = (*theElementVector)[i]->GetZa << 271 G4int iz = G4int(Z);
325 G4double f = 1.0; 272 G4double f = 1.0;
326 G4double Z2= (Z-0.3)*(Z-0.3); 273 G4double Z2= (Z-0.3)*(Z-0.3);
327 if(1 == iz) { 274 if(1 == iz) {
328 f = 0.5; 275 f = 0.5;
329 Z2 = 1.0; 276 Z2 = 1.0;
330 } 277 }
331 const G4double eta = ba2/Z2; << 278 G4double eta = ba2/Z2;
332 const G4double tet = (11 < iz) ? sThetaK-> << 279 G4double tet = Z2*(1. + Z2*0.25*alpha2);
>> 280 if(11 < iz) { tet = ThetaK->Value(Z); }
333 term += f*atomDensity[i]*KShell(tet,eta)/Z 281 term += f*atomDensity[i]*KShell(tet,eta)/Z;
334 } 282 }
335 283
336 term /= material->GetTotNbOfAtomsPerVolume() 284 term /= material->GetTotNbOfAtomsPerVolume();
337 285
338 return term; 286 return term;
339 } 287 }
340 288
341 //....oooOO0OOooo........oooOO0OOooo........oo 289 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
342 290
343 G4double G4EmCorrections:: LShellCorrection(co 291 G4double G4EmCorrections:: LShellCorrection(const G4ParticleDefinition* p,
344 co << 292 const G4Material* mat,
345 co << 293 G4double e)
346 { 294 {
347 SetupKinematics(p, mat, e); 295 SetupKinematics(p, mat, e);
348 G4double term = 0.0; 296 G4double term = 0.0;
349 for (G4int i = 0; i<numberOfElements; ++i) { 297 for (G4int i = 0; i<numberOfElements; ++i) {
350 298
351 const G4double Z = (*theElementVector)[i]- << 299 G4double Z = (*theElementVector)[i]->GetZ();
352 const G4int iz = (*theElementVector)[i]->G << 300 G4int iz = G4int(Z);
353 if(2 < iz) { 301 if(2 < iz) {
354 const G4double Zeff = (iz < 10) ? Z - ZD << 302 G4double Zeff = Z - ZD[10];
355 const G4double Z2= Zeff*Zeff; << 303 if(iz < 10) { Zeff = Z - ZD[iz]; }
356 const G4double eta = ba2/Z2; << 304 G4double Z2= Zeff*Zeff;
357 G4double tet = sThetaL->Value(Z); << 305 G4double f = 0.125;
358 G4int nmax = std::min(4, G4AtomicShells: << 306 G4double eta = ba2/Z2;
359 for (G4int j=1; j<nmax; ++j) { << 307 G4double tet = ThetaL->Value(Z);
>> 308 G4int nmax = std::min(4,G4AtomicShells::GetNumberOfShells(iz));
>> 309 for(G4int j=1; j<nmax; ++j) {
360 G4int ne = G4AtomicShells::GetNumberOf 310 G4int ne = G4AtomicShells::GetNumberOfElectrons(iz,j);
361 if (15 >= iz) { << 311 if(15 >= iz) {
362 tet = (3 > j) ? 0.25*Z2*(1.0 + 5*Z2* << 312 if(3 > j) { tet = 0.25*Z2*(1.0 + 5*Z2*alpha2/16.); }
363 0.25*Z2*(1.0 + Z2*alpha2/16.); << 313 else { tet = 0.25*Z2*(1.0 + Z2*alpha2/16.); }
364 } << 314 }
365 //G4cout << " LShell: j= " << j << " n << 315 //G4cout << " LShell: j= " << j << " ne= " << ne << " e(eV)= " << e/eV
366 // << " ThetaL= " << tet << G4en << 316 // << " ThetaL= " << tet << G4endl;
367 term += 0.125*ne*atomDensity[i]*LShell << 317 term += f*ne*atomDensity[i]*LShell(tet,eta)/Z;
368 } 318 }
369 } 319 }
370 } 320 }
371 321
372 term /= material->GetTotNbOfAtomsPerVolume() 322 term /= material->GetTotNbOfAtomsPerVolume();
373 323
374 return term; 324 return term;
375 } 325 }
376 326
377 //....oooOO0OOooo........oooOO0OOooo........oo 327 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
378 328
379 G4double G4EmCorrections::KShell(const G4doubl << 329 G4double G4EmCorrections::KShell(G4double tet, G4double eta)
380 { 330 {
381 G4double corr = 0.0; 331 G4double corr = 0.0;
382 332
383 static const G4double TheK[20] = <<
384 {0.64, 0.65, 0.66, 0.68, 0.70, 0.72, 0.74, <<
385 0.80, 0.82, 0.84, 0.85, 0.86, 0.88, 0.90, <<
386 <<
387 <<
388 G4double x = tet; 333 G4double x = tet;
389 G4int itet = 0; 334 G4int itet = 0;
390 G4int ieta = 0; 335 G4int ieta = 0;
391 if(tet < TheK[0]) { 336 if(tet < TheK[0]) {
392 x = TheK[0]; 337 x = TheK[0];
393 } else if(tet > TheK[nK-1]) { 338 } else if(tet > TheK[nK-1]) {
394 x = TheK[nK-1]; 339 x = TheK[nK-1];
395 itet = nK-2; 340 itet = nK-2;
396 } else { 341 } else {
397 itet = Index(x, TheK, nK); 342 itet = Index(x, TheK, nK);
398 } 343 }
399 // asymptotic case << 344 // assimptotic case
400 if(eta >= Eta[nEtaK-1]) { 345 if(eta >= Eta[nEtaK-1]) {
401 corr = << 346 corr = (Value(x, TheK[itet], TheK[itet+1], UK[itet], UK[itet+1]) +
402 (Value(x, TheK[itet], TheK[itet+1], UK[i << 347 Value(x, TheK[itet], TheK[itet+1], VK[itet], VK[itet+1])/eta +
403 Value(x, TheK[itet], TheK[itet+1], VK[i << 348 Value(x, TheK[itet], TheK[itet+1], ZK[itet], ZK[itet+1])/(eta*eta))/eta;
404 Value(x, TheK[itet], TheK[itet+1], ZK[i <<
405 } else { 349 } else {
406 G4double y = eta; 350 G4double y = eta;
407 if(eta < Eta[0]) { 351 if(eta < Eta[0]) {
408 y = Eta[0]; << 352 y = Eta[0];
409 } else { 353 } else {
410 ieta = Index(y, Eta, nEtaK); 354 ieta = Index(y, Eta, nEtaK);
411 } 355 }
412 corr = Value2(x, y, TheK[itet], TheK[itet+ 356 corr = Value2(x, y, TheK[itet], TheK[itet+1], Eta[ieta], Eta[ieta+1],
413 CK[itet][ieta], CK[itet+1][i 357 CK[itet][ieta], CK[itet+1][ieta],
414 CK[itet][ieta+1], CK[itet+1] << 358 CK[itet][ieta+1], CK[itet+1][ieta+1]);
415 //G4cout << " x= " <<x<<" y= "<<y<<" tet 359 //G4cout << " x= " <<x<<" y= "<<y<<" tet= " <<TheK[itet]
416 // <<" "<< TheK[itet+1]<<" eta= << 360 // <<" "<< TheK[itet+1]<<" eta= "<< Eta[ieta]<<" "<< Eta[ieta+1]
417 // <<" CK= " << CK[itet][ieta]<< << 361 // <<" CK= " << CK[itet][ieta]<<" "<< CK[itet+1][ieta]
418 // <<" "<< CK[itet][ieta+1]<<" " << 362 // <<" "<< CK[itet][ieta+1]<<" "<< CK[itet+1][ieta+1]<<G4endl;
419 } 363 }
420 //G4cout << "Kshell: tet= " << tet << " eta 364 //G4cout << "Kshell: tet= " << tet << " eta= " << eta << " C= " << corr
421 // << " itet= " << itet << " ieta= << 365 // << " itet= " << itet << " ieta= " << ieta <<G4endl;
422 return corr; 366 return corr;
423 } 367 }
424 368
425 //....oooOO0OOooo........oooOO0OOooo........oo 369 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
426 370
427 G4double G4EmCorrections::LShell(const G4doubl << 371 G4double G4EmCorrections::LShell(G4double tet, G4double eta)
428 { 372 {
429 G4double corr = 0.0; 373 G4double corr = 0.0;
430 374
431 static const G4double TheL[26] = <<
432 {0.24, 0.26, 0.28, 0.30, 0.32, 0.34, 0.3 <<
433 0.42, 0.44, 0.45, 0.46, 0.48, 0.50, 0.5 <<
434 0.58, 0.60, 0.62, 0.64, 0.65, 0.66}; <<
435 <<
436 G4double x = tet; 375 G4double x = tet;
437 G4int itet = 0; 376 G4int itet = 0;
438 G4int ieta = 0; 377 G4int ieta = 0;
439 if(tet < TheL[0]) { 378 if(tet < TheL[0]) {
440 x = TheL[0]; 379 x = TheL[0];
441 } else if(tet > TheL[nL-1]) { 380 } else if(tet > TheL[nL-1]) {
442 x = TheL[nL-1]; 381 x = TheL[nL-1];
443 itet = nL-2; 382 itet = nL-2;
444 } else { 383 } else {
445 itet = Index(x, TheL, nL); 384 itet = Index(x, TheL, nL);
446 } 385 }
447 386
448 // asymptotic case << 387 // assimptotic case
449 if(eta >= Eta[nEtaL-1]) { 388 if(eta >= Eta[nEtaL-1]) {
450 corr = (Value(x, TheL[itet], TheL[itet+1], 389 corr = (Value(x, TheL[itet], TheL[itet+1], UL[itet], UL[itet+1])
451 + Value(x, TheL[itet], TheL[itet << 390 + Value(x, TheL[itet], TheL[itet+1], VL[itet], VL[itet+1])/eta
452 )/eta; << 391 )/eta;
453 } else { 392 } else {
454 G4double y = eta; 393 G4double y = eta;
455 if(eta < Eta[0]) { 394 if(eta < Eta[0]) {
456 y = Eta[0]; 395 y = Eta[0];
457 } else { 396 } else {
458 ieta = Index(y, Eta, nEtaL); 397 ieta = Index(y, Eta, nEtaL);
459 } 398 }
460 corr = Value2(x, y, TheL[itet], TheL[itet+ 399 corr = Value2(x, y, TheL[itet], TheL[itet+1], Eta[ieta], Eta[ieta+1],
461 CL[itet][ieta], CL[itet+1][i << 400 CL[itet][ieta], CL[itet+1][ieta],
462 CL[itet][ieta+1], CL[itet+1] << 401 CL[itet][ieta+1], CL[itet+1][ieta+1]);
463 //G4cout << " x= " <<x<<" y= "<<y<<" tet 402 //G4cout << " x= " <<x<<" y= "<<y<<" tet= " <<TheL[itet]
464 // <<" "<< TheL[itet+1]<<" eta= << 403 // <<" "<< TheL[itet+1]<<" eta= "<< Eta[ieta]<<" "<< Eta[ieta+1]
465 // <<" CL= " << CL[itet][ieta]<< << 404 // <<" CL= " << CL[itet][ieta]<<" "<< CL[itet+1][ieta]
466 // <<" "<< CL[itet][ieta+1]<<" " << 405 // <<" "<< CL[itet][ieta+1]<<" "<< CL[itet+1][ieta+1]<<G4endl;
467 } 406 }
468 //G4cout<<"Lshell: tet= "<<tet<<" eta= "<<e 407 //G4cout<<"Lshell: tet= "<<tet<<" eta= "<<eta<<" itet= "<<itet
469 // <<" ieta= "<<ieta<<" Corr= "<<cor << 408 // <<" ieta= "<<ieta<<" Corr= "<<corr<<G4endl;
470 return corr; 409 return corr;
471 } 410 }
472 411
473 //....oooOO0OOooo........oooOO0OOooo........oo 412 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
474 413
475 G4double G4EmCorrections::ShellCorrectionSTD(c 414 G4double G4EmCorrections::ShellCorrectionSTD(const G4ParticleDefinition* p,
476 c << 415 const G4Material* mat,
477 c << 416 G4double e)
478 { 417 {
479 SetupKinematics(p, mat, e); 418 SetupKinematics(p, mat, e);
480 G4double taulim= 8.0*MeV/mass; 419 G4double taulim= 8.0*MeV/mass;
481 G4double bg2lim= taulim * (taulim+2.0); 420 G4double bg2lim= taulim * (taulim+2.0);
482 421
483 G4double* shellCorrectionVector = 422 G4double* shellCorrectionVector =
484 material->GetIonisation()->GetShel 423 material->GetIonisation()->GetShellCorrectionVector();
485 G4double sh = 0.0; 424 G4double sh = 0.0;
486 G4double x = 1.0; 425 G4double x = 1.0;
487 G4double taul = material->GetIonisation()-> 426 G4double taul = material->GetIonisation()->GetTaul();
488 427
489 if ( bg2 >= bg2lim ) { 428 if ( bg2 >= bg2lim ) {
490 for (G4int k=0; k<3; ++k) { << 429 for (G4int k=0; k<3; k++) {
491 x *= bg2 ; 430 x *= bg2 ;
492 sh += shellCorrectionVector[k]/x; 431 sh += shellCorrectionVector[k]/x;
493 } 432 }
494 433
495 } else { 434 } else {
496 for (G4int k=0; k<3; ++k) { << 435 for (G4int k=0; k<3; k++) {
497 x *= bg2lim ; 436 x *= bg2lim ;
498 sh += shellCorrectionVector[k]/x; 437 sh += shellCorrectionVector[k]/x;
499 } 438 }
500 sh *= G4Log(tau/taul)/G4Log(taulim/taul); << 439 sh *= std::log(tau/taul)/std::log(taulim/taul);
501 } 440 }
502 sh *= 0.5; 441 sh *= 0.5;
503 return sh; 442 return sh;
504 } 443 }
505 444
506 445
507 //....oooOO0OOooo........oooOO0OOooo........oo 446 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
508 447
509 G4double G4EmCorrections::ShellCorrection(cons 448 G4double G4EmCorrections::ShellCorrection(const G4ParticleDefinition* p,
510 cons << 449 const G4Material* mat,
511 cons << 450 G4double ekin)
512 { 451 {
513 SetupKinematics(p, mat, ekin); 452 SetupKinematics(p, mat, ekin);
>> 453
514 G4double term = 0.0; 454 G4double term = 0.0;
515 //G4cout << "### G4EmCorrections::ShellCorre 455 //G4cout << "### G4EmCorrections::ShellCorrection " << mat->GetName()
516 // << " " << ekin/MeV << " MeV " << << 456 // << " " << ekin/MeV << " MeV " << G4endl;
517 for (G4int i = 0; i<numberOfElements; ++i) { 457 for (G4int i = 0; i<numberOfElements; ++i) {
518 458
519 G4double res = 0.0; 459 G4double res = 0.0;
520 G4double res0 = 0.0; 460 G4double res0 = 0.0;
521 const G4double Z = (*theElementVector)[i]- << 461 G4double Z = (*theElementVector)[i]->GetZ();
522 const G4int iz = (*theElementVector)[i]->G << 462 G4int iz = G4int(Z);
523 G4double Z2= (Z-0.3)*(Z-0.3); 463 G4double Z2= (Z-0.3)*(Z-0.3);
524 G4double f = 1.0; 464 G4double f = 1.0;
525 if(1 == iz) { 465 if(1 == iz) {
526 f = 0.5; 466 f = 0.5;
527 Z2 = 1.0; 467 Z2 = 1.0;
528 } 468 }
529 G4double eta = ba2/Z2; 469 G4double eta = ba2/Z2;
530 G4double tet = (11 < iz) ? sThetaK->Value( << 470 G4double tet = Z2*(1. + Z2*0.25*alpha2);
>> 471 if(11 < iz) { tet = ThetaK->Value(Z); }
531 res0 = f*KShell(tet,eta); 472 res0 = f*KShell(tet,eta);
532 res += res0; 473 res += res0;
533 //G4cout << " Z= " << iz << " Shell 0" << 474 //G4cout << " Z= " << iz << " Shell 0" << " tet= " << tet
534 // << " eta= " << eta << " resK= " 475 // << " eta= " << eta << " resK= " << res0 << G4endl;
535 <<
536 if(2 < iz) { 476 if(2 < iz) {
537 const G4double Zeff = (iz < 10) ? Z - ZD << 477 G4double Zeff = Z - ZD[10];
>> 478 if(iz < 10) { Zeff = Z - ZD[iz]; }
538 Z2= Zeff*Zeff; 479 Z2= Zeff*Zeff;
539 eta = ba2/Z2; 480 eta = ba2/Z2;
540 tet = sThetaL->Value(Z); <<
541 f = 0.125; 481 f = 0.125;
542 const G4int ntot = G4AtomicShells::GetNu << 482 tet = ThetaL->Value(Z);
543 const G4int nmax = std::min(4, ntot); << 483 G4int ntot = G4AtomicShells::GetNumberOfShells(iz);
>> 484 G4int nmax = std::min(4, ntot);
544 G4double norm = 0.0; 485 G4double norm = 0.0;
545 G4double eshell = 0.0; 486 G4double eshell = 0.0;
546 for(G4int j=1; j<nmax; ++j) { 487 for(G4int j=1; j<nmax; ++j) {
547 G4int ne = G4AtomicShells::GetNumberOf 488 G4int ne = G4AtomicShells::GetNumberOfElectrons(iz,j);
548 if(15 >= iz) { << 489 if(15 >= iz) {
549 tet = (3 > j) ? 0.25*Z2*(1.0 + 5*Z2* << 490 if(3 > j) { tet = 0.25*Z2*(1.0 + 5*Z2*alpha2/16.); }
550 0.25*Z2*(1.0 + Z2*alpha2/16.); << 491 else { tet = 0.25*Z2*(1.0 + Z2*alpha2/16.); }
551 } << 492 }
552 norm += ne; << 493 norm += ne;
553 eshell += tet*ne; << 494 eshell += tet*ne;
554 res0 = f*ne*LShell(tet,eta); 495 res0 = f*ne*LShell(tet,eta);
555 res += res0; 496 res += res0;
556 //G4cout << " Zeff= " << Zeff << " She << 497 //G4cout << " Z= " << iz << " Shell " << j << " Ne= " << ne
557 // << " tet= " << tet << " eta= << 498 // << " tet= " << tet << " eta= " << eta
558 // << " resL= " << res0 << G4en << 499 // << " resL= " << res0 << G4endl;
559 } 500 }
560 if(ntot > nmax) { 501 if(ntot > nmax) {
561 if (norm > 0.0) { norm = 1.0/norm; } << 502 eshell /= norm;
562 eshell *= norm; << 503 // Add M-shell
563 <<
564 static const G4double HM[53] = { <<
565 12.0, 12.0, 12.0, 12.0, 11.9, 11.7, 11.5, <<
566 10.0, 9.51, 8.97, 8.52, 8.03, 7.46, 6.95, <<
567 5.61, 5.39, 5.19, 5.01, 4.86, 4.72, 4.62, <<
568 4.32, 4.26, 4.20, 4.15, 4.1, 4.04, 4.00, <<
569 3.90, 3.89, 3.89, 3.88, 3.88, 3.88, 3.88, <<
570 3.90, 3.92, 3.93 }; <<
571 static const G4double HN[31] = { <<
572 75.5, 61.9, 52.2, 45.1, 39.6, 35.4, 31.9, <<
573 24.5, 23.6, 22.7, 22.0, 21.4, 20.9, 20.5, <<
574 19.5, 19.3, 19.2, 19.1, 18.4, 18.8, 18.7, <<
575 18.2}; <<
576 <<
577 // Add M-shell <<
578 if(28 > iz) { 504 if(28 > iz) {
579 res += f*(iz - 10)*LShell(eshell,HM[ 505 res += f*(iz - 10)*LShell(eshell,HM[iz-11]*eta);
580 } else if(63 > iz) { << 506 } else if(63 > iz) {
581 res += f*18*LShell(eshell,HM[iz-11]* 507 res += f*18*LShell(eshell,HM[iz-11]*eta);
582 } else { << 508 } else {
583 res += f*18*LShell(eshell,HM[52]*eta 509 res += f*18*LShell(eshell,HM[52]*eta);
584 } << 510 }
585 // Add N-shell << 511 // Add N-shell
586 if(32 < iz) { 512 if(32 < iz) {
587 if(60 > iz) { 513 if(60 > iz) {
588 res += f*(iz - 28)*LShell(eshell,H << 514 res += f*(iz - 28)*LShell(eshell,HN[iz-33]*eta);
589 } else if(63 > iz) { << 515 } else if(63 > iz) {
590 res += 4*LShell(eshell,HN[iz-33]*e << 516 res += 4*LShell(eshell,HN[iz-33]*eta);
591 } else { << 517 } else {
592 res += 4*LShell(eshell,HN[30]*eta) << 518 res += 4*LShell(eshell,HN[30]*eta);
>> 519 }
>> 520 // Add O-P-shells
>> 521 if(60 < iz) {
>> 522 res += f*(iz - 60)*LShell(eshell,150*eta);
593 } 523 }
594 // Add O-P-shells << 524 }
595 if(60 < iz) { <<
596 res += f*(iz - 60)*LShell(eshell,1 <<
597 } <<
598 } <<
599 } 525 }
600 } 526 }
601 term += res*atomDensity[i]/Z; 527 term += res*atomDensity[i]/Z;
602 } 528 }
603 529
604 term /= material->GetTotNbOfAtomsPerVolume() 530 term /= material->GetTotNbOfAtomsPerVolume();
605 //G4cout << "##Shell Correction=" << term << << 531 //G4cout << "# Shell Correction= " << term << G4endl;
606 return term; 532 return term;
607 } 533 }
608 534
609 //....oooOO0OOooo........oooOO0OOooo........oo 535 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
610 536
611 G4double G4EmCorrections::DensityCorrection(co 537 G4double G4EmCorrections::DensityCorrection(const G4ParticleDefinition* p,
612 co << 538 const G4Material* mat,
613 co << 539 G4double e)
614 { 540 {
615 SetupKinematics(p, mat, e); 541 SetupKinematics(p, mat, e);
616 542
617 G4double cden = material->GetIonisation()-> 543 G4double cden = material->GetIonisation()->GetCdensity();
618 G4double mden = material->GetIonisation()-> 544 G4double mden = material->GetIonisation()->GetMdensity();
619 G4double aden = material->GetIonisation()-> 545 G4double aden = material->GetIonisation()->GetAdensity();
620 G4double x0den = material->GetIonisation()-> 546 G4double x0den = material->GetIonisation()->GetX0density();
621 G4double x1den = material->GetIonisation()-> 547 G4double x1den = material->GetIonisation()->GetX1density();
622 548
>> 549 G4double twoln10 = 2.0*std::log(10.0);
623 G4double dedx = 0.0; 550 G4double dedx = 0.0;
624 551
625 // density correction 552 // density correction
626 static const G4double twoln10 = 2.0*G4Log(10 << 553 G4double x = std::log(bg2)/twoln10;
627 G4double x = G4Log(bg2)/twoln10; <<
628 if ( x >= x0den ) { 554 if ( x >= x0den ) {
629 dedx = twoln10*x - cden ; 555 dedx = twoln10*x - cden ;
630 if ( x < x1den ) dedx += aden*G4Exp(G4Log( << 556 if ( x < x1den ) dedx += aden*std::pow((x1den-x),mden) ;
631 } 557 }
632 558
633 return dedx; 559 return dedx;
634 } 560 }
635 561
636 //....oooOO0OOooo........oooOO0OOooo........oo 562 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
637 563
638 G4double G4EmCorrections::BarkasCorrection(con 564 G4double G4EmCorrections::BarkasCorrection(const G4ParticleDefinition* p,
639 con << 565 const G4Material* mat,
640 con << 566 G4double e)
641 con <<
642 { 567 {
643 // . Z^3 Barkas effect in the stopping power 568 // . Z^3 Barkas effect in the stopping power of matter for charged particles
644 // J.C Ashley and R.H.Ritchie 569 // J.C Ashley and R.H.Ritchie
645 // Physical review B Vol.5 No.7 1 April 19 570 // Physical review B Vol.5 No.7 1 April 1972 pp. 2393-2397
646 // valid for kineticEnergy > 0.5 MeV 571 // valid for kineticEnergy > 0.5 MeV
647 572
648 if (!isInitialized) { SetupKinematics(p, mat << 573 SetupKinematics(p, mat, e);
649 G4double BarkasTerm = 0.0; 574 G4double BarkasTerm = 0.0;
650 575
651 for (G4int i = 0; i<numberOfElements; ++i) { 576 for (G4int i = 0; i<numberOfElements; ++i) {
652 577
653 const G4int iz = (*theElementVector)[i]->G << 578 G4double Z = (*theElementVector)[i]->GetZ();
>> 579 G4int iz = G4int(Z);
654 if(iz == 47) { 580 if(iz == 47) {
655 BarkasTerm += atomDensity[i]*0.006812*G4 << 581 BarkasTerm += atomDensity[i]*0.006812*std::pow(beta,-0.9);
656 } else if(iz >= 64) { 582 } else if(iz >= 64) {
657 BarkasTerm += atomDensity[i]*0.002833*G4 << 583 BarkasTerm += atomDensity[i]*0.002833*std::pow(beta,-1.2);
658 } else { 584 } else {
659 585
660 const G4double Z = (*theElementVector)[i << 586 G4double X = ba2 / Z;
661 const G4double X = ba2 / Z; <<
662 G4double b = 1.3; 587 G4double b = 1.3;
663 if(1 == iz) { b = (material->GetName() = << 588 if(1 == iz) {
>> 589 if(material->GetName() == "G4_lH2") { b = 0.6; }
>> 590 else { b = 1.8; }
>> 591 }
664 else if(2 == iz) { b = 0.6; } 592 else if(2 == iz) { b = 0.6; }
665 else if(10 >= iz) { b = 1.8; } 593 else if(10 >= iz) { b = 1.8; }
666 else if(17 >= iz) { b = 1.4; } 594 else if(17 >= iz) { b = 1.4; }
667 else if(18 == iz) { b = 1.8; } 595 else if(18 == iz) { b = 1.8; }
668 else if(25 >= iz) { b = 1.4; } 596 else if(25 >= iz) { b = 1.4; }
669 else if(50 >= iz) { b = 1.35;} 597 else if(50 >= iz) { b = 1.35;}
670 598
671 const G4double W = b/std::sqrt(X); << 599 G4double W = b/std::sqrt(X);
672 600
673 G4double val = sBarkasCorr->Value(W, idx << 601 G4double val = BarkasCorr->Value(W);
674 if (W > sWmaxBarkas) { val *= (sWmaxBark << 602 if(W > BarkasCorr->Energy(46)) {
>> 603 val *= BarkasCorr->Energy(46)/W;
>> 604 }
675 // G4cout << "i= " << i << " b= " << 605 // G4cout << "i= " << i << " b= " << b << " W= " << W
676 // << " Z= " << Z << " X= " << X << " va 606 // << " Z= " << Z << " X= " << X << " val= " << val<< G4endl;
677 BarkasTerm += val*atomDensity[i] / (std: 607 BarkasTerm += val*atomDensity[i] / (std::sqrt(Z*X)*X);
678 } 608 }
679 } 609 }
680 610
681 BarkasTerm *= 1.29*charge/material->GetTotNb 611 BarkasTerm *= 1.29*charge/material->GetTotNbOfAtomsPerVolume();
682 612
>> 613 // temporary protection
>> 614 //if(charge < -7.0 ) { BarkasTerm *= (-7.0/charge); }
>> 615
683 return BarkasTerm; 616 return BarkasTerm;
684 } 617 }
685 618
686 //....oooOO0OOooo........oooOO0OOooo........oo 619 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
687 620
688 G4double G4EmCorrections::BlochCorrection(cons 621 G4double G4EmCorrections::BlochCorrection(const G4ParticleDefinition* p,
689 cons << 622 const G4Material* mat,
690 cons << 623 G4double e)
691 cons <<
692 { 624 {
693 if (!isInitialized) { SetupKinematics(p, mat << 625 SetupKinematics(p, mat, e);
694 626
695 G4double y2 = q2/ba2; 627 G4double y2 = q2/ba2;
696 628
697 G4double term = 1.0/(1.0 + y2); 629 G4double term = 1.0/(1.0 + y2);
698 G4double del; 630 G4double del;
699 G4double j = 1.0; 631 G4double j = 1.0;
700 do { 632 do {
701 j += 1.0; 633 j += 1.0;
702 del = 1.0/(j* (j*j + y2)); 634 del = 1.0/(j* (j*j + y2));
703 term += del; 635 term += del;
704 // Loop checking, 03-Aug-2015, Vladimir Iv <<
705 } while (del > 0.01*term); 636 } while (del > 0.01*term);
706 637
707 return -y2*term; << 638 G4double res = -y2*term;
>> 639 // temporary protection
>> 640 //if(q2 > 49. && res < -0.2) { res = -0.2; }
>> 641
>> 642 return res;
708 } 643 }
709 644
710 //....oooOO0OOooo........oooOO0OOooo........oo 645 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
711 646
712 G4double G4EmCorrections::MottCorrection(const 647 G4double G4EmCorrections::MottCorrection(const G4ParticleDefinition* p,
713 const << 648 const G4Material* mat,
714 const << 649 G4double e)
715 const <<
716 { 650 {
717 if (!isInitialized) { SetupKinematics(p, mat << 651 SetupKinematics(p, mat, e);
718 return CLHEP::pi*CLHEP::fine_structure_const << 652 G4double mterm = CLHEP::pi*fine_structure_const*beta*charge;
>> 653 return mterm;
719 } 654 }
720 655
721 //....oooOO0OOooo........oooOO0OOooo........oo 656 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
722 657
723 G4double << 658 G4double G4EmCorrections::NuclearDEDX(const G4ParticleDefinition* p,
724 G4EmCorrections::EffectiveChargeCorrection(con << 659 const G4Material* mat,
725 con << 660 G4double e,
726 con << 661 G4bool fluct)
727 { 662 {
728 G4double factor = 1.0; << 663 G4double nloss = 0.0;
729 if(p->GetPDGCharge() <= 2.5*CLHEP::eplus || << 664 if(e <= 0.0) return nloss;
>> 665 SetupKinematics(p, mat, e);
>> 666
>> 667 lossFlucFlag = fluct;
>> 668
>> 669 // Projectile nucleus
>> 670 G4double z1 = std::fabs(particle->GetPDGCharge()/eplus);
>> 671 G4double m1 = mass/amu_c2;
>> 672
>> 673 // loop for the elements in the material
>> 674 for (G4int iel=0; iel<numberOfElements; iel++) {
>> 675 const G4Element* element = (*theElementVector)[iel] ;
>> 676 G4double z2 = element->GetZ();
>> 677 G4double m2 = element->GetA()*mole/g ;
>> 678 nloss += (NuclearStoppingPower(kinEnergy, z1, z2, m1, m2))
>> 679 * atomDensity[iel] ;
>> 680 }
>> 681 nloss *= theZieglerFactor;
>> 682 return nloss;
>> 683 }
>> 684
>> 685 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 686
>> 687 G4double G4EmCorrections::NuclearStoppingPower(G4double kineticEnergy,
>> 688 G4double z1, G4double z2,
>> 689 G4double m1, G4double m2)
>> 690 {
>> 691 G4double energy = kineticEnergy/keV ; // energy in keV
>> 692 G4double nloss = 0.0;
730 693
>> 694 G4double rm;
>> 695 if(z1 > 1.5) rm = (m1 + m2) * ( Z23[G4int(z1)] + Z23[G4int(z2)] ) ;
>> 696 else rm = (m1 + m2) * nist->GetZ13(G4int(z2));
>> 697
>> 698 G4double er = 32.536 * m2 * energy / ( z1 * z2 * rm ) ; // reduced energy
>> 699
>> 700 if (er >= ed[0]) nloss = a[0];
>> 701 else {
>> 702 // the table is inverse in energy
>> 703 for (G4int i=102; i>=0; i--)
>> 704 {
>> 705 if (er <= ed[i]) {
>> 706 nloss = (a[i] - a[i+1])*(er - ed[i+1])/(ed[i] - ed[i+1]) + a[i+1];
>> 707 break;
>> 708 }
>> 709 }
>> 710 }
>> 711
>> 712 // Stragling
>> 713 if(lossFlucFlag) {
>> 714 // G4double sig = 4.0 * m1 * m2 / ((m1 + m2)*(m1 + m2)*
>> 715 // (4.0 + 0.197*std::pow(er,-1.6991)+6.584*std::pow(er,-1.0494))) ;
>> 716 G4double sig = 4.0 * m1 * m2 / ((m1 + m2)*(m1 + m2)*
>> 717 (4.0 + 0.197/(er*er) + 6.584/er));
>> 718
>> 719 nloss *= G4RandGauss::shoot(1.0,sig) ;
>> 720 }
>> 721
>> 722 nloss *= 8.462 * z1 * z2 * m1 / rm ; // Return to [ev/(10^15 atoms/cm^2]
>> 723
>> 724 if ( nloss < 0.0) nloss = 0.0 ;
>> 725
>> 726 return nloss;
>> 727 }
>> 728
>> 729 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 730
>> 731 G4double G4EmCorrections::EffectiveChargeCorrection(const G4ParticleDefinition* p,
>> 732 const G4Material* mat,
>> 733 G4double ekin)
>> 734 {
>> 735 G4double factor = 1.0;
>> 736 if(p->GetPDGCharge() <= 2.5*eplus || nIons <= 0) return factor;
>> 737 /*
731 if(verbose > 1) { 738 if(verbose > 1) {
732 G4cout << "EffectiveChargeCorrection: " << 739 G4cout << "EffectiveChargeCorrection: " << p->GetParticleName()
733 << " in " << mat->GetName() << 740 << " in " << mat->GetName()
734 << " ekin(MeV)= " << ekin/MeV << G4 741 << " ekin(MeV)= " << ekin/MeV << G4endl;
735 } 742 }
736 << 743 */
737 if(p != curParticle || mat != curMaterial) { 744 if(p != curParticle || mat != curMaterial) {
738 curParticle = p; 745 curParticle = p;
739 curMaterial = mat; 746 curMaterial = mat;
740 curVector = nullptr; << 747 curVector = 0;
741 currentZ = p->GetAtomicNumber(); 748 currentZ = p->GetAtomicNumber();
742 if(verbose > 1) { 749 if(verbose > 1) {
743 G4cout << "G4EmCorrections::EffectiveCha 750 G4cout << "G4EmCorrections::EffectiveChargeCorrection: Zion= "
744 << currentZ << " Aion= " << p->Ge << 751 << currentZ << " Aion= " << p->GetPDGMass()/amu_c2 << G4endl;
745 } 752 }
746 massFactor = CLHEP::proton_mass_c2/p->GetP << 753 massFactor = proton_mass_c2/p->GetPDGMass();
747 idx = -1; 754 idx = -1;
748 755
749 for(G4int i=0; i<nIons; ++i) { 756 for(G4int i=0; i<nIons; ++i) {
750 if(materialList[i] == mat && currentZ == 757 if(materialList[i] == mat && currentZ == Zion[i]) {
751 idx = i; 758 idx = i;
752 break; 759 break;
753 } 760 }
754 } 761 }
755 //G4cout << " idx= " << idx << " dz= " << << 762 // G4cout << " idx= " << idx << " dz= " << G4endl;
756 if(idx >= 0) { 763 if(idx >= 0) {
757 if(nullptr == ionList[idx]) { BuildCorre << 764 if(!ionList[idx]) BuildCorrectionVector();
758 curVector = stopData[idx]; << 765 if(ionList[idx]) curVector = stopData[idx];
759 } else { << 766 } else { return factor; }
760 return factor; <<
761 } <<
762 } 767 }
763 if(nullptr != curVector) { << 768 if(curVector) {
764 factor = curVector->Value(ekin*massFactor) 769 factor = curVector->Value(ekin*massFactor);
765 if(verbose > 1) { 770 if(verbose > 1) {
766 G4cout << "E= " << ekin << " factor= " < 771 G4cout << "E= " << ekin << " factor= " << factor << " massfactor= "
767 << massFactor << G4endl; << 772 << massFactor << G4endl;
768 } 773 }
769 } 774 }
770 return factor; 775 return factor;
771 } 776 }
772 777
773 //....oooOO0OOooo........oooOO0OOooo........oo 778 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
774 779
775 void G4EmCorrections::AddStoppingData(const G4 << 780 void G4EmCorrections::AddStoppingData(G4int Z, G4int A,
776 const G4 << 781 const G4String& mname,
777 G4Physic << 782 G4PhysicsVector* dVector)
778 { 783 {
779 G4int i = 0; 784 G4int i = 0;
780 for(; i<nIons; ++i) { 785 for(; i<nIons; ++i) {
781 if(Z == Zion[i] && A == Aion[i] && mname = 786 if(Z == Zion[i] && A == Aion[i] && mname == materialName[i]) break;
782 } 787 }
783 if(i == nIons) { 788 if(i == nIons) {
784 Zion.push_back(Z); 789 Zion.push_back(Z);
785 Aion.push_back(A); 790 Aion.push_back(A);
786 materialName.push_back(mname); 791 materialName.push_back(mname);
787 materialList.push_back(nullptr); << 792 materialList.push_back(0);
788 ionList.push_back(nullptr); << 793 ionList.push_back(0);
789 stopData.push_back(dVector); 794 stopData.push_back(dVector);
790 nIons++; 795 nIons++;
791 if(verbose > 1) { << 796 if(verbose>1) {
792 G4cout << "AddStoppingData Z= " << Z << 797 G4cout << "AddStoppingData Z= " << Z << " A= " << A << " " << mname
793 << " idx= " << i << G4endl; << 798 << " idx= " << i << G4endl;
794 } 799 }
795 } 800 }
796 } 801 }
797 802
798 //....oooOO0OOooo........oooOO0OOooo........oo 803 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
799 804
800 void G4EmCorrections::BuildCorrectionVector() 805 void G4EmCorrections::BuildCorrectionVector()
801 { 806 {
802 if(nullptr == ionLEModel || nullptr == ionHE << 807 if(!ionLEModel || !ionHEModel) {
803 return; 808 return;
804 } 809 }
805 810
806 const G4ParticleDefinition* ion = curParticl 811 const G4ParticleDefinition* ion = curParticle;
807 const G4ParticleDefinition* gion = G4Generic <<
808 G4int Z = Zion[idx]; 812 G4int Z = Zion[idx];
809 G4double A = Aion[idx]; << 813 if(currentZ != Z) {
>> 814 ion = G4ParticleTable::GetParticleTable()->FindIon(Z, Aion[idx], 0, Z);
>> 815 }
>> 816 //G4cout << "BuildCorrectionVector: idx= " << idx << " Z= " << Z
>> 817 // << " curZ= " << currentZ << G4endl;
>> 818
>> 819 // G4double A = nist->GetAtomicMassAmu(Z);
>> 820 G4double A = G4double(ion->GetBaryonNumber());
810 G4PhysicsVector* v = stopData[idx]; 821 G4PhysicsVector* v = stopData[idx];
811 << 822
812 if(verbose > 1) { << 823 const G4ParticleDefinition* p = G4GenericIon::GenericIon();
>> 824 G4double massRatio = proton_mass_c2/ion->GetPDGMass();
>> 825
>> 826 if(verbose>1) {
813 G4cout << "BuildCorrectionVector: Stopping 827 G4cout << "BuildCorrectionVector: Stopping for "
814 << curParticle->GetParticleName() < << 828 << curParticle->GetParticleName() << " in "
815 << materialName[idx] << " Ion Z= " << 829 << materialName[idx] << " Ion Z= " << Z << " A= " << A
816 << " massFactor= " << massFactor << << 830 << " massRatio= " << massRatio << G4endl;
817 G4cout << " Nbins=" << nbinCorr << " Em <<
818 << " Emax(MeV)=" << eCorrMax << " ion: " <<
819 << ion->GetParticleName() << G4endl <<
820 } 831 }
821 832
822 auto vv = new G4PhysicsLogVector(eCorrMin,eC << 833 G4PhysicsLogVector* vv =
823 const G4double eth0 = v->Energy(0); << 834 new G4PhysicsLogVector(eCorrMin,eCorrMax,nbinCorr);
824 const G4double escal = eth/massFactor; << 835 vv->SetSpline(true);
>> 836 G4double e, eion, dedx, dedx1;
>> 837 G4double eth0 = v->Energy(0);
>> 838 G4double escal = eth/massRatio;
825 G4double qe = 839 G4double qe =
826 effCharge.EffectiveChargeSquareRatio(curPa << 840 effCharge.EffectiveChargeSquareRatio(ion, curMaterial, escal);
827 const G4double dedxt = << 841 G4double dedxt =
828 ionLEModel->ComputeDEDXPerVolume(curMateri << 842 ionLEModel->ComputeDEDXPerVolume(curMaterial, p, eth, eth)*qe;
829 const G4double dedx1t = << 843 G4double dedx1t =
830 ionHEModel->ComputeDEDXPerVolume(curMateri << 844 ionHEModel->ComputeDEDXPerVolume(curMaterial, p, eth, eth)*qe
831 + ComputeIonCorrections(curParticle, curMa 845 + ComputeIonCorrections(curParticle, curMaterial, escal);
832 const G4double rest = escal*(dedxt - dedx1t) << 846 G4double rest = escal*(dedxt - dedx1t);
833 if(verbose > 1) { << 847 //G4cout << "Escal(MeV)= "<<escal<<" dedxt0= " <<dedxt
834 G4cout << "Escal(MeV)= " << escal << " qe= << 848 // << " dedxt1= " << dedx1t << G4endl;
835 << " dedxt= " << dedxt << " dedx1t= << 849
836 } <<
837 for(G4int i=0; i<=nbinCorr; ++i) { 850 for(G4int i=0; i<=nbinCorr; ++i) {
838 // energy in the local table (GenericIon) << 851 e = vv->Energy(i);
839 G4double e = vv->Energy(i); << 852 escal = e/massRatio;
840 // energy of the real ion << 853 eion = escal/A;
841 G4double eion = e/massFactor; << 854 if(eion <= eth0) {
842 // energy in the imput stopping data vecto << 855 dedx = v->Value(eth0)*std::sqrt(eion/eth0);
843 G4double e1 = eion/A; << 856 } else {
844 G4double dedx = (e1 < eth0) << 857 dedx = v->Value(eion);
845 ? v->Value(eth0)*std::sqrt(e1/eth0) : v- << 858 }
846 qe = effCharge.EffectiveChargeSquareRatio( << 859 qe = effCharge.EffectiveChargeSquareRatio(curParticle,curMaterial,escal);
847 G4double dedx1 = (e <= eth) << 860 if(e <= eth) {
848 ? ionLEModel->ComputeDEDXPerVolume(curMa << 861 dedx1 = ionLEModel->ComputeDEDXPerVolume(curMaterial, p, e, e)*qe;
849 : ionHEModel->ComputeDEDXPerVolume(curMa << 862 } else {
850 ComputeIonCorrections(curParticle, curMa << 863 dedx1 = ionHEModel->ComputeDEDXPerVolume(curMaterial, p, e, e)*qe +
>> 864 ComputeIonCorrections(curParticle, curMaterial, escal) + rest/escal;
>> 865 }
851 vv->PutValue(i, dedx/dedx1); 866 vv->PutValue(i, dedx/dedx1);
852 if(verbose > 1) { << 867 if(verbose>1) {
853 G4cout << "E(MeV)=" << e/CLHEP::MeV << " << 868 G4cout << " E(meV)= " << e/MeV << " Correction= " << dedx/dedx1
854 << " e1=" << e1 << " dedxRatio= " << de << 869 << " " << dedx << " " << dedx1
855 << " dedx=" << dedx << " dedx1=" << 870 << " massF= " << massFactor << G4endl;
856 << " qe=" << qe << " rest/eion=" <<
857 } 871 }
858 } 872 }
859 delete v; 873 delete v;
860 ionList[idx] = ion; 874 ionList[idx] = ion;
861 stopData[idx] = vv; 875 stopData[idx] = vv;
862 if(verbose > 1) { G4cout << "End data set " << 876 if(verbose>1) { G4cout << "End data set " << G4endl; }
863 } 877 }
864 878
865 //....oooOO0OOooo........oooOO0OOooo........oo 879 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
866 880
867 void G4EmCorrections::InitialiseForNewRun() 881 void G4EmCorrections::InitialiseForNewRun()
868 { 882 {
869 G4ProductionCutsTable* tb = G4ProductionCuts 883 G4ProductionCutsTable* tb = G4ProductionCutsTable::GetProductionCutsTable();
870 ncouples = tb->GetTableSize(); 884 ncouples = tb->GetTableSize();
871 if(currmat.size() != ncouples) { 885 if(currmat.size() != ncouples) {
872 currmat.resize(ncouples); 886 currmat.resize(ncouples);
873 for(auto it = thcorr.begin(); it != thcorr << 887 size_t i;
874 (it->second).clear(); << 888 for(i=0; i<100; ++i) {thcorr[i].clear();}
875 } << 889 for(i=0; i<ncouples; ++i) {
876 thcorr.clear(); << 890 currmat[i] = tb->GetMaterialCutsCouple(i)->GetMaterial();
877 for(std::size_t i=0; i<ncouples; ++i) { <<
878 currmat[i] = tb->GetMaterialCutsCouple(( <<
879 G4String nam = currmat[i]->GetName(); 891 G4String nam = currmat[i]->GetName();
880 for(G4int j=0; j<nIons; ++j) { 892 for(G4int j=0; j<nIons; ++j) {
881 if(nam == materialName[j]) { materialL << 893 if(nam == materialName[j]) { materialList[j] = currmat[i]; }
882 } 894 }
883 } 895 }
884 } 896 }
885 } 897 }
886 898
887 //....oooOO0OOooo........oooOO0OOooo........oo 899 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
888 900
889 void G4EmCorrections::Initialise() 901 void G4EmCorrections::Initialise()
890 { 902 {
891 // Z^3 Barkas effect in the stopping power o << 903 // . Z^3 Barkas effect in the stopping power of matter for charged particles
892 // J.C Ashley and R.H.Ritchie << 904 // J.C Ashley and R.H.Ritchie
893 // Physical review B Vol.5 No.7 1 April 1972 << 905 // 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; 906 G4int i, j;
898 static const G4double fTable[47][2] = { << 907 const G4double fTable[47][2] = {
899 { 0.02, 21.5}, 908 { 0.02, 21.5},
900 { 0.03, 20.0}, 909 { 0.03, 20.0},
901 { 0.04, 18.0}, 910 { 0.04, 18.0},
902 { 0.05, 15.6}, 911 { 0.05, 15.6},
903 { 0.06, 15.0}, 912 { 0.06, 15.0},
904 { 0.07, 14.0}, 913 { 0.07, 14.0},
905 { 0.08, 13.5}, 914 { 0.08, 13.5},
906 { 0.09, 13.}, 915 { 0.09, 13.},
907 { 0.1, 12.2}, 916 { 0.1, 12.2},
908 { 0.2, 9.25}, 917 { 0.2, 9.25},
909 { 0.3, 7.0}, 918 { 0.3, 7.0},
910 { 0.4, 6.0}, 919 { 0.4, 6.0},
911 { 0.5, 4.5}, 920 { 0.5, 4.5},
912 { 0.6, 3.5}, 921 { 0.6, 3.5},
913 { 0.7, 3.0}, 922 { 0.7, 3.0},
914 { 0.8, 2.5}, 923 { 0.8, 2.5},
915 { 0.9, 2.0}, 924 { 0.9, 2.0},
916 { 1.0, 1.7}, 925 { 1.0, 1.7},
917 { 1.2, 1.2}, 926 { 1.2, 1.2},
918 { 1.3, 1.0}, 927 { 1.3, 1.0},
919 { 1.4, 0.86}, 928 { 1.4, 0.86},
920 { 1.5, 0.7}, 929 { 1.5, 0.7},
921 { 1.6, 0.61}, 930 { 1.6, 0.61},
922 { 1.7, 0.52}, 931 { 1.7, 0.52},
923 { 1.8, 0.5}, 932 { 1.8, 0.5},
924 { 1.9, 0.43}, 933 { 1.9, 0.43},
925 { 2.0, 0.42}, 934 { 2.0, 0.42},
926 { 2.1, 0.3}, 935 { 2.1, 0.3},
927 { 2.4, 0.2}, 936 { 2.4, 0.2},
928 { 3.0, 0.13}, 937 { 3.0, 0.13},
929 { 3.08, 0.1}, 938 { 3.08, 0.1},
930 { 3.1, 0.09}, 939 { 3.1, 0.09},
931 { 3.3, 0.08}, 940 { 3.3, 0.08},
932 { 3.5, 0.07}, 941 { 3.5, 0.07},
933 { 3.8, 0.06}, 942 { 3.8, 0.06},
934 { 4.0, 0.051}, 943 { 4.0, 0.051},
935 { 4.1, 0.04}, 944 { 4.1, 0.04},
936 { 4.8, 0.03}, 945 { 4.8, 0.03},
937 { 5.0, 0.024}, 946 { 5.0, 0.024},
938 { 5.1, 0.02}, 947 { 5.1, 0.02},
939 { 6.0, 0.013}, 948 { 6.0, 0.013},
940 { 6.5, 0.01}, 949 { 6.5, 0.01},
941 { 7.0, 0.009}, 950 { 7.0, 0.009},
942 { 7.1, 0.008}, 951 { 7.1, 0.008},
943 { 8.0, 0.006}, 952 { 8.0, 0.006},
944 { 9.0, 0.0032}, 953 { 9.0, 0.0032},
945 { 10.0, 0.0025} }; 954 { 10.0, 0.0025} };
946 955
947 sBarkasCorr = new G4PhysicsFreeVector(47, fa << 956 BarkasCorr = new G4LPhysicsFreeVector(47, 0.02, 10.);
948 for(i=0; i<47; ++i) { sBarkasCorr->PutValues << 957 for(i=0; i<47; ++i) { BarkasCorr->PutValues(i, fTable[i][0], fTable[i][1]); }
>> 958 BarkasCorr->SetSpline(true);
>> 959
>> 960 const G4double nuca[104][2] = {
>> 961 { 1.0E+8, 5.831E-8},
>> 962 { 8.0E+7, 7.288E-8},
>> 963 { 6.0E+7, 9.719E-8},
>> 964 { 5.0E+7, 1.166E-7},
>> 965 { 4.0E+7, 1.457E-7},
>> 966 { 3.0E+7, 1.942E-7},
>> 967 { 2.0E+7, 2.916E-7},
>> 968 { 1.5E+7, 3.887E-7},
>> 969
>> 970 { 1.0E+7, 5.833E-7},
>> 971 { 8.0E+6, 7.287E-7},
>> 972 { 6.0E+6, 9.712E-7},
>> 973 { 5.0E+6, 1.166E-6},
>> 974 { 4.0E+6, 1.457E-6},
>> 975 { 3.0E+6, 1.941E-6},
>> 976 { 2.0E+6, 2.911E-6},
>> 977 { 1.5E+6, 3.878E-6},
>> 978
>> 979 { 1.0E+6, 5.810E-6},
>> 980 { 8.0E+5, 7.262E-6},
>> 981 { 6.0E+5, 9.663E-6},
>> 982 { 5.0E+5, 1.157E-5},
>> 983 { 4.0E+5, 1.442E-5},
>> 984 { 3.0E+5, 1.913E-5},
>> 985 { 2.0E+5, 2.845E-5},
>> 986 { 1.5E+5, 3.762E-5},
>> 987
>> 988 { 1.0E+5, 5.554E-5},
>> 989 { 8.0E+4, 6.866E-5},
>> 990 { 6.0E+4, 9.020E-5},
>> 991 { 5.0E+4, 1.070E-4},
>> 992 { 4.0E+4, 1.319E-4},
>> 993 { 3.0E+4, 1.722E-4},
>> 994 { 2.0E+4, 2.499E-4},
>> 995 { 1.5E+4, 3.248E-4},
>> 996
>> 997 { 1.0E+4, 4.688E-4},
>> 998 { 8.0E+3, 5.729E-4},
>> 999 { 6.0E+3, 7.411E-4},
>> 1000 { 5.0E+3, 8.718E-4},
>> 1001 { 4.0E+3, 1.063E-3},
>> 1002 { 3.0E+3, 1.370E-3},
>> 1003 { 2.0E+3, 1.955E-3},
>> 1004 { 1.5E+3, 2.511E-3},
>> 1005
>> 1006 { 1.0E+3, 3.563E-3},
>> 1007 { 8.0E+2, 4.314E-3},
>> 1008 { 6.0E+2, 5.511E-3},
>> 1009 { 5.0E+2, 6.430E-3},
>> 1010 { 4.0E+2, 7.756E-3},
>> 1011 { 3.0E+2, 9.855E-3},
>> 1012 { 2.0E+2, 1.375E-2},
>> 1013 { 1.5E+2, 1.736E-2},
>> 1014
>> 1015 { 1.0E+2, 2.395E-2},
>> 1016 { 8.0E+1, 2.850E-2},
>> 1017 { 6.0E+1, 3.552E-2},
>> 1018 { 5.0E+1, 4.073E-2},
>> 1019 { 4.0E+1, 4.802E-2},
>> 1020 { 3.0E+1, 5.904E-2},
>> 1021 { 1.5E+1, 9.426E-2},
>> 1022
>> 1023 { 1.0E+1, 1.210E-1},
>> 1024 { 8.0E+0, 1.377E-1},
>> 1025 { 6.0E+0, 1.611E-1},
>> 1026 { 5.0E+0, 1.768E-1},
>> 1027 { 4.0E+0, 1.968E-1},
>> 1028 { 3.0E+0, 2.235E-1},
>> 1029 { 2.0E+0, 2.613E-1},
>> 1030 { 1.5E+0, 2.871E-1},
>> 1031
>> 1032 { 1.0E+0, 3.199E-1},
>> 1033 { 8.0E-1, 3.354E-1},
>> 1034 { 6.0E-1, 3.523E-1},
>> 1035 { 5.0E-1, 3.609E-1},
>> 1036 { 4.0E-1, 3.693E-1},
>> 1037 { 3.0E-1, 3.766E-1},
>> 1038 { 2.0E-1, 3.803E-1},
>> 1039 { 1.5E-1, 3.788E-1},
>> 1040
>> 1041 { 1.0E-1, 3.711E-1},
>> 1042 { 8.0E-2, 3.644E-1},
>> 1043 { 6.0E-2, 3.530E-1},
>> 1044 { 5.0E-2, 3.444E-1},
>> 1045 { 4.0E-2, 3.323E-1},
>> 1046 { 3.0E-2, 3.144E-1},
>> 1047 { 2.0E-2, 2.854E-1},
>> 1048 { 1.5E-2, 2.629E-1},
>> 1049
>> 1050 { 1.0E-2, 2.298E-1},
>> 1051 { 8.0E-3, 2.115E-1},
>> 1052 { 6.0E-3, 1.883E-1},
>> 1053 { 5.0E-3, 1.741E-1},
>> 1054 { 4.0E-3, 1.574E-1},
>> 1055 { 3.0E-3, 1.372E-1},
>> 1056 { 2.0E-3, 1.116E-1},
>> 1057 { 1.5E-3, 9.559E-2},
>> 1058
>> 1059 { 1.0E-3, 7.601E-2},
>> 1060 { 8.0E-4, 6.668E-2},
>> 1061 { 6.0E-4, 5.605E-2},
>> 1062 { 5.0E-4, 5.008E-2},
>> 1063 { 4.0E-4, 4.352E-2},
>> 1064 { 3.0E-4, 3.617E-2},
>> 1065 { 2.0E-4, 2.768E-2},
>> 1066 { 1.5E-4, 2.279E-2},
>> 1067
>> 1068 { 1.0E-4, 1.723E-2},
>> 1069 { 8.0E-5, 1.473E-2},
>> 1070 { 6.0E-5, 1.200E-2},
>> 1071 { 5.0E-5, 1.052E-2},
>> 1072 { 4.0E-5, 8.950E-3},
>> 1073 { 3.0E-5, 7.246E-3},
>> 1074 { 2.0E-5, 5.358E-3},
>> 1075 { 1.5E-5, 4.313E-3},
>> 1076 { 0.0, 3.166E-3}
>> 1077 };
>> 1078
>> 1079 for(i=0; i<104; ++i) {
>> 1080 ed[i] = nuca[i][0];
>> 1081 a[i] = nuca[i][1];
>> 1082 }
949 1083
950 const G4double SK[20] = {1.9477, 1.9232, 1.8 << 1084 // Constants
>> 1085 theZieglerFactor = eV*cm2*1.0e-15 ;
>> 1086 alpha2 = fine_structure_const*fine_structure_const;
>> 1087 lossFlucFlag = true;
>> 1088
>> 1089 // G.S. Khandelwal Nucl. Phys. A116(1968)97 - 111.
>> 1090 // "Shell corrections for K- and L- electrons
>> 1091
>> 1092 nK = 20;
>> 1093 nL = 26;
>> 1094 nEtaK = 29;
>> 1095 nEtaL = 28;
>> 1096
>> 1097 const G4double d[11] = {0., 0., 0., 1.72, 2.09, 2.48, 2.82, 3.16, 3.53, 3.84, 4.15};
>> 1098 const G4double thek[20] = {0.64, 0.65, 0.66, 0.68, 0.70, 0.72, 0.74, 0.75, 0.76, 0.78,
>> 1099 0.80, 0.82, 0.84, 0.85, 0.86, 0.88, 0.90, 0.92, 0.94, 0.95};
>> 1100 const G4double sk[20] = {1.9477, 1.9232, 1.8996, 1.8550, 1.8137,
951 1.7754, 1.7396, 1.7 1101 1.7754, 1.7396, 1.7223, 1.7063, 1.6752,
952 1.6461, 1.6189, 1.5 1102 1.6461, 1.6189, 1.5933, 1.5811, 1.5693,
953 1.5467, 1.5254, 1.5 1103 1.5467, 1.5254, 1.5053, 1.4863, 1.4772};
954 const G4double TK[20] = {2.5222, 2.5125, 2.5 << 1104 const G4double tk[20] = {2.5222, 2.5125, 2.5026, 2.4821, 2.4608,
955 2.4388, 2.4163, 2.4 1105 2.4388, 2.4163, 2.4044, 2.3933, 2.3701,
956 2.3466, 2.3229, 2.2 1106 2.3466, 2.3229, 2.2992, 2.2872, 2.2753,
957 2.2515, 2.2277, 2.2 1107 2.2515, 2.2277, 2.2040, 2.1804, 2.1686};
>> 1108 const G4double uk[20] = {1.9999, 2.0134, 2.0258, 2.0478, 2.0662,
>> 1109 2.0817, 2.0945, 2.0999, 2.1049, 2.1132,
>> 1110 2.1197, 2.1246, 2.1280, 2.1292, 2.1301,
>> 1111 2.1310, 2.1310, 2.1300, 2.1283, 2.1271};
>> 1112 const G4double vk[20] = {8.3410, 8.3373, 8.3340, 8.3287, 8.3247,
>> 1113 8.3219, 8.3201, 8.3194, 8.3191, 8.3188,
>> 1114 8.3191, 8.3199, 8.3211, 8.3218, 8.3226,
>> 1115 8.3244, 8.3264, 8.3285, 8.3308, 8.3320};
>> 1116
>> 1117 for(i=0; i<11; ++i) { ZD[i] = d[i];}
>> 1118
>> 1119 for(i=0; i<nK; ++i) {
>> 1120 TheK[i] = thek[i];
>> 1121 SK[i] = sk[i];
>> 1122 TK[i] = tk[i];
>> 1123 UK[i] = uk[i];
>> 1124 VK[i] = vk[i];
>> 1125 }
958 1126
959 const G4double SL[26] = {15.3343, 13.9389, 1 << 1127 const G4double thel[26] = {0.24, 0.26, 0.28, 0.30, 0.32, 0.34, 0.35, 0.36, 0.38, 0.40,
>> 1128 0.42, 0.44, 0.45, 0.46, 0.48, 0.50, 0.52, 0.54, 0.55, 0.56,
>> 1129 0.58, 0.60, 0.62, 0.64, 0.65, 0.66};
>> 1130 const G4double sl[26] = {15.3343, 13.9389, 12.7909, 11.8343, 11.0283,
960 10.3424, 10.0371, 1131 10.3424, 10.0371, 9.7537, 9.2443, 8.8005,
961 8.4114, 8.0683, << 1132 8.4114, 8.0683, 7.9117, 7.7641, 7.4931,
962 7.2506, 7.0327, 1133 7.2506, 7.0327, 6.8362, 6.7452, 6.6584,
963 6.4969, 6.3498, 1134 6.4969, 6.3498, 6.2154, 6.0923, 6.0345, 5.9792};
964 const G4double TL[26] = {35.0669, 33.4344, 3 << 1135 const G4double tl[26] = {35.0669, 33.4344, 32.0073, 30.7466, 29.6226,
965 28.6128, 28.1449, 2 1136 28.6128, 28.1449, 27.6991, 26.8674, 26.1061,
966 25.4058, 24.7587, 2 1137 25.4058, 24.7587, 24.4531, 24.1583, 23.5992,
967 23.0771, 22.5880, 2 1138 23.0771, 22.5880, 22.1285, 21.9090, 21.6958,
968 21.2872, 20.9006, 2 1139 21.2872, 20.9006, 20.5341, 20.1859, 20.0183, 19.8546};
>> 1140 const G4double ul[26] = {0.1215, 0.5265, 0.8411, 1.0878, 1.2828,
>> 1141 1.4379, 1.5032, 1.5617, 1.6608, 1.7401,
>> 1142 1.8036, 1.8543, 1.8756, 1.8945, 1.9262,
>> 1143 1.9508, 1.9696, 1.9836, 1.9890, 1.9935,
>> 1144 2.0001, 2.0039, 2.0053, 2.0049, 2.0040, 2.0028};
>> 1145 for(i=0; i<nL; ++i) {
>> 1146 TheL[i] = thel[i];
>> 1147 SL[i] = sl[i];
>> 1148 TL[i] = tl[i];
>> 1149 UL[i] = ul[i];
>> 1150 }
>> 1151
>> 1152 const G4double eta[29] = {0.005, 0.007, 0.01, 0.015, 0.02,
>> 1153 0.03, 0.04, 0.05, 0.06, 0.08,
>> 1154 0.1, 0.15, 0.2, 0.3, 0.4,
>> 1155 0.5, 0.6, 0.7, 0.8, 1.0,
>> 1156 1.2, 1.4, 1.5, 1.7, 2.0, 3.0, 5.0, 7.0, 10.};
969 1157
970 const G4double bk1[29][11] = { 1158 const G4double bk1[29][11] = {
971 {0.005, 1.34782E-8, 1.46132E-8, 1.72179E-8, 1159 {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, 1160 {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 1161 {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, 1162 {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 1163 {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 1164 {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 1165 {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 1166 {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 1167 {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 1168 {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. 1169 {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 1170 {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. 1171 {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. 1172 {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. 1173 {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. 1174 {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. 1175 {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. 1176 {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. 1177 {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. 1178 {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. 1179 {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. 1180 {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. 1181 {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. 1182 {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. 1183 {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. 1184 {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. 1185 {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. 1186 {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 1187 {10.0, 5.32590, 5.35848, 5.42560, 5.49547, 5.56830, 5.64429, 5.68353, 5.72366, 5.80666, 5.89359}
1000 }; 1188 };
1001 1189
1002 const G4double bk2[29][11] = { 1190 const G4double bk2[29][11] = {
1003 {0.005, 5.98040E-8, 7.25636E-8, 8.00602E-8, 1191 {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, 1192 {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, 1193 {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, 1194 {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, 1195 {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, 1196 {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, 1197 {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, 1198 {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, 1199 {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, 1200 {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 1201 {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, 1202 {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 1203 {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 1204 {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 1205 {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 1206 {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 1207 {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 1208 {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 1209 {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 1210 {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 1211 {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 1212 {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 1213 {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 1214 {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 1215 {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 1216 {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 1217 {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 1218 {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, 1219 {10.0, 5.98474, 6.08046, 6.13015, 6.18112, 6.28715, 6.39903, 6.51728, 6.64249, 6.70792, 6.77535}
1032 }; 1220 };
1033 1221
1034 const G4double bls1[28][10] = { 1222 const G4double bls1[28][10] = {
1035 {0.005, 2.4111E-4, 2.5612E-4, 2.7202E-4, 3. 1223 {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. 1224 {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 1225 {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. 1226 {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 1227 {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 1228 {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 1229 {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 1230 {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 1231 {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 1232 {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 1233 {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 1234 {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 1235 {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 1236 {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 1237 {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 1238 {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 1239 {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 1240 {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 1241 {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 1242 {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 1243 {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 1244 {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 1245 {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 1246 {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 1247 {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 1248 {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 1249 {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 1250 {7.0, 6.8634, 6.9194, 6.9767, 7.0957, 7.2208, 7.3526, 7.4915, 7.5639, 7.6384}
1063 }; 1251 };
1064 1252
>> 1253
1065 const G4double bls2[28][10] = { 1254 const G4double bls2[28][10] = {
1066 {0.005, 5.4561E-4, 6.0905E-4, 6.7863E-4, 7. 1255 {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. 1256 {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 1257 {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. 1258 {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 1259 {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 1260 {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 1261 {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 1262 {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 1263 {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 1264 {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 1265 {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 1266 {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 1267 {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 1268 {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 1269 {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 1270 {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 1271 {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 1272 {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 1273 {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 1274 {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 1275 {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 1276 {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 1277 {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 1278 {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 1279 {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 1280 {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 1281 {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 1282 {7.0, 7.7938, 7.9588, 8.1342, 8.3211, 8.4193, 8.5209, 8.7350, 8.9651, 9.2133}
1094 }; 1283 };
1095 1284
1096 const G4double bls3[28][9] = { 1285 const G4double bls3[28][9] = {
1097 {0.005, 1.2895E-3, 1.3670E-3, 1.4524E-3, 1. 1286 {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. 1287 {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 1288 {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 1289 {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 1290 {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 1291 {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 1292 {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 1293 {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 1294 {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 1295 {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 1296 {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 1297 {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 1298 {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 1299 {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 1300 {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 1301 {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 1302 {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 1303 {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 1304 {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 1305 {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 1306 {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 1307 {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 1308 {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 1309 {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 1310 {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 1311 {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 1312 {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 1313 {7.0, 9.4819, 9.6248, 9.7739, 10.0926, 10.4423, 10.8282, 11.2565, 11.7356}
1125 }; 1314 };
1126 1315
1127 const G4double bll1[28][10] = { 1316 const G4double bll1[28][10] = {
1128 {0.005, 3.6324E-5, 4.0609E-5, 4.5430E-5, 5. 1317 {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. 1318 {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 1319 {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. 1320 {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 1321 {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 1322 {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 1323 {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 1324 {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 1325 {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 1326 {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 1327 {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 1328 {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 1329 {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 1330 {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 1331 {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 1332 {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 1333 {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 1334 {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 1335 {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 1336 {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 1337 {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 1338 {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 1339 {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 1340 {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 1341 {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 1342 {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 1343 {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 1344 {7.0, 7.7362, 7.8079, 7.8821, 8.0383, 8.2061, 8.3866, 8.5816, 8.6850, 8.7927}
1156 }; 1345 };
1157 1346
1158 const G4double bll2[28][10] = { 1347 const G4double bll2[28][10] = {
1159 {0.005, 1.8339E-4, 2.3330E-4, 2.9738E-4, 3. 1348 {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. 1349 {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 1350 {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. 1351 {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 1352 {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 1353 {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 1354 {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 1355 {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 1356 {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 1357 {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 1358 {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 1359 {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 1360 {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 1361 {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 1362 {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 1363 {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 1364 {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 1365 {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 1366 {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 1367 {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 1368 {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 1369 {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 1370 {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 1371 {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 1372 {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 1373 {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 1374 {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 1375 {7.0, 9.0221, 9.2724, 9.5464, 9.8477, 10.0099, 10.1805, 10.5499, 10.9622, 11.4250}
1187 }; 1376 };
1188 1377
1189 const G4double bll3[28][9] = { 1378 const G4double bll3[28][9] = {
1190 {0.005, 1.3190E-3, 1.4961E-3, 1.6974E-3, 2. 1379 {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. 1380 {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 1381 {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. 1382 {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 1383 {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 1384 {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 1385 {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 1386 {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 1387 {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 1388 {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 1389 {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 1390 {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 1391 {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 1392 {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 1393 {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 1394 {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 1395 {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 1396 {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 1397 {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 1398 {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 1399 {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 1400 {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 1401 {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 1402 {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 1403 {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 1404 {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 1405 {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 1406 {7.0, 11.9480, 12.2357, 12.5432, 13.2260, 14.0164, 14.9404, 16.0330, 17.3420}
1218 }; 1407 };
1219 1408
1220 G4double b, bs; 1409 G4double b, bs;
1221 for(i=0; i<nEtaK; ++i) { 1410 for(i=0; i<nEtaK; ++i) {
1222 1411
1223 G4double et = Eta[i]; << 1412 G4double et = eta[i];
1224 G4double loget = G4Log(et); << 1413 G4double loget = std::log(et);
>> 1414 Eta[i] = et;
>> 1415 //G4cout << "### eta["<<i<<"]="<<et<<" KShell: tet= "<<TheK[0]<<" - "<<TheK[nK-1]<<G4endl;
1225 1416
1226 for(j=0; j<nK; ++j) { 1417 for(j=0; j<nK; ++j) {
1227 1418
1228 if(j < 10) { b = bk2[i][10-j]; } 1419 if(j < 10) { b = bk2[i][10-j]; }
1229 else { b = bk1[i][20-j]; } 1420 else { b = bk1[i][20-j]; }
1230 1421
1231 CK[j][i] = SK[j]*loget + TK[j] - b; 1422 CK[j][i] = SK[j]*loget + TK[j] - b;
1232 1423
1233 if(i == nEtaK-1) { 1424 if(i == nEtaK-1) {
1234 ZK[j] = et*(et*et*CK[j][i] - et*UK[j] << 1425 ZK[j] = et*(et*et*CK[j][i] - et*UK[j] - VK[j]);
1235 //G4cout << "i= " << i << " j= " << j << 1426 //G4cout << "i= " << i << " j= " << j
1236 // << " CK[j][i]= " << CK[j][i << 1427 // << " CK[j][i]= " << CK[j][i]
1237 // << " ZK[j]= " << ZK[j] << " << 1428 // << " ZK[j]= " << ZK[j] << " b= " << b << G4endl;
1238 } 1429 }
1239 } 1430 }
1240 //G4cout << G4endl; 1431 //G4cout << G4endl;
1241 if(i < nEtaL) { 1432 if(i < nEtaL) {
1242 //G4cout << " LShell:" <<G4endl; 1433 //G4cout << " LShell:" <<G4endl;
1243 for(j=0; j<nL; ++j) { 1434 for(j=0; j<nL; ++j) {
1244 1435
1245 if(j < 8) { 1436 if(j < 8) {
1246 bs = bls3[i][8-j]; 1437 bs = bls3[i][8-j];
1247 b = bll3[i][8-j]; 1438 b = bll3[i][8-j];
1248 } else if(j < 17) { 1439 } else if(j < 17) {
1249 bs = bls2[i][17-j]; 1440 bs = bls2[i][17-j];
1250 b = bll2[i][17-j]; 1441 b = bll2[i][17-j];
1251 } else { 1442 } else {
1252 bs = bls1[i][26-j]; 1443 bs = bls1[i][26-j];
1253 b = bll1[i][26-j]; 1444 b = bll1[i][26-j];
1254 } << 1445 }
1255 G4double c = SL[j]*loget + TL[j]; << 1446 G4double c = SL[j]*loget + TL[j];
1256 CL[j][i] = c - bs - 3.0*b; 1447 CL[j][i] = c - bs - 3.0*b;
1257 if(i == nEtaL-1) { 1448 if(i == nEtaL-1) {
1258 VL[j] = et*(et*CL[j][i] - UL[j]); << 1449 VL[j] = et*(et*CL[j][i] - UL[j]);
1259 //G4cout << "i= " << i << " j= " << << 1450 //G4cout << "i= " << i << " j= " << j
1260 // << " CL[j][i]= " << CL[ << 1451 // << " CL[j][i]= " << CL[j][i]
1261 // << " VL[j]= " << VL[j] < << 1452 // << " VL[j]= " << VL[j] << " b= " << b << " bs= " << bs
1262 // << " et= " << et << G4en << 1453 // << " et= " << et << G4endl;
1263 //" UL= " << UL[j] << " TL= " << TL << 1454 //" UL= " << UL[j] << " TL= " << TL[j] << " SL= " << SL[j] <<G4endl;
1264 } << 1455 }
1265 } 1456 }
1266 //G4cout << G4endl; 1457 //G4cout << G4endl;
1267 } 1458 }
1268 } 1459 }
>> 1460 const G4double hm[53] = {
>> 1461 12.0, 12.0, 12.0, 12.0, 11.9, 11.7, 11.5, 11.2, 10.8, 10.4,
>> 1462 10.0, 9.51, 8.97, 8.52, 8.03, 7.46, 6.95, 6.53, 6.18, 5.87,
>> 1463 5.61, 5.39, 5.19, 5.01, 4.86, 4.72, 4.62, 4.53, 4.44, 4.38,
>> 1464 4.32, 4.26, 4.20, 4.15, 4.1, 4.04, 4.00, 3.95, 3.93, 3.91,
>> 1465 3.90, 3.89, 3.89, 3.88, 3.88, 3.88, 3.88, 3.88, 3.89, 3.89,
>> 1466 3.90, 3.92, 3.93 };
>> 1467 const G4double hn[31] = {
>> 1468 75.5, 61.9, 52.2, 45.1, 39.6, 35.4, 31.9, 29.1, 27.2, 25.8,
>> 1469 24.5, 23.6, 22.7, 22.0, 21.4, 20.9, 20.5, 20.2, 19.9, 19.7,
>> 1470 19.5, 19.3, 19.2, 19.1, 18.4, 18.8, 18.7, 18.6, 18.5, 18.4,
>> 1471 18.2
>> 1472 };
>> 1473 for(i=0; i<53; ++i) {HM[i] = hm[i];}
>> 1474 for(i=0; i<31; ++i) {HN[i] = hn[i];}
1269 1475
1270 const G4double xzk[34] = { 11.7711, 1476 const G4double xzk[34] = { 11.7711,
1271 13.3669, 15.5762, 17.1715, 18.7667, 20.85 1477 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 1478 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 1479 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 << 1480 93.4549, 96.2753, 99.709};
1275 const G4double yzk[34] = { 0.70663, 1481 const G4double yzk[34] = { 0.70663,
1276 0.72033, 0.73651, 0.74647, 0.75518, 0.763 1482 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 1483 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 1484 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 << 1485 0.86891, 0.87011, 0.87381};
1280 1486
1281 const G4double xzl[36] = { 15.5102, 1487 const G4double xzl[36] = { 15.5102,
1282 16.7347, 17.9592, 19.551, 21.0204, 22.612 1488 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 1489 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 1490 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, << 1491 91.551, 94.2449, 96.449, 98.4082, 99.7551};
1286 const G4double yzl[36] = { 0.29875, 1492 const G4double yzl[36] = { 0.29875,
1287 0.31746, 0.33368, 0.35239, 0.36985, 0.387 1493 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 1494 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 1495 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 << 1496 0.6368, 0.64054, 0.64304, 0.64428, 0.64678};
>> 1497
>> 1498 ThetaK = new G4LPhysicsFreeVector(34, xzk[0], xzk[33]);
>> 1499 ThetaL = new G4LPhysicsFreeVector(36, xzl[0], xzl[35]);
>> 1500 for(i=0; i<34; ++i) { ThetaK->PutValues(i, xzk[i], yzk[i]); }
>> 1501 for(i=0; i<36; ++i) { ThetaL->PutValues(i, xzl[i], yzl[i]); }
>> 1502 ThetaK->SetSpline(true);
>> 1503 ThetaL->SetSpline(true);
>> 1504
>> 1505 const G4double coseb[14] = {0.0,0.05,0.1,0.15,0.2,0.3,0.4,0.5,0.6,0.8,
>> 1506 1.0,1.2,1.5,2.0};
>> 1507 const G4double cosxi[14] = {1.0000, 0.9905, 0.9631, 0.9208, 0.8680,
>> 1508 0.7478, 0.6303, 0.5290, 0.4471, 0.3323,
>> 1509 0.2610, 0.2145, 0.1696, 0.1261};
>> 1510 for(i=0; i<14; ++i) {
>> 1511 COSEB[i] = coseb[i];
>> 1512 COSXI[i] = cosxi[i];
>> 1513 }
1291 1514
1292 sThetaK = new G4PhysicsFreeVector(34, false << 1515 for(i=1; i<100; ++i) {
1293 for(i=0; i<34; ++i) { sThetaK->PutValues(i, << 1516 Z23[i] = std::pow(G4double(i),0.23);
1294 sThetaL = new G4PhysicsFreeVector(36, false << 1517 }
1295 for(i=0; i<36; ++i) { sThetaL->PutValues(i, <<
1296 } 1518 }
1297 1519
1298 //....oooOO0OOooo........oooOO0OOooo........o 1520 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1299 1521
1300 1522