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