source: trunk/source/processes/electromagnetic/lowenergy/src/G4eIonisationSpectrum.cc@ 966

Last change on this file since 966 was 961, checked in by garnier, 17 years ago

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