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