| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
>> 23 // $Id: G4eIonisationSpectrum.cc,v 1.12 2001/12/04 11:34:16 vnivanch Exp $
>> 24 // GEANT4 tag $Name: geant4-04-00 $
26 // 25 //
27 // ------------------------------------------- 26 // -------------------------------------------------------------------
28 // 27 //
29 // GEANT4 Class file 28 // GEANT4 Class file
30 // 29 //
31 // 30 //
32 // File name: G4eIonisationSpectrum 31 // File name: G4eIonisationSpectrum
33 // 32 //
34 // Author: V.Ivanchenko (Vladimir.Ivanc 33 // Author: V.Ivanchenko (Vladimir.Ivanchenko@cern.ch)
35 // << 34 //
36 // Creation date: 29 September 2001 35 // Creation date: 29 September 2001
37 // 36 //
38 // Modifications: << 37 // Modifications:
39 // 10.10.2001 MGP Revision to improv << 38 // 10.10.2001 MGP Revision to improve code quality and
40 // consistency with d 39 // consistency with design
41 // 02.11.2001 VI Optimize sampling << 40 // 02.11.2001 VI Optimize sampling of energy
42 // 29.11.2001 VI New parametrisatio << 41 // 29.11.2001 VI New parametrisation
43 // 19.04.2002 VI Add protection in <<
44 // 30.05.2002 VI Update to 24-param <<
45 // 11.07.2002 VI Fix in integration <<
46 // 23.03.2009 LP Added protection a <<
47 // faulty database fi <<
48 // 42 //
49 // ------------------------------------------- 43 // -------------------------------------------------------------------
50 // 44 //
51 45
52 #include "G4eIonisationSpectrum.hh" 46 #include "G4eIonisationSpectrum.hh"
53 #include "G4AtomicTransitionManager.hh" 47 #include "G4AtomicTransitionManager.hh"
54 #include "G4AtomicShell.hh" 48 #include "G4AtomicShell.hh"
>> 49 #include "G4eIonisationParameters.hh"
55 #include "G4DataVector.hh" 50 #include "G4DataVector.hh"
56 #include "Randomize.hh" 51 #include "Randomize.hh"
57 #include "G4PhysicalConstants.hh" << 52
58 #include "G4SystemOfUnits.hh" <<
59 #include "G4Exp.hh" <<
60 53
61 G4eIonisationSpectrum::G4eIonisationSpectrum() 54 G4eIonisationSpectrum::G4eIonisationSpectrum():G4VEnergySpectrum(),
62 lowestE(0.1*eV), 55 lowestE(0.1*eV),
63 factor(1.3), <<
64 iMax(24), <<
65 verbose(0) 56 verbose(0)
66 { 57 {
67 theParam = new G4eIonisationParameters(); 58 theParam = new G4eIonisationParameters();
68 } 59 }
69 60
70 61
71 G4eIonisationSpectrum::~G4eIonisationSpectrum( << 62 G4eIonisationSpectrum::~G4eIonisationSpectrum()
72 { 63 {
73 delete theParam; 64 delete theParam;
74 } 65 }
75 66
76 67
77 G4double G4eIonisationSpectrum::Probability(G4 << 68 G4double G4eIonisationSpectrum::Probability(G4int Z,
78 G4double tMin, << 69 G4double tMin,
79 G4double tMax, << 70 G4double tMax,
80 G4double e, 71 G4double e,
81 G4int shell, 72 G4int shell,
82 const G4ParticleDefinition* ) co << 73 const G4ParticleDefinition* part) const
83 { 74 {
84 // Please comment what Probability does and << 75 // Please comment what Probability does and what are the three
85 // functions mentioned below 76 // functions mentioned below
86 // Describe the algorithms used 77 // Describe the algorithms used
87 78
88 G4double eMax = MaxEnergyOfSecondaries(e); 79 G4double eMax = MaxEnergyOfSecondaries(e);
89 G4double t0 = std::max(tMin, lowestE); << 80 G4double t0 = G4std::max(tMin, lowestE);
90 G4double tm = std::min(tMax, eMax); << 81 G4double tm = G4std::min(tMax, eMax);
91 if(t0 >= tm) return 0.0; 82 if(t0 >= tm) return 0.0;
92 83
93 G4double bindingEnergy = (G4AtomicTransition 84 G4double bindingEnergy = (G4AtomicTransitionManager::Instance())->
94 Shell(Z, shell)->Bi 85 Shell(Z, shell)->BindingEnergy();
95 86
96 if(e <= bindingEnergy) return 0.0; << 87 G4double x1 = G4std::min(0.5,(t0 + bindingEnergy)/(e + bindingEnergy));
97 << 88 G4double x2 = G4std::min(0.5,(tm + bindingEnergy)/(e + bindingEnergy));
98 G4double energy = e + bindingEnergy; <<
99 <<
100 G4double x1 = std::min(0.5,(t0 + bindingEner <<
101 G4double x2 = std::min(0.5,(tm + bindingEner <<
102 89
103 if(verbose > 1 || (Z==4 && e>= 1.0 && e<= 0. << 90 if(verbose > 1) {
104 G4cout << "G4eIonisationSpectrum::Probabil 91 G4cout << "G4eIonisationSpectrum::Probability: Z= " << Z
105 << "; shell= " << shell 92 << "; shell= " << shell
106 << "; E(keV)= " << e/keV 93 << "; E(keV)= " << e/keV
107 << "; Eb(keV)= " << bindingEnergy/k << 94 << "; x1= " << x1
108 << "; x1= " << x1 << 95 << "; x2= " << x2
109 << "; x2= " << x2 <<
110 << G4endl; 96 << G4endl;
111 <<
112 } 97 }
113 98
>> 99 G4int iMax = 7;
114 G4DataVector p; 100 G4DataVector p;
115 101
116 // Access parameters 102 // Access parameters
117 for (G4int i=0; i<iMax; i++) << 103 for (G4int i=0; i<iMax; i++)
118 { 104 {
119 G4double x = theParam->Parameter(Z, shell, << 105 p.push_back(theParam->Parameter(Z, shell, i, e));
120 if(i<4) x /= energy; <<
121 p.push_back(x); <<
122 } 106 }
123 107
124 if(p[3] > 0.5) p[3] = 0.5; << 108 G4double g = (e + bindingEnergy)/electron_mass_c2 + 1.;
125 << 109 p.push_back((2.0*g - 1.0)/(g*g));
126 G4double gLocal = energy/electron_mass_c2 + << 110
127 p.push_back((2.0*gLocal - 1.0)/(gLocal*gLoca <<
128 <<
129 //Add protection against division by zero: a <<
130 //parameter p[3] appears in the denominator. <<
131 if (p[3] > 0) <<
132 p[iMax-1] = Function(p[3], p); <<
133 else <<
134 { <<
135 G4cout << "WARNING: G4eIonisationSpectru <<
136 << "parameter p[3] <= 0. G4LEDATA dabat <<
137 << Z << ". Please check and/or update i <<
138 } <<
139 <<
140 if(e >= 1. && e <= 0. && Z == 4) p.push_back <<
141 <<
142 <<
143 G4double val = IntSpectrum(x1, x2, p); 111 G4double val = IntSpectrum(x1, x2, p);
144 G4double x0 = (lowestE + bindingEnergy)/ene << 112 G4double x0 = (lowestE + bindingEnergy)/(e + bindingEnergy);
145 G4double nor = IntSpectrum(x0, 0.5, p); 113 G4double nor = IntSpectrum(x0, 0.5, p);
146 << 114
147 if(verbose > 1 || (Z==4 && e>= 1.0 && e<= 0. << 115 if(verbose > 1) {
148 G4cout << "tcut= " << tMin << 116 G4cout << "tcut= " << tMin
149 << "; tMax= " << tMax << 117 << "; tMax= " << tMax
150 << "; x0= " << x0 << 118 << "; x0= " << x0
151 << "; x1= " << x1 << 119 << "; x1= " << x1
152 << "; x2= " << x2 << 120 << "; x2= " << x2
153 << "; val= " << val << 121 << "; val= " << val
154 << "; nor= " << nor << 122 << "; nor= " << nor
155 << "; sum= " << p[0] << 123 << "; sum= " << p[0]
156 << "; a= " << p[1] << 124 << "; a= " << p[1]
157 << "; b= " << p[2] << 125 << "; b= " << p[2]
158 << "; c= " << p[3] << 126 << "; c= " << p[3]
159 << G4endl; 127 << G4endl;
160 if(shell == 1) G4cout << "============" << <<
161 } 128 }
162 129
163 p.clear(); 130 p.clear();
164 131
165 if(nor > 0.0) val /= nor; 132 if(nor > 0.0) val /= nor;
166 else val = 0.0; 133 else val = 0.0;
>> 134 if(val < 0.0) val = 0.0;
167 135
168 return val; << 136 return val;
169 } 137 }
170 138
171 139
172 G4double G4eIonisationSpectrum::AverageEnergy( 140 G4double G4eIonisationSpectrum::AverageEnergy(G4int Z,
173 G4double tMin, << 141 G4double tMin,
174 G4double tMax, << 142 G4double tMax,
175 G4double e, 143 G4double e,
176 G4int shell, 144 G4int shell,
177 const G4ParticleDefinition* ) << 145 const G4ParticleDefinition* part) const
178 { 146 {
179 // Please comment what AverageEnergy does an << 147 // Please comment what AverageEnergy does and what are the three
180 // functions mentioned below 148 // functions mentioned below
181 // Describe the algorithms used 149 // Describe the algorithms used
182 150
183 G4double eMax = MaxEnergyOfSecondaries(e); 151 G4double eMax = MaxEnergyOfSecondaries(e);
184 G4double t0 = std::max(tMin, lowestE); << 152 G4double t0 = G4std::max(tMin, lowestE);
185 G4double tm = std::min(tMax, eMax); << 153 G4double tm = G4std::min(tMax, eMax);
186 if(t0 >= tm) return 0.0; 154 if(t0 >= tm) return 0.0;
187 155
188 G4double bindingEnergy = (G4AtomicTransition 156 G4double bindingEnergy = (G4AtomicTransitionManager::Instance())->
189 Shell(Z, shell)->Bi 157 Shell(Z, shell)->BindingEnergy();
190 158
191 if(e <= bindingEnergy) return 0.0; << 159 G4double x1 = G4std::min(0.5,(t0 + bindingEnergy)/(e + bindingEnergy));
192 << 160 G4double x2 = G4std::min(0.5,(tm + bindingEnergy)/(e + bindingEnergy));
193 G4double energy = e + bindingEnergy; <<
194 <<
195 G4double x1 = std::min(0.5,(t0 + bindingEner <<
196 G4double x2 = std::min(0.5,(tm + bindingEner <<
197 161
198 if(verbose > 1) { 162 if(verbose > 1) {
199 G4cout << "G4eIonisationSpectrum::AverageE 163 G4cout << "G4eIonisationSpectrum::AverageEnergy: Z= " << Z
200 << "; shell= " << shell 164 << "; shell= " << shell
201 << "; E(keV)= " << e/keV 165 << "; E(keV)= " << e/keV
202 << "; bindingE(keV)= " << bindingEn 166 << "; bindingE(keV)= " << bindingEnergy/keV
203 << "; x1= " << x1 << 167 << "; x1= " << x1
204 << "; x2= " << x2 << 168 << "; x2= " << x2
205 << G4endl; 169 << G4endl;
206 } 170 }
207 171
>> 172 G4int iMax = 7;
208 G4DataVector p; 173 G4DataVector p;
209 174
210 // Access parameters 175 // Access parameters
211 for (G4int i=0; i<iMax; i++) << 176 for (G4int i=0; i<iMax; i++)
212 { 177 {
213 G4double x = theParam->Parameter(Z, shell, << 178 p.push_back(theParam->Parameter(Z, shell, i, e));
214 if(i<4) x /= energy; <<
215 p.push_back(x); <<
216 } 179 }
217 180
218 if(p[3] > 0.5) p[3] = 0.5; << 181 G4double g = (e + bindingEnergy)/electron_mass_c2 + 1.;
219 << 182 p.push_back((2.0*g - 1.0)/(g*g));
220 G4double gLocal2 = energy/electron_mass_c2 + << 183
221 p.push_back((2.0*gLocal2 - 1.0)/(gLocal2*gLo <<
222 <<
223 <<
224 //Add protection against division by zero: a <<
225 //parameter p[3] appears in the denominator. <<
226 if (p[3] > 0) <<
227 p[iMax-1] = Function(p[3], p); <<
228 else <<
229 { <<
230 G4cout << "WARNING: G4eIonisationSpectru <<
231 << "parameter p[3] <= 0. G4LEDATA dabat <<
232 << Z << ". Please check and/or update i <<
233 } <<
234 <<
235 G4double val = AverageValue(x1, x2, p); 184 G4double val = AverageValue(x1, x2, p);
236 G4double x0 = (lowestE + bindingEnergy)/ene << 185 G4double x0 = (lowestE + bindingEnergy)/(e + bindingEnergy);
237 G4double nor = IntSpectrum(x0, 0.5, p); 186 G4double nor = IntSpectrum(x0, 0.5, p);
238 val *= energy; << 187 val *= (e + bindingEnergy);
239 188
240 if(verbose > 1) { 189 if(verbose > 1) {
241 G4cout << "tcut(MeV)= " << tMin/MeV << 190 G4cout << "tcut(MeV)= " << tMin/MeV
242 << "; tMax(MeV)= " << tMax/MeV << 191 << "; tMax(MeV)= " << tMax/MeV
243 << "; x0= " << x0 << 192 << "; x0= " << x0
244 << "; x1= " << x1 << 193 << "; x1= " << x1
245 << "; x2= " << x2 << 194 << "; x2= " << x2
246 << "; val= " << val << 195 << "; val= " << val
247 << "; nor= " << nor << 196 << "; nor= " << nor
248 << "; sum= " << p[0] << 197 << "; sum= " << p[0]
249 << "; a= " << p[1] << 198 << "; a= " << p[1]
250 << "; b= " << p[2] << 199 << "; b= " << p[2]
251 << "; c= " << p[3] << 200 << "; c= " << p[3]
252 << G4endl; 201 << G4endl;
253 } 202 }
254 203
255 p.clear(); 204 p.clear();
256 205
257 if(nor > 0.0) val /= nor; 206 if(nor > 0.0) val /= nor;
258 else val = 0.0; 207 else val = 0.0;
>> 208 if(val < 0.0) val = 0.0;
259 209
260 return val; << 210 return val;
261 } 211 }
262 212
263 213
264 G4double G4eIonisationSpectrum::SampleEnergy(G 214 G4double G4eIonisationSpectrum::SampleEnergy(G4int Z,
265 G4double tMin, << 215 G4double tMin,
266 G4double tMax, << 216 G4double tMax,
267 G4double e, 217 G4double e,
268 G4int shell, 218 G4int shell,
269 const G4ParticleDefinition* ) c << 219 const G4ParticleDefinition* part) const
270 { 220 {
271 // Please comment what SampleEnergy does 221 // Please comment what SampleEnergy does
272 G4double tDelta = 0.0; 222 G4double tDelta = 0.0;
273 G4double t0 = std::max(tMin, lowestE); << 223 G4double t0 = G4std::max(tMin, lowestE);
274 G4double tm = std::min(tMax, MaxEnergyOfSeco << 224 G4double tm = G4std::min(tMax, MaxEnergyOfSecondaries(e));
275 if(t0 > tm) return tDelta; 225 if(t0 > tm) return tDelta;
276 226
277 G4double bindingEnergy = (G4AtomicTransition 227 G4double bindingEnergy = (G4AtomicTransitionManager::Instance())->
278 Shell(Z, shell)->Bi 228 Shell(Z, shell)->BindingEnergy();
279 229
280 if(e <= bindingEnergy) return 0.0; << 230 G4double x1 = G4std::min(0.5,(t0 + bindingEnergy)/(e + bindingEnergy));
281 << 231 G4double x2 = G4std::min(0.5,(tm + bindingEnergy)/(e + bindingEnergy));
282 G4double energy = e + bindingEnergy; <<
283 <<
284 G4double x1 = std::min(0.5,(t0 + bindingEner <<
285 G4double x2 = std::min(0.5,(tm + bindingEner <<
286 if(x1 >= x2) return tDelta; 232 if(x1 >= x2) return tDelta;
287 233
288 if(verbose > 1) { 234 if(verbose > 1) {
289 G4cout << "G4eIonisationSpectrum::SampleEn 235 G4cout << "G4eIonisationSpectrum::SampleEnergy: Z= " << Z
290 << "; shell= " << shell 236 << "; shell= " << shell
291 << "; E(keV)= " << e/keV 237 << "; E(keV)= " << e/keV
292 << G4endl; 238 << G4endl;
293 } 239 }
294 240
295 // Access parameters 241 // Access parameters
>> 242 G4int iMax = 7;
296 G4DataVector p; 243 G4DataVector p;
297 244
298 // Access parameters 245 // Access parameters
299 for (G4int i=0; i<iMax; i++) << 246 for (G4int i=0; i<iMax; i++)
300 { 247 {
301 G4double x = theParam->Parameter(Z, shell, << 248 p.push_back(theParam->Parameter(Z, shell, i, e));
302 if(i<4) x /= energy; <<
303 p.push_back(x); <<
304 } 249 }
305 250
306 if(p[3] > 0.5) p[3] = 0.5; << 251 G4double g = (e + bindingEnergy)/electron_mass_c2 + 1.;
307 << 252 p.push_back((2.0*g - 1.0)/(g*g));
308 G4double gLocal3 = energy/electron_mass_c2 + <<
309 p.push_back((2.0*gLocal3 - 1.0)/(gLocal3*gLo <<
310 253
311 254
312 //Add protection against division by zero: a <<
313 //parameter p[3] appears in the denominator. <<
314 if (p[3] > 0) <<
315 p[iMax-1] = Function(p[3], p); <<
316 else <<
317 { <<
318 G4cout << "WARNING: G4eIonisationSpectru <<
319 << "parameter p[3] <= 0. G4LEDATA dabat <<
320 << Z << ". Please check and/or update i <<
321 } <<
322 <<
323 G4double aria1 = 0.0; 255 G4double aria1 = 0.0;
324 G4double a1 = std::max(x1,p[1]); << 256 G4double a1 = G4std::min(x1,p[6]);
325 G4double a2 = std::min(x2,p[3]); << 257 G4double a2 = G4std::min(x2,p[6]);
326 if(a1 < a2) aria1 = IntSpectrum(a1, a2, p); 258 if(a1 < a2) aria1 = IntSpectrum(a1, a2, p);
327 G4double aria2 = 0.0; 259 G4double aria2 = 0.0;
328 G4double a3 = std::max(x1,p[3]); << 260 G4double a3 = G4std::max(x1,p[6]);
329 G4double a4 = x2; << 261 G4double a4 = G4std::max(x2,p[6]);
330 if(a3 < a4) aria2 = IntSpectrum(a3, a4, p); 262 if(a3 < a4) aria2 = IntSpectrum(a3, a4, p);
331 263
332 G4double aria = (aria1 + aria2)*G4UniformRan 264 G4double aria = (aria1 + aria2)*G4UniformRand();
333 G4double amaj, fun, q, x, z1, z2, dx, dx1; << 265 G4double amaj, fun, q, x;
334 266
335 //======= First aria to sample ===== 267 //======= First aria to sample =====
336 268
337 if(aria <= aria1) { << 269 if(aria <= aria1) {
338 270
339 amaj = p[4]; 271 amaj = p[4];
340 for (G4int j=5; j<iMax; j++) { <<
341 if(p[j] > amaj) amaj = p[j]; <<
342 } <<
343 <<
344 a1 = 1./a1; 272 a1 = 1./a1;
345 a2 = 1./a2; 273 a2 = 1./a2;
346 274
347 G4int i; <<
348 do { <<
349 <<
350 x = 1./(a2 + G4UniformRand()*(a1 - a2)); <<
351 z1 = p[1]; <<
352 z2 = p[3]; <<
353 dx = (p[2] - p[1]) / 3.0; <<
354 dx1= G4Exp(std::log(p[3]/p[2]) / 16.0); <<
355 for (i=4; i<iMax-1; i++) { <<
356 <<
357 if (i < 7) { <<
358 z2 = z1 + dx; <<
359 } else if(iMax-2 == i) { <<
360 z2 = p[3]; <<
361 break; <<
362 } else { <<
363 z2 = z1*dx1; <<
364 } <<
365 if(x >= z1 && x <= z2) break; <<
366 z1 = z2; <<
367 } <<
368 fun = p[i] + (x - z1) * (p[i+1] - p[i])/ <<
369 <<
370 if(fun > amaj) { <<
371 G4cout << "WARNING in G4eIonisationS <<
372 << " Majoranta " << amaj <<
373 << " < " << fun <<
374 << " in the first aria at x= <<
375 << G4endl; <<
376 } <<
377 <<
378 q = amaj*G4UniformRand(); <<
379 <<
380 } while (q >= fun); <<
381 <<
382 //======= Second aria to sample ===== 275 //======= Second aria to sample =====
383 276
384 } else { 277 } else {
385 278
386 amaj = std::max(p[iMax-1], Function(0.5, p << 279 amaj = p[5];
387 a1 = 1./a3; 280 a1 = 1./a3;
388 a2 = 1./a4; 281 a2 = 1./a4;
>> 282 }
>> 283 amaj *= 1.25;
389 284
390 do { << 285 do {
391 286
392 x = 1./(a2 + G4UniformRand()*(a1 - a2)); << 287 x = 1./(a2 + G4UniformRand()*(a1 - a2));
393 fun = Function(x, p); << 288 fun = Function(x, p);
394 289
395 if(fun > amaj) { << 290 if(fun > amaj) {
396 G4cout << "WARNING in G4eIonisationS << 291 G4cout << "WARNING in G4eIonisationSpectrum::SampleEnergy:"
397 << " Majoranta " << amaj << 292 << " Majoranta " << amaj
398 << " < " << fun 293 << " < " << fun
399 << " in the second aria at x= <<
400 << G4endl; 294 << G4endl;
401 } << 295 }
402 296
403 q = amaj*G4UniformRand(); << 297 q = amaj*G4UniformRand();
404 298
405 } while (q >= fun); << 299 } while (q >= fun);
406 <<
407 } <<
408 300
409 p.clear(); 301 p.clear();
410 302
411 tDelta = x*energy - bindingEnergy; << 303 tDelta = x*(e + bindingEnergy) - bindingEnergy;
412 304
413 if(verbose > 1) { 305 if(verbose > 1) {
414 G4cout << "tcut(MeV)= " << tMin/MeV << 306 G4cout << "tcut(MeV)= " << tMin/MeV
415 << "; tMax(MeV)= " << tMax/MeV << 307 << "; tMax(MeV)= " << tMax/MeV
416 << "; x1= " << x1 << 308 << "; x1= " << x1
417 << "; x2= " << x2 << 309 << "; x2= " << x2
418 << "; a1= " << a1 << 310 << "; a1= " << a1
419 << "; a2= " << a2 << 311 << "; a2= " << a2
420 << "; x= " << x << 312 << "; x= " << x
421 << "; be= " << bindingEnergy << 313 << "; be= " << bindingEnergy
422 << "; e= " << e << 314 << "; e= " << e
423 << "; tDelta= " << tDelta << 315 << "; tDelta= " << tDelta
424 << G4endl; 316 << G4endl;
425 } 317 }
426 318
427 319
428 return tDelta; << 320 return tDelta;
429 } 321 }
430 322
431 323
432 G4double G4eIonisationSpectrum::IntSpectrum(G4 << 324 G4double G4eIonisationSpectrum::IntSpectrum(G4double xMin,
433 G4double xMax, 325 G4double xMax,
434 const G4DataVector& p) const 326 const G4DataVector& p) const
435 { 327 {
436 // Please comment what IntSpectrum does 328 // Please comment what IntSpectrum does
437 G4double sum = 0.0; << 329 G4double x1 = 1./xMin;
438 if(xMin >= xMax) return sum; << 330 G4double x2 = 1./xMax;
439 << 331 G4double x = x1 - x2 - p[7]*log(xMax/xMin) + (1. - p[7])*(xMax - xMin)
440 G4double x1, x2, xs1, xs2, y1, y2, ys1, ys2, << 332 + 1./(1. - xMax) - 1./(1. - xMin)
441 << 333 + p[7]*log((1. - xMax)/(1. - xMin))
442 // Integral over interpolation aria << 334 + 0.5*p[1]*p[3]*(x1*x1 - x2*x2);
443 if(xMin < p[3]) { << 335
444 << 336 if(x < 0.0) x = 0.0;
445 x1 = p[1]; << 337 return x;
446 y1 = p[4]; << 338 }
447 <<
448 G4double dx = (p[2] - p[1]) / 3.0; <<
449 G4double dx1= G4Exp(std::log(p[3]/p[2]) / <<
450 <<
451 for (size_t i=0; i<19; i++) { <<
452 <<
453 q = 0.0; <<
454 if (i < 3) { <<
455 x2 = x1 + dx; <<
456 } else if(18 == i) { <<
457 x2 = p[3]; <<
458 } else { <<
459 x2 = x1*dx1; <<
460 } <<
461 <<
462 y2 = p[5 + i]; <<
463 <<
464 if (xMax <= x1) { <<
465 break; <<
466 } else if (xMin < x2) { <<
467 <<
468 xs1 = x1; <<
469 xs2 = x2; <<
470 ys1 = y1; <<
471 ys2 = y2; <<
472 <<
473 if (x2 > x1) { <<
474 if (xMin > x1) { <<
475 xs1 = xMin; <<
476 ys1 += (xs1 - x1)*(y2 - y1)/(x2 - <<
477 } <<
478 if (xMax < x2) { <<
479 xs2 = xMax; <<
480 ys2 += (xs2 - x2)*(y1 - y2)/(x1 - <<
481 } <<
482 if (xs2 > xs1) { <<
483 q = (ys1*xs2 - ys2*xs1)/(xs1*xs2) <<
484 + std::log(xs2/xs1)*(ys2 - ys1) <<
485 sum += q; <<
486 if(p.size() == 26) G4cout << "i= " <<
487 } <<
488 } <<
489 } <<
490 x1 = x2; <<
491 y1 = y2; <<
492 } <<
493 } <<
494 <<
495 // Integral over aria with parametrised form <<
496 339
497 x1 = std::max(xMin, p[3]); <<
498 if(x1 >= xMax) return sum; <<
499 x2 = xMax; <<
500 <<
501 xs1 = 1./x1; <<
502 xs2 = 1./x2; <<
503 q = (xs1 - xs2)*(1.0 - p[0]) <<
504 - p[iMax]*std::log(x2/x1) <<
505 + (1. - p[iMax])*(x2 - x1) <<
506 + 1./(1. - x2) - 1./(1. - x1) <<
507 + p[iMax]*std::log((1. - x2)/(1. - x1)) <<
508 + 0.25*p[0]*(xs1*xs1 - xs2*xs2); <<
509 sum += q; <<
510 if(p.size() == 26) G4cout << "param... q= " <<
511 340
512 return sum; << 341 G4double G4eIonisationSpectrum::AverageValue(G4double xMin,
513 } <<
514 <<
515 <<
516 G4double G4eIonisationSpectrum::AverageValue(G <<
517 G4double xMax, 342 G4double xMax,
518 const G4DataVector& p) const 343 const G4DataVector& p) const
519 { 344 {
520 G4double sum = 0.0; <<
521 if(xMin >= xMax) return sum; <<
522 345
523 G4double x1, x2, xs1, xs2, y1, y2, ys1, ys2; << 346 // G4double x1 = 1.;
>> 347 // G4double x2 = 1.;
>> 348 G4double x = log(xMax/xMin)
>> 349 + 0.5*(1. - p[7])*(xMax*xMax - xMin*xMin)
>> 350 + 1./(1. - xMax) - 1./(1. - xMin)
>> 351 + (1. + p[7])*log((1. - xMax)/(1. - xMin))
>> 352 + p[1]*p[3]*(1./xMin - 1./xMax);
>> 353
>> 354 if(x < 0.0) x = 0.0;
>> 355 return x;
>> 356 }
524 357
525 // Integral over interpolation aria <<
526 if(xMin < p[3]) { <<
527 358
528 x1 = p[1]; << 359 G4double G4eIonisationSpectrum::Function(G4double x,
529 y1 = p[4]; << 360 const G4DataVector& p) const
530 << 361 {
531 G4double dx = (p[2] - p[1]) / 3.0; << 362 // Please comment what Function does
532 G4double dx1= G4Exp(std::log(p[3]/p[2]) / <<
533 <<
534 for (size_t i=0; i<19; i++) { <<
535 <<
536 if (i < 3) { <<
537 x2 = x1 + dx; <<
538 } else if(18 == i) { <<
539 x2 = p[3]; <<
540 } else { <<
541 x2 = x1*dx1; <<
542 } <<
543 <<
544 y2 = p[5 + i]; <<
545 <<
546 if (xMax <= x1) { <<
547 break; <<
548 } else if (xMin < x2) { <<
549 <<
550 xs1 = x1; <<
551 xs2 = x2; <<
552 ys1 = y1; <<
553 ys2 = y2; <<
554 <<
555 if (x2 > x1) { <<
556 if (xMin > x1) { <<
557 xs1 = xMin; <<
558 ys1 += (xs1 - x1)*(y2 - y1)/(x2 - <<
559 } <<
560 if (xMax < x2) { <<
561 xs2 = xMax; <<
562 ys2 += (xs2 - x2)*(y1 - y2)/(x1 - <<
563 } <<
564 if (xs2 > xs1) { <<
565 sum += std::log(xs2/xs1)*(ys1*xs2 <<
566 + ys2 - ys1; <<
567 } <<
568 } <<
569 } <<
570 x1 = x2; <<
571 y1 = y2; <<
572 363
573 } <<
574 } <<
575 364
576 // Integral over aria with parametrised form << 365 // G4double x1 = 1.0;
>> 366 G4double f = 1.0 - p[7]*x + x*x*(1.0 - p[7]
>> 367 + (1.0/(1.0 - x) - p[7])/(1.0 - x) )
>> 368 + p[1]*p[3]/x;
577 369
578 x1 = std::max(xMin, p[3]); << 370 if(f < 0.0) f = 0.0;
579 if(x1 >= xMax) return sum; << 371 return f;
580 x2 = xMax; << 372 }
581 <<
582 xs1 = 1./x1; <<
583 xs2 = 1./x2; <<
584 <<
585 sum += std::log(x2/x1)*(1.0 - p[0]) <<
586 + 0.5*(1. - p[iMax])*(x2*x2 - x1*x1) <<
587 + 1./(1. - x2) - 1./(1. - x1) <<
588 + (1. + p[iMax])*std::log((1. - x2)/(1 <<
589 + 0.5*p[0]*(xs1 - xs2); <<
590 373
591 return sum; << 374 G4double G4eIonisationSpectrum::Excitation(G4int Z, G4double e) const
>> 375 {
>> 376 return theParam->Excitation(Z, e);
592 } 377 }
593 378
594 << 379 void G4eIonisationSpectrum::PrintData() const
595 void G4eIonisationSpectrum::PrintData() const <<
596 { 380 {
597 theParam->PrintData(); 381 theParam->PrintData();
598 } 382 }
599 383
600 G4double G4eIonisationSpectrum::MaxEnergyOfSec << 384
601 G4int, // Z = 0, << 385
602 const G4ParticleDefinition* << 386
603 { << 387
604 return 0.5 * kineticEnergy; <<
605 } <<
606 388