| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26
27 //....oooOO0OOooo........oooOO0OOooo........oo
28
29 #include <cmath>
30 #include <iostream>
31 #include "G4ecpssrBaseKxsModel.hh"
32 #include "globals.hh"
33 #include "G4PhysicalConstants.hh"
34 #include "G4SystemOfUnits.hh"
35 #include "G4AtomicTransitionManager.hh"
36 #include "G4NistManager.hh"
37 #include "G4Proton.hh"
38 #include "G4Alpha.hh"
39 #include "G4SemiLogInterpolation.hh"
40 #include "G4Exp.hh"
41
42 //....oooOO0OOooo........oooOO0OOooo........oo
43
44 G4ecpssrBaseKxsModel::G4ecpssrBaseKxsModel()
45 {
46 verboseLevel=0;
47
48 // Storing C coefficients for high velocit
49 G4String fileC1("pixe/uf/c1");
50 tableC1 = new G4CrossSectionDataSet(new G4
51
52 G4String fileC2("pixe/uf/c2");
53 tableC2 = new G4CrossSectionDataSet(new G4
54
55 G4String fileC3("pixe/uf/c3");
56 tableC3 = new G4CrossSectionDataSet(new G4
57
58 // Storing FK data needed for medium veloc
59 const char* path = G4FindDataDir("G4LEDATA
60
61 if (!path) {
62 G4Exception("G4ecpssrBaseKxsModel::G4ecp
63 return;
64 }
65
66 std::ostringstream fileName;
67 fileName << path << "/pixe/uf/FK.dat";
68 std::ifstream FK(fileName.str().c_str());
69
70 if (!FK)
71 G4Exception("G4ecpssrBaseKxsModel::G4ecp
72
73 dummyVec.push_back(0.);
74
75 while(!FK.eof())
76 {
77 double x;
78 double y;
79
80 FK>>x>>y;
81
82 // Mandatory vector initialization
83 if (x != dummyVec.back())
84 {
85 dummyVec.push_back(x);
86 aVecMap[x].push_back(-1.);
87 }
88
89 FK>>FKData[x][y];
90
91 if (y != aVecMap[x].back()) aVecMap[x]
92
93 }
94
95 tableC1->LoadData(fileC1);
96 tableC2->LoadData(fileC2);
97 tableC3->LoadData(fileC3);
98
99 }
100
101 //....oooOO0OOooo........oooOO0OOooo........oo
102
103 void print (G4double elem)
104 {
105 G4cout << elem << " ";
106 }
107 //....oooOO0OOooo........oooOO0OOooo........oo
108
109 G4ecpssrBaseKxsModel::~G4ecpssrBaseKxsModel()
110 {
111 delete tableC1;
112 delete tableC2;
113 delete tableC3;
114 }
115
116 //....oooOO0OOooo........oooOO0OOooo........oo
117
118 G4double G4ecpssrBaseKxsModel::ExpIntFunction(
119
120 {
121 // this "ExpIntFunction" function allows fast
122 G4int i;
123 G4int ii;
124 G4int nm1;
125 G4double a;
126 G4double b;
127 G4double c;
128 G4double d;
129 G4double del;
130 G4double fact;
131 G4double h;
132 G4double psi;
133 G4double ans = 0;
134 const G4double euler= 0.5772156649;
135 const G4int maxit= 100;
136 const G4double fpmin = 1.0e-30;
137 const G4double eps = 1.0e-7;
138 nm1=n-1;
139 if (n<0 || x<0.0 || (x==0.0 && (n==0 || n==1
140 G4cout << "*** WARNING in G4ecpssrBaseKxsM
141 G4cout << n << ", " << x << G4endl;
142 }
143 else {
144 if (n==0) ans=G4Exp(-x)/x;
145 else {
146 if (x==0.0) ans=1.0/nm1;
147 else {
148 if (x > 1.0) {
149 b=x+n;
150 c=1.0/fpmin;
151 d=1.0/b;
152 h=d;
153 for (i=1;i<=maxit;i++) {
154 a=-i*(nm1+i);
155 b +=2.0;
156 d=1.0/(a*d+b);
157 c=b+a/c;
158 del=c*d;
159 h *=del;
160 if (std::fabs(del-1.0) < eps) {
161 ans=h*G4Exp(-x);
162 return ans;
163 }
164 }
165 } else {
166 ans = (nm1!=0 ? 1.0/nm1 : -std::log(x
167 fact=1.0;
168 for (i=1;i<=maxit;i++) {
169 fact *=-x/i;
170 if (i !=nm1) del = -fact/(i-nm1);
171 else {
172 psi = -euler;
173 for (ii=1;ii<=nm1;ii++) psi +=1.0/ii;
174 del=fact*(-std::log(x)+psi);
175 }
176 ans += del;
177 if (std::fabs(del) < std::fabs(ans)
178 }
179 }
180 }
181 }
182 }
183 return ans;
184 }
185
186 //....oooOO0OOooo........oooOO0OOooo........oo
187
188 G4double G4ecpssrBaseKxsModel::CalculateCrossS
189
190 {
191 // this K-CrossSection calculation method is
192 G4NistManager* massManager = G4NistManager::
193
194 G4AtomicTransitionManager* transitionManage
195
196 G4double zIncident = 0;
197 G4Proton* aProtone = G4Proton::Proton();
198 G4Alpha* aAlpha = G4Alpha::Alpha();
199
200 if (massIncident == aProtone->GetPDGMass() )
201 {
202 zIncident = (aProtone->GetPDGCharge())/eplu
203 }
204 else
205 {
206 if (massIncident == aAlpha->GetPDGMass()
207 {
208 zIncident = (aAlpha->GetPDGCharge())/eplu
209 }
210 else
211 {
212 G4cout << "*** WARNING in G4ecpssrBaseKxsM
213 return 0;
214 }
215 }
216
217 if (verboseLevel>0) G4cout << " massIncid
218
219 G4double kBindingEnergy = transitionManager-
220
221 if (verboseLevel>0) G4cout << " kBindingE
222
223 G4double massTarget = (massManager->GetAtomi
224
225 if (verboseLevel>0) G4cout << " massTarge
226
227 G4double systemMass =((massIncident*massTarg
228
229 if (verboseLevel>0) G4cout << " systemMas
230
231 constexpr G4double zkshell= 0.3;
232 // *** see Brandt, Phys Rev A23, p 1727
233
234 G4double screenedzTarget = zTarget-zkshell;
235 // *** see Brandt, Phys Rev A23, p 1727
236
237 constexpr G4double rydbergMeV= 13.6056923e-6
238
239 G4double tetaK = kBindingEnergy/((screenedzT
240 // *** see Rice, ADANDT 20, p 504, f 2
241
242 if (verboseLevel>0) G4cout << " tetaK=" <
243
244 G4double velocity =(2./(tetaK*screenedzTarge
245 // *** also called xiK
246 // *** see Brandt, Phys Rev A23, p 1727
247 // *** see Basbas, Phys Rev A17, p 1656, f4
248
249 if (verboseLevel>0) G4cout << " velocity=
250
251 const G4double bohrPow2Barn=(Bohr_radius*Boh
252
253 if (verboseLevel>0) G4cout << " bohrPow2B
254
255 G4double sigma0 = 8.*pi*(zIncident*zIncident
256 // *** see Benka, ADANDT 22, p 220, f2, for
257 // *** see Basbas, Phys Rev A7, p 1000
258
259 if (verboseLevel>0) G4cout << " sigma0=" <<
260
261 const G4double kAnalyticalApproximation= 1.5
262 G4double x = kAnalyticalApproximation/veloci
263 // *** see Brandt, Phys Rev A23, p 1727
264 // *** see Brandt, Phys Rev A20, p 469, f16
265
266 if (verboseLevel>0) G4cout << " x=" << x<
267
268 G4double electrIonizationEnergy;
269 // *** see Basbas, Phys Rev A17, p1665, f27
270 // *** see Brandt, Phys Rev A20, p469
271 // *** see Liu, Comp Phys Comm 97, p325, f A
272
273 if ((0.< x) && (x <= 0.035))
274 {
275 electrIonizationEnergy= 0.75*pi*(std::lo
276 }
277 else
278 {
279 if ( (0.035 < x) && (x <=3.))
280 {
281 electrIonizationEnergy =G4Exp(-2.*x)/(0.03
282 }
283
284 else
285 {
286 if ( (3.< x) && (x<=11.))
287 {
288 electrIonizationEnergy =2.*G4Exp(-2.*x)
289 }
290
291 else electrIonizationEnergy =0.;
292 }
293 }
294
295 if (verboseLevel>0) G4cout << " electrIon
296
297 G4double hFunction =(electrIonizationEnergy*
298 // *** see Brandt, Phys Rev A20, p 469, f16
299
300 if (verboseLevel>0) G4cout << " hFunction
301
302 G4double gFunction = (1.+(9.*velocity)+(31.*
303 +(4.2*std::pow(velocity,6.))+(0.515*std:
304 // *** see Brandt, Phys Rev A20, p 469, f19
305
306 if (verboseLevel>0) G4cout << " gFunction
307
308 //------------------------------------------
309
310 G4double sigmaPSS = 1.+(((2.*zIncident)/(scr
311 // *** also called dzeta
312 // *** also called epsilon
313 // *** see Basbas, Phys Rev A17, p1667, f45
314
315 if (verboseLevel>0) G4cout << " sigmaPSS=
316
317 if (verboseLevel>0) G4cout << " sigmaPSS*
318
319 //------------------------------------------
320
321 const G4double cNaturalUnit= 1/fine_structur
322
323 if (verboseLevel>0) G4cout << " cNaturalU
324
325 G4double ykFormula=0.4*(screenedzTarget/cNat
326 // *** also called yS
327 // *** see Brandt, Phys Rev A20, p467, f6
328 // *** see Brandt, Phys Rev A23, p1728
329
330 if (verboseLevel>0) G4cout << " ykFormula
331
332 G4double relativityCorrection = std::pow((1.
333 // *** also called mRS
334 // *** see Brandt, Phys Rev A20, p467, f6
335
336 if (verboseLevel>0) G4cout << " relativit
337
338 G4double reducedVelocity = velocity*std::pow
339 // *** also called xiR
340 // *** see Brandt, Phys Rev A20, p468, f7
341 // *** see Brandt, Phys Rev A23, p1728
342
343 if (verboseLevel>0) G4cout << " reducedVe
344
345 G4double etaOverTheta2 = (energyIncident*ele
346 /(sigmaPSS*tetaK)/(
347 // *** see Benka, ADANDT 22, p220, f4 for et
348 // then we use sigmaPSS*tetaK == epsilon*tet
349
350 if (verboseLevel>0) G4cout << " etaOverTh
351
352 G4double universalFunction = 0;
353
354 // low velocity formula
355 // *****************
356 if ( velocity < 1. )
357 // OR
358 //if ( reducedVelocity/sigmaPSS < 1.)
359 // *** see Brandt, Phys Rev A23, p1727
360 // *** reducedVelocity/sigmaPSS is also call
361 // *****************
362 {
363 if (verboseLevel>0) G4cout << " Notice
364
365 universalFunction = (std::pow(2.,9.)/45.
366 // *** see Brandt, Phys Rev A23, p1728
367
368 if (verboseLevel>0) G4cout << " univers
369
370 }
371 else
372 {
373
374 if ( etaOverTheta2 > 86.6 && (sigmaPSS*tet
375 {
376 // High and medium energies. Method from
377
378 if (verboseLevel>0) G4cout << " Notice
379
380 if (verboseLevel>0) G4cout << " sigmaPS
381
382 G4double C1= tableC1->FindValue(sigmaPSS
383 G4double C2= tableC2->FindValue(sigmaPSS
384 G4double C3= tableC3->FindValue(sigmaPSS
385
386 if (verboseLevel>0) G4cout << " C1="
387 if (verboseLevel>0) G4cout << " C2="
388 if (verboseLevel>0) G4cout << " C3="
389
390 G4double etaK = (energyIncident*electron
391 // *** see Benka, ADANDT 22, p220, f4 fo
392
393 if (verboseLevel>0) G4cout << " etaK=
394
395 G4double etaT = (sigmaPSS*tetaK)*(sigmaP
396 // *** see Rice, ADANDT 20, p506
397
398 if (verboseLevel>0) G4cout << " etaT=
399
400 G4double fKT = FunctionFK((sigmaPSS*teta
401 // *** see Rice, ADANDT 20, p506
402
403 if (FunctionFK((sigmaPSS*tetaK),86.6)<=0
404 {
405 G4cout <<
406 "*** WARNING in G4ecpssrBaseKxsModel::
407 return 0;
408 }
409
410 if (verboseLevel>0) G4cout << " Funct
411
412 if (verboseLevel>0) G4cout << " fKT="
413
414 G4double GK = C2/(4*etaK) + C3/(32*etaK*
415
416 if (verboseLevel>0) G4cout << " GK=" << GK
417
418 G4double GT = C2/(4*etaT) + C3/(32*etaT*
419
420 if (verboseLevel>0) G4cout << " GT="
421
422 G4double DT = fKT - C1*std::log(etaT) +
423
424 if (verboseLevel>0) G4cout << " DT="
425
426 G4double fKK = C1*std::log(etaK) + DT -
427
428 if (verboseLevel>0) G4cout << " fKK="
429
430 G4double universalFunction3= fKK/(etaK/t
431 // *** see Rice, ADANDT 20, p505, f7
432
433 if (verboseLevel>0) G4cout << " unive
434
435 universalFunction=universalFunction3;
436
437 }
438 else if ( etaOverTheta2 >= 1.e-3 && etaOve
439 {
440 // From Benka 1978
441
442 if (verboseLevel>0) G4cout << " Notice
443
444 G4double universalFunction2 = FunctionFK
445
446 if (universalFunction2<=0)
447 {
448 G4cout <<
449 "*** WARNING : G4ecpssrBaseKxsModel::C
450 return 0;
451 }
452
453 if (verboseLevel>0) G4cout << " univers
454
455 universalFunction=universalFunction2;
456 }
457
458 }
459
460 //------------------------------------------
461
462 G4double sigmaPSSR = (sigma0/(sigmaPSS*tetaK
463 // *** see Benka, ADANDT 22, p220, f1
464
465 if (verboseLevel>0) G4cout << " sigmaPSSR
466
467 //------------------------------------------
468
469 G4double pssDeltaK = (4./(systemMass*sigmaPS
470 // *** also called dzetaK*deltaK
471 // *** see Brandt, Phys Rev A23, p1727, f B2
472
473 if (verboseLevel>0) G4cout << " pssDeltaK
474
475 if (pssDeltaK>1) return 0.;
476
477 G4double energyLoss = std::pow(1-pssDeltaK,0
478 // *** also called zK
479 // *** see Brandt, Phys Rev A23, p1727, afte
480
481 if (verboseLevel>0) G4cout << " energyLos
482
483 G4double energyLossFunction = (std::pow(2.,-
484 // *** also called fs
485 // *** see Brandt, Phys Rev A23, p1718, f7
486
487 if (verboseLevel>0) G4cout << " energyLos
488
489 //------------------------------------------
490
491 G4double coulombDeflection = (4.*pi*zInciden
492 // *** see Brandt, Phys Rev A23, p1727, f B3
493
494 if (verboseLevel>0) G4cout << " cParamete
495
496 G4double cParameter = 2.*coulombDeflection/(
497 // *** see Brandt, Phys Rev A23, p1727, f B4
498
499 if (verboseLevel>0) G4cout << " cParamete
500
501 G4double coulombDeflectionFunction = 9.*ExpI
502 // *** see Brandt, Phys Rev A23, p1727
503
504 if (verboseLevel>0) G4cout << " ExpIntFun
505
506 if (verboseLevel>0) G4cout << " coulombDe
507
508 //------------------------------------------
509
510 G4double crossSection = 0;
511
512 crossSection = energyLossFunction* coulombDe
513
514
515
516 //------------------------------------------
517
518 if (crossSection >= 0) {
519 return crossSection * barn;
520 }
521 else {return 0;}
522
523 }
524
525 //....oooOO0OOooo........oooOO0OOooo........oo
526
527 G4double G4ecpssrBaseKxsModel::FunctionFK(G4do
528 {
529
530 G4double sigma = 0.;
531 G4double valueT1 = 0;
532 G4double valueT2 = 0;
533 G4double valueE21 = 0;
534 G4double valueE22 = 0;
535 G4double valueE12 = 0;
536 G4double valueE11 = 0;
537 G4double xs11 = 0;
538 G4double xs12 = 0;
539 G4double xs21 = 0;
540 G4double xs22 = 0;
541
542 // PROTECTION TO ALLOW INTERPOLATION AT MINI
543 // (in particular for FK computation at 8.66
544
545 if (
546 theta==8.66e-3 ||
547 theta==8.66e-2 ||
548 theta==8.66e-1 ||
549 theta==8.66e+0 ||
550 theta==8.66e+1
551 ) theta=theta-1e-12;
552
553 if (
554 theta==1.e-3 ||
555 theta==1.e-2 ||
556 theta==1.e-1 ||
557 theta==1.e+00 ||
558 theta==1.e+01
559 ) theta=theta+1e-12;
560
561 // END PROTECTION
562
563 auto t2 = std::upper_bound(dummyVec.begin(),
564 auto t1 = t2-1;
565
566 auto e12 = std::upper_bound(aVecMap[(*t1)].b
567 auto e11 = e12-1;
568
569 auto e22 = std::upper_bound(aVecMap[(*t2)].b
570 auto e21 = e22-1;
571
572 valueT1 =*t1;
573 valueT2 =*t2;
574 valueE21 =*e21;
575 valueE22 =*e22;
576 valueE12 =*e12;
577 valueE11 =*e11;
578
579 xs11 = FKData[valueT1][valueE11];
580 xs12 = FKData[valueT1][valueE12];
581 xs21 = FKData[valueT2][valueE21];
582 xs22 = FKData[valueT2][valueE22];
583
584 G4double xsProduct = xs11 * xs12 * xs21 * xs
585
586 if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0)
587
588 if (xsProduct != 0.)
589 {
590 sigma = QuadInterpolator( valueE11, value
591 valueE21, valueE22,
592 xs11, xs12,
593 xs21, xs22,
594 valueT1, valueT2,
595 k, theta );
596 }
597
598 return sigma;
599 }
600
601 //....oooOO0OOooo........oooOO0OOooo........oo
602
603 G4double G4ecpssrBaseKxsModel::LinLogInterpola
604 G4double e2,
605 G4double e,
606 G4double xs1,
607 G4double xs2)
608 {
609 G4double d1 = std::log(xs1);
610 G4double d2 = std::log(xs2);
611 G4double value = G4Exp(d1 + (d2 - d1)*(e - e
612 return value;
613 }
614
615 //....oooOO0OOooo........oooOO0OOooo........oo
616
617 G4double G4ecpssrBaseKxsModel::LogLogInterpola
618 G4double e2,
619 G4double e,
620 G4double xs1,
621 G4double xs2)
622 {
623 G4double a = (std::log10(xs2)-std::log10(xs1
624 G4double b = std::log10(xs2) - a*std::log10(
625 G4double sigma = a*std::log10(e) + b;
626 G4double value = (std::pow(10.,sigma));
627 return value;
628 }
629
630 //....oooOO0OOooo........oooOO0OOooo........oo
631
632 G4double G4ecpssrBaseKxsModel::QuadInterpolato
633 G4double e21, G4double e22,
634 G4double xs11, G4double xs1
635 G4double xs21, G4double xs2
636 G4double t1, G4double t2,
637 G4double t, G4double e)
638 {
639 // Log-Log
640 G4double interpolatedvalue1 = LogLogInterpol
641 G4double interpolatedvalue2 = LogLogInterpol
642 G4double value = LogLogInterpolate(t1, t2, t
643
644 return value;
645 }
646