Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4eIonisationSpectrum.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 //
 27 // -------------------------------------------------------------------
 28 //
 29 // GEANT4 Class file
 30 //
 31 //
 32 // File name:     G4eIonisationSpectrum
 33 //
 34 // Author:        V.Ivanchenko (Vladimir.Ivanchenko@cern.ch)
 35 //
 36 // Creation date: 29 September 2001
 37 //
 38 // Modifications:
 39 // 10.10.2001 MGP           Revision to improve code quality and
 40 //                          consistency with design
 41 // 02.11.2001 VI            Optimize sampling of energy
 42 // 29.11.2001 VI            New parametrisation
 43 // 19.04.2002 VI            Add protection in case of energy below binding
 44 // 30.05.2002 VI            Update to 24-parameters data
 45 // 11.07.2002 VI            Fix in integration over spectrum
 46 // 23.03.2009 LP            Added protection against division by zero (for
 47 //                          faulty database files), Bug Report 1042
 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 #include "G4PhysicalConstants.hh"
 58 #include "G4SystemOfUnits.hh"
 59 #include "G4Exp.hh"
 60 
 61 G4eIonisationSpectrum::G4eIonisationSpectrum():G4VEnergySpectrum(),
 62   lowestE(0.1*eV),
 63   factor(1.3),
 64   iMax(24),
 65   verbose(0)
 66 {
 67   theParam = new G4eIonisationParameters();
 68 }
 69 
 70 
 71 G4eIonisationSpectrum::~G4eIonisationSpectrum()
 72 {
 73   delete theParam;
 74 }
 75 
 76 
 77 G4double G4eIonisationSpectrum::Probability(G4int Z,
 78                     G4double tMin,
 79               G4double tMax,
 80               G4double e,
 81               G4int shell,
 82               const G4ParticleDefinition* ) const
 83 {
 84   // Please comment what Probability does and what are the three
 85   // functions mentioned below
 86   // Describe the algorithms used
 87 
 88   G4double eMax = MaxEnergyOfSecondaries(e);
 89   G4double t0 = std::max(tMin, lowestE);
 90   G4double tm = std::min(tMax, eMax);
 91   if(t0 >= tm) return 0.0;
 92 
 93   G4double bindingEnergy = (G4AtomicTransitionManager::Instance())->
 94                            Shell(Z, shell)->BindingEnergy();
 95 
 96   if(e <= bindingEnergy) return 0.0;
 97 
 98   G4double energy = e + bindingEnergy;
 99 
100   G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
101   G4double x2 = std::min(0.5,(tm + bindingEnergy)/energy);
102 
103   if(verbose > 1 || (Z==4 && e>= 1.0 && e<= 0.0)) {
104     G4cout << "G4eIonisationSpectrum::Probability: Z= " << Z
105            << "; shell= " << shell
106            << "; E(keV)= " << e/keV
107            << "; Eb(keV)= " << bindingEnergy/keV
108            << "; x1= " << x1
109            << "; x2= " << x2
110            << G4endl;
111 
112   }
113 
114   G4DataVector p;
115 
116   // Access parameters
117   for (G4int i=0; i<iMax; i++)
118   {
119     G4double x = theParam->Parameter(Z, shell, i, e);
120     if(i<4) x /= energy;
121     p.push_back(x);
122   }
123 
124   if(p[3] > 0.5) p[3] = 0.5;
125 
126   G4double gLocal = energy/electron_mass_c2 + 1.;
127   p.push_back((2.0*gLocal - 1.0)/(gLocal*gLocal));
128 
129   //Add protection against division by zero: actually in Function()
130   //parameter p[3] appears in the denominator. Bug report 1042
131   if (p[3] > 0)
132     p[iMax-1] = Function(p[3], p);
133   else
134     {
135       G4cout << "WARNING: G4eIonisationSpectrum::Probability "
136        << "parameter p[3] <= 0. G4LEDATA dabatase might be corrupted for Z = "
137        << Z << ". Please check and/or update it " << G4endl;
138     }
139 
140   if(e >= 1. && e <= 0. && Z == 4) p.push_back(0.0);
141 
142 
143   G4double val = IntSpectrum(x1, x2, p);
144   G4double x0  = (lowestE + bindingEnergy)/energy;
145   G4double nor = IntSpectrum(x0, 0.5, p);
146 
147   if(verbose > 1 || (Z==4 && e>= 1.0 && e<= 0.0)) {
148     G4cout << "tcut= " << tMin
149            << "; tMax= " << tMax
150            << "; x0= " << x0
151            << "; x1= " << x1
152            << "; x2= " << x2
153            << "; val= " << val
154            << "; nor= " << nor
155            << "; sum= " << p[0]
156            << "; a= " << p[1]
157            << "; b= " << p[2]
158            << "; c= " << p[3]
159            << G4endl;
160     if(shell == 1) G4cout << "============" << G4endl;
161   }
162 
163   p.clear();
164 
165   if(nor > 0.0) val /= nor;
166   else          val  = 0.0;
167 
168   return val;
169 }
170 
171 
172 G4double G4eIonisationSpectrum::AverageEnergy(G4int Z,
173                       G4double tMin,
174                 G4double tMax,
175                 G4double e,
176                 G4int shell,
177                 const G4ParticleDefinition* ) const
178 {
179   // Please comment what AverageEnergy does and what are the three
180   // functions mentioned below
181   // Describe the algorithms used
182 
183   G4double eMax = MaxEnergyOfSecondaries(e);
184   G4double t0 = std::max(tMin, lowestE);
185   G4double tm = std::min(tMax, eMax);
186   if(t0 >= tm) return 0.0;
187 
188   G4double bindingEnergy = (G4AtomicTransitionManager::Instance())->
189                            Shell(Z, shell)->BindingEnergy();
190 
191   if(e <= bindingEnergy) return 0.0;
192 
193   G4double energy = e + bindingEnergy;
194 
195   G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
196   G4double x2 = std::min(0.5,(tm + bindingEnergy)/energy);
197 
198   if(verbose > 1) {
199     G4cout << "G4eIonisationSpectrum::AverageEnergy: Z= " << Z
200            << "; shell= " << shell
201            << "; E(keV)= " << e/keV
202            << "; bindingE(keV)= " << bindingEnergy/keV
203            << "; x1= " << x1
204            << "; x2= " << x2
205            << G4endl;
206   }
207 
208   G4DataVector p;
209 
210   // Access parameters
211   for (G4int i=0; i<iMax; i++)
212   {
213     G4double x = theParam->Parameter(Z, shell, i, e);
214     if(i<4) x /= energy;
215     p.push_back(x);
216   }
217 
218   if(p[3] > 0.5) p[3] = 0.5;
219 
220   G4double gLocal2 = energy/electron_mass_c2 + 1.;
221   p.push_back((2.0*gLocal2 - 1.0)/(gLocal2*gLocal2));
222 
223 
224   //Add protection against division by zero: actually in Function()
225   //parameter p[3] appears in the denominator. Bug report 1042
226   if (p[3] > 0)
227     p[iMax-1] = Function(p[3], p);
228   else
229     {
230       G4cout << "WARNING: G4eIonisationSpectrum::AverageEnergy "
231        << "parameter p[3] <= 0. G4LEDATA dabatase might be corrupted for Z = "
232        << Z << ". Please check and/or update it " << G4endl;
233     }
234 
235   G4double val = AverageValue(x1, x2, p);
236   G4double x0  = (lowestE + bindingEnergy)/energy;
237   G4double nor = IntSpectrum(x0, 0.5, p);
238   val *= energy;
239 
240   if(verbose > 1) {
241     G4cout << "tcut(MeV)= " << tMin/MeV
242            << "; tMax(MeV)= " << tMax/MeV
243            << "; x0= " << x0
244            << "; x1= " << x1
245            << "; x2= " << x2
246            << "; val= " << val
247            << "; nor= " << nor
248            << "; sum= " << p[0]
249            << "; a= " << p[1]
250            << "; b= " << p[2]
251            << "; c= " << p[3]
252            << G4endl;
253   }
254 
255   p.clear();
256 
257   if(nor > 0.0) val /= nor;
258   else          val  = 0.0;
259 
260   return val;
261 }
262 
263 
264 G4double G4eIonisationSpectrum::SampleEnergy(G4int Z,
265                G4double tMin,
266                G4double tMax,
267                      G4double e,
268                      G4int shell,
269                const G4ParticleDefinition* ) const
270 {
271   // Please comment what SampleEnergy does
272   G4double tDelta = 0.0;
273   G4double t0 = std::max(tMin, lowestE);
274   G4double tm = std::min(tMax, MaxEnergyOfSecondaries(e));
275   if(t0 > tm) return tDelta;
276 
277   G4double bindingEnergy = (G4AtomicTransitionManager::Instance())->
278                            Shell(Z, shell)->BindingEnergy();
279 
280   if(e <= bindingEnergy) return 0.0;
281 
282   G4double energy = e + bindingEnergy;
283 
284   G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
285   G4double x2 = std::min(0.5,(tm + bindingEnergy)/energy);
286   if(x1 >= x2) return tDelta;
287 
288   if(verbose > 1) {
289     G4cout << "G4eIonisationSpectrum::SampleEnergy: Z= " << Z
290            << "; shell= " << shell
291            << "; E(keV)= " << e/keV
292            << G4endl;
293   }
294 
295   // Access parameters
296   G4DataVector p;
297 
298   // Access parameters
299   for (G4int i=0; i<iMax; i++)
300   {
301     G4double x = theParam->Parameter(Z, shell, i, e);
302     if(i<4) x /= energy;
303     p.push_back(x);
304   }
305 
306   if(p[3] > 0.5) p[3] = 0.5;
307 
308   G4double gLocal3 = energy/electron_mass_c2 + 1.;
309   p.push_back((2.0*gLocal3 - 1.0)/(gLocal3*gLocal3));
310 
311 
312   //Add protection against division by zero: actually in Function()
313   //parameter p[3] appears in the denominator. Bug report 1042
314   if (p[3] > 0)
315     p[iMax-1] = Function(p[3], p);
316   else
317     {
318       G4cout << "WARNING: G4eIonisationSpectrum::SampleSpectrum "
319        << "parameter p[3] <= 0. G4LEDATA dabatase might be corrupted for Z = "
320        << Z << ". Please check and/or update it " << G4endl;
321     }
322 
323   G4double aria1 = 0.0;
324   G4double a1 = std::max(x1,p[1]);
325   G4double a2 = std::min(x2,p[3]);
326   if(a1 < a2) aria1 = IntSpectrum(a1, a2, p);
327   G4double aria2 = 0.0;
328   G4double a3 = std::max(x1,p[3]);
329   G4double a4 = x2;
330   if(a3 < a4) aria2 = IntSpectrum(a3, a4, p);
331 
332   G4double aria = (aria1 + aria2)*G4UniformRand();
333   G4double amaj, fun, q, x, z1, z2, dx, dx1;
334 
335   //======= First aria to sample =====
336 
337   if(aria <= aria1) {
338 
339     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;
345     a2 = 1./a2;
346 
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])/(z2 - z1);
369 
370       if(fun > amaj) {
371           G4cout << "WARNING in G4eIonisationSpectrum::SampleEnergy:"
372                  << " Majoranta " << amaj
373                  << " < " << fun
374                  << " in the first aria at x= " << x
375                  << G4endl;
376       }
377 
378       q = amaj*G4UniformRand();
379 
380     } while (q >= fun);
381 
382   //======= Second aria to sample =====
383 
384   } else {
385 
386     amaj = std::max(p[iMax-1], Function(0.5, p)) * factor;
387     a1 = 1./a3;
388     a2 = 1./a4;
389 
390     do {
391 
392       x = 1./(a2 + G4UniformRand()*(a1 - a2));
393       fun = Function(x, p);
394 
395       if(fun > amaj) {
396           G4cout << "WARNING in G4eIonisationSpectrum::SampleEnergy:"
397                  << " Majoranta " << amaj
398                  << " < " << fun
399                  << " in the second aria at x= " << x
400                  << G4endl;
401       }
402 
403       q = amaj*G4UniformRand();
404 
405     } while (q >= fun);
406 
407   }
408 
409   p.clear();
410 
411   tDelta = x*energy - bindingEnergy;
412 
413   if(verbose > 1) {
414     G4cout << "tcut(MeV)= " << tMin/MeV
415            << "; tMax(MeV)= " << tMax/MeV
416            << "; x1= " << x1
417            << "; x2= " << x2
418            << "; a1= " << a1
419            << "; a2= " << a2
420            << "; x= " << x
421            << "; be= " << bindingEnergy
422            << "; e= " << e
423            << "; tDelta= " << tDelta
424            << G4endl;
425   }
426 
427 
428   return tDelta;
429 }
430 
431 
432 G4double G4eIonisationSpectrum::IntSpectrum(G4double xMin,
433               G4double xMax,
434               const G4DataVector& p) const
435 {
436   // Please comment what IntSpectrum does
437   G4double sum = 0.0;
438   if(xMin >= xMax) return sum;
439 
440   G4double x1, x2, xs1, xs2, y1, y2, ys1, ys2, q;
441 
442   // Integral over interpolation aria
443   if(xMin < p[3]) {
444 
445     x1 = p[1];
446     y1 = p[4];
447 
448     G4double dx = (p[2] - p[1]) / 3.0;
449     G4double dx1= G4Exp(std::log(p[3]/p[2]) / 16.0);
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 - x1);
477     }
478           if (xMax < x2) {
479             xs2 = xMax;
480             ys2 += (xs2 - x2)*(y1 - y2)/(x1 - x2);
481     }
482           if (xs2 > xs1) {
483             q = (ys1*xs2 - ys2*xs1)/(xs1*xs2)
484               +  std::log(xs2/xs1)*(ys2 - ys1)/(xs2 - xs1);
485             sum += q;
486             if(p.size() == 26) G4cout << "i= " << i << "  q= " << q << " sum= " << sum << G4endl;
487     }
488   }
489       }
490       x1 = x2;
491       y1 = y2;
492     }
493   }
494 
495   // Integral over aria with parametrised formula
496 
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= " << q <<  " sum= " << sum << G4endl;
511 
512   return sum;
513 }
514 
515 
516 G4double G4eIonisationSpectrum::AverageValue(G4double xMin,
517                      G4double xMax,
518                const G4DataVector& p) const
519 {
520   G4double sum = 0.0;
521   if(xMin >= xMax) return sum;
522 
523   G4double x1, x2, xs1, xs2, y1, y2, ys1, ys2;
524 
525   // Integral over interpolation aria
526   if(xMin < p[3]) {
527 
528     x1 = p[1];
529     y1 = p[4];
530 
531     G4double dx = (p[2] - p[1]) / 3.0;
532     G4double dx1= G4Exp(std::log(p[3]/p[2]) / 16.0);
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 - x1);
559     }
560           if (xMax < x2) {
561             xs2 = xMax;
562             ys2 += (xs2 - x2)*(y1 - y2)/(x1 - x2);
563     }
564           if (xs2 > xs1) {
565             sum += std::log(xs2/xs1)*(ys1*xs2 - ys2*xs1)/(xs2 - xs1)
566                 +  ys2 - ys1;
567     }
568   }
569       }
570       x1 = x2;
571       y1 = y2;
572 
573     }
574   }
575 
576   // Integral over aria with parametrised formula
577 
578   x1 = std::max(xMin, p[3]);
579   if(x1 >= xMax) return sum;
580   x2 = xMax;
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. - x1))
589         + 0.5*p[0]*(xs1 - xs2);
590 
591   return sum;
592 }
593 
594 
595 void G4eIonisationSpectrum::PrintData() const
596 {
597   theParam->PrintData();
598 }
599 
600 G4double G4eIonisationSpectrum::MaxEnergyOfSecondaries(G4double kineticEnergy,
601                    G4int, // Z = 0,
602                    const G4ParticleDefinition* ) const
603 {
604   return 0.5 * kineticEnergy;
605 }
606