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

Last change on this file since 989 was 961, checked in by garnier, 15 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.