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