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