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