Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/eRosita/physics/src/G4RDeIonisationSpectrum.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:     G4RDeIonisationSpectrum
 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 //
 47 // -------------------------------------------------------------------
 48 //
 49 
 50 #include "G4RDeIonisationSpectrum.hh"
 51 #include "G4RDAtomicTransitionManager.hh"
 52 #include "G4RDAtomicShell.hh"
 53 #include "G4PhysicalConstants.hh"
 54 #include "G4SystemOfUnits.hh"
 55 #include "G4DataVector.hh"
 56 #include "Randomize.hh"
 57 
 58 
 59 G4RDeIonisationSpectrum::G4RDeIonisationSpectrum():G4RDVEnergySpectrum(),
 60   lowestE(0.1*eV),
 61   factor(1.3),
 62   iMax(24),
 63   verbose(0)
 64 {
 65   theParam = new G4RDeIonisationParameters();
 66 }
 67 
 68 
 69 G4RDeIonisationSpectrum::~G4RDeIonisationSpectrum() 
 70 {
 71   delete theParam;
 72 }
 73 
 74 
 75 G4double G4RDeIonisationSpectrum::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 = (G4RDAtomicTransitionManager::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 << "G4RDeIonisationSpectrum::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 
161 G4double G4RDeIonisationSpectrum::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 = (G4RDAtomicTransitionManager::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 << "G4RDeIonisationSpectrum::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 
243 G4double G4RDeIonisationSpectrum::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 = (G4RDAtomicTransitionManager::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 << "G4RDeIonisationSpectrum::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 G4RDeIonisationSpectrum::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 G4RDeIonisationSpectrum::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 
401 G4double G4RDeIonisationSpectrum::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 
485 G4double G4RDeIonisationSpectrum::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 
564 void G4RDeIonisationSpectrum::PrintData() const 
565 {
566   theParam->PrintData();
567 }
568 
569 G4double G4RDeIonisationSpectrum::MaxEnergyOfSecondaries(G4double kineticEnergy,
570                    G4int, // Z = 0,
571                    const G4ParticleDefinition* ) const
572 {
573   return 0.5 * kineticEnergy;
574 }
575