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