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