Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4ecpssrBaseLixsModel.cc
Go to the documentation of this file.
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#include <cmath>
28#include <iostream>
29
31
32#include "globals.hh"
34#include "G4SystemOfUnits.hh"
36#include "G4NistManager.hh"
37#include "G4Proton.hh"
38#include "G4Alpha.hh"
40
41//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
42
44{
45 verboseLevel=0;
46
47 // Storing FLi data needed for 0.2 to 3.0 velocities region
48
49 char *path = getenv("G4LEDATA");
50
51 if (!path) {
52 G4Exception("G4ecpssrLCrossSection::G4ecpssrBaseLixsModel()","em0006", FatalException ,"G4LEDATA environment variable not set");
53 return;
54 }
55 std::ostringstream fileName1;
56 std::ostringstream fileName2;
57
58 fileName1 << path << "/pixe/uf/FL1.dat";
59 fileName2 << path << "/pixe/uf/FL2.dat";
60
61 // Reading of FL1.dat
62
63 std::ifstream FL1(fileName1.str().c_str());
64 if (!FL1) G4Exception("G4ecpssrLCrossSection::G4ecpssrBaseLixsModel()","em0003",FatalException, "error opening FL1 data file");
65
66 dummyVec1.push_back(0.);
67
68 while(!FL1.eof())
69 {
70 double x1;
71 double y1;
72
73 FL1>>x1>>y1;
74
75 // Mandatory vector initialization
76 if (x1 != dummyVec1.back())
77 {
78 dummyVec1.push_back(x1);
79 aVecMap1[x1].push_back(-1.);
80 }
81
82 FL1>>FL1Data[x1][y1];
83
84 if (y1 != aVecMap1[x1].back()) aVecMap1[x1].push_back(y1);
85 }
86
87 // Reading of FL2.dat
88
89 std::ifstream FL2(fileName2.str().c_str());
90 if (!FL2) G4Exception("G4ecpssrLCrossSection::G4ecpssrBaseLixsModel()","em0003", FatalException," error opening FL2 data file");
91
92 dummyVec2.push_back(0.);
93
94 while(!FL2.eof())
95 {
96 double x2;
97 double y2;
98
99 FL2>>x2>>y2;
100
101 // Mandatory vector initialization
102 if (x2 != dummyVec2.back())
103 {
104 dummyVec2.push_back(x2);
105 aVecMap2[x2].push_back(-1.);
106 }
107
108 FL2>>FL2Data[x2][y2];
109
110 if (y2 != aVecMap2[x2].back()) aVecMap2[x2].push_back(y2);
111 }
112
113}
114
115//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
116
118{}
119
120//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
121
123
124{
125// this function allows fast evaluation of the n order exponential integral function En(x)
126
127 G4int i;
128 G4int ii;
129 G4int nm1;
130 G4double a;
131 G4double b;
132 G4double c;
133 G4double d;
134 G4double del;
135 G4double fact;
136 G4double h;
137 G4double psi;
138 G4double ans = 0;
139 const G4double euler= 0.5772156649;
140 const G4int maxit= 100;
141 const G4double fpmin = 1.0e-30;
142 const G4double eps = 1.0e-7;
143 nm1=n-1;
144 if (n<0 || x<0.0 || (x==0.0 && (n==0 || n==1)))
145 G4cout << "*** WARNING in G4ecpssrBaseLixsModel::ExpIntFunction: bad arguments in ExpIntFunction"
146 << G4endl;
147 else {
148 if (n==0) ans=std::exp(-x)/x;
149 else {
150 if (x==0.0) ans=1.0/nm1;
151 else {
152 if (x > 1.0) {
153 b=x+n;
154 c=1.0/fpmin;
155 d=1.0/b;
156 h=d;
157 for (i=1;i<=maxit;i++) {
158 a=-i*(nm1+i);
159 b +=2.0;
160 d=1.0/(a*d+b);
161 c=b+a/c;
162 del=c*d;
163 h *=del;
164 if (std::fabs(del-1.0) < eps) {
165 ans=h*std::exp(-x);
166 return ans;
167 }
168 }
169 } else {
170 ans = (nm1!=0 ? 1.0/nm1 : -std::log(x)-euler);
171 fact=1.0;
172 for (i=1;i<=maxit;i++) {
173 fact *=-x/i;
174 if (i !=nm1) del = -fact/(i-nm1);
175 else {
176 psi = -euler;
177 for (ii=1;ii<=nm1;ii++) psi +=1.0/ii;
178 del=fact*(-std::log(x)+psi);
179 }
180 ans += del;
181 if (std::fabs(del) < std::fabs(ans)*eps) return ans;
182 }
183 }
184 }
185 }
186 }
187return ans;
188}
189
190//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
191
193{
194
195 if (zTarget <=4) return 0.;
196
197 //this L1-CrossSection calculation method is done according to Werner Brandt and Grzegorz Lapicki, Phys.Rev.A20 N2 (1979),
198 //and using data tables of O. Benka et al. At.Data Nucl.Data Tables Vol.22 No.3 (1978).
199
200 G4NistManager* massManager = G4NistManager::Instance();
201
203
204 G4double zIncident = 0;
205 G4Proton* aProtone = G4Proton::Proton();
206 G4Alpha* aAlpha = G4Alpha::Alpha();
207
208 if (massIncident == aProtone->GetPDGMass() )
209
210 zIncident = (aProtone->GetPDGCharge())/eplus;
211
212 else
213 {
214 if (massIncident == aAlpha->GetPDGMass())
215
216 zIncident = (aAlpha->GetPDGCharge())/eplus;
217
218 else
219 {
220 G4cout << "*** WARNING in G4ecpssrBaseLixsModel::CalculateL1CrossSection : Proton or Alpha incident particles only. " << G4endl;
221 G4cout << massIncident << ", " << aAlpha->GetPDGMass() << " (alpha)" << aProtone->GetPDGMass() << " (proton)" << G4endl;
222 return 0;
223 }
224 }
225
226 G4double l1BindingEnergy = transitionManager->Shell(zTarget,1)->BindingEnergy(); //Observed binding energy of L1-subshell
227
228 G4double massTarget = (massManager->GetAtomicMassAmu(zTarget))*amu_c2;
229
230 G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2; //Mass of the system (projectile, target)
231
232 const G4double zlshell= 4.15;
233 // *** see Benka, ADANDT 22, p 223
234
235 G4double screenedzTarget = zTarget-zlshell; //Effective nuclear charge as seen by electrons in L1-sub shell
236
237 const G4double rydbergMeV= 13.6056923e-6;
238
239 const G4double nl= 2.;
240 // *** see Benka, ADANDT 22, p 220, f3
241
242 G4double tetal1 = (l1BindingEnergy*nl*nl)/((screenedzTarget*screenedzTarget)*rydbergMeV); //Screening parameter
243 // *** see Benka, ADANDT 22, p 220, f3
244
245 if (verboseLevel>0) G4cout << " tetal1=" << tetal1<< G4endl;
246
247 G4double reducedEnergy = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
248 // *** also called etaS
249 // *** see Benka, ADANDT 22, p 220, f3
250
251 const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/barn ; //Bohr radius of hydrogen
252
253 G4double sigma0 = 8.*pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.);
254 // *** see Benka, ADANDT 22, p 220, f2, for protons
255 // *** see Basbas, Phys Rev A7, p 1000
256
257 G4double velocityl1 = CalculateVelocity(1, zTarget, massIncident, energyIncident); // Scaled velocity
258
259 if (verboseLevel>0) G4cout << " velocityl1=" << velocityl1<< G4endl;
260
261 const G4double l1AnalyticalApproximation= 1.5;
262 G4double x1 =(nl*l1AnalyticalApproximation)/velocityl1;
263 // *** 1.5 is cK = cL1 (it is 1.25 for L2 & L3)
264 // *** see Brandt, Phys Rev A20, p 469, f16 in expression of h
265
266 if (verboseLevel>0) G4cout << " x1=" << x1<< G4endl;
267
268 G4double electrIonizationEnergyl1=0.;
269 // *** see Basbas, Phys Rev A17, p1665, f27
270 // *** see Brandt, Phys Rev A20, p469
271 // *** see Liu, Comp Phys Comm 97, p325, f A5
272
273 if ( x1<=0.035) electrIonizationEnergyl1= 0.75*pi*(std::log(1./(x1*x1))-1.);
274 else
275 {
276 if ( x1<=3.)
277 electrIonizationEnergyl1 =std::exp(-2.*x1)/(0.031+(0.213*std::pow(x1,0.5))+(0.005*x1)-(0.069*std::pow(x1,3./2.))+(0.324*x1*x1));
278 else
279 {if ( x1<=11.) electrIonizationEnergyl1 =2.*std::exp(-2.*x1)/std::pow(x1,1.6);}
280 }
281
282 G4double hFunctionl1 =(electrIonizationEnergyl1*2.*nl)/(tetal1*std::pow(velocityl1,3)); //takes into account the polarization effect
283 // *** see Brandt, Phys Rev A20, p 469, f16
284
285 if (verboseLevel>0) G4cout << " hFunctionl1=" << hFunctionl1<< G4endl;
286
287 G4double gFunctionl1 = (1.+(9.*velocityl1)+(31.*velocityl1*velocityl1)+(49.*std::pow(velocityl1,3.))+(162.*std::pow(velocityl1,4.))+(63.*std::pow(velocityl1,5.))+(18.*std::pow(velocityl1,6.))+(1.97*std::pow(velocityl1,7.)))/std::pow(1.+velocityl1,9.);//takes into account the reduced binding effect
288 // *** see Brandt, Phys Rev A20, p 469, f19
289
290 if (verboseLevel>0) G4cout << " gFunctionl1=" << gFunctionl1<< G4endl;
291
292 G4double sigmaPSS_l1 = 1.+(((2.*zIncident)/(screenedzTarget*tetal1))*(gFunctionl1-hFunctionl1)); //Binding-polarization factor
293 // *** also called dzeta
294 // *** also called epsilon
295 // *** see Basbas, Phys Rev A17, p1667, f45
296
297 if (verboseLevel>0) G4cout << "sigmaPSS_l1 =" << sigmaPSS_l1<< G4endl;
298
299 const G4double cNaturalUnit= 137.;
300
301 G4double yl1Formula=0.4*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(nl*velocityl1/sigmaPSS_l1);
302 // *** also called yS
303 // *** see Brandt, Phys Rev A20, p467, f6
304 // *** see Brandt, Phys Rev A23, p1728
305
306 G4double l1relativityCorrection = std::pow((1.+(1.1*yl1Formula*yl1Formula)),0.5)+yl1Formula; // Relativistic correction parameter
307 // *** also called mRS
308 // *** see Brandt, Phys Rev A20, p467, f6
309
310 //G4double reducedVelocity_l1 = velocityl1*std::pow(l1relativityCorrection,0.5); //Reduced velocity parameter
311
312 G4double L1etaOverTheta2;
313
314 G4double universalFunction_l1 = 0.;
315
316 G4double sigmaPSSR_l1;
317
318 // low velocity formula
319 // *****************
320 if ( velocityl1 <20. )
321 {
322
323 L1etaOverTheta2 =(reducedEnergy* l1relativityCorrection)/((tetal1*sigmaPSS_l1)*(tetal1*sigmaPSS_l1));
324 // *** 1) RELATIVISTIC CORRECTION ADDED
325 // *** 2) sigma_PSS_l1 ADDED
326 // *** reducedEnergy is etaS, l1relativityCorrection is mRS
327 // *** see Phys Rev A20, p468, top
328
329 if ( ((tetal1*sigmaPSS_l1) >=0.2) && ((tetal1*sigmaPSS_l1) <=2.6670) && (L1etaOverTheta2>=0.1e-3) && (L1etaOverTheta2<=0.866e2) )
330
331 universalFunction_l1 = FunctionFL1((tetal1*sigmaPSS_l1),L1etaOverTheta2);
332
333 if (verboseLevel>0) G4cout << "at low velocity range, universalFunction_l1 =" << universalFunction_l1 << G4endl;
334
335 sigmaPSSR_l1 = (sigma0/(tetal1*sigmaPSS_l1))*universalFunction_l1;// Plane-wave Born -Aproximation L1-subshell ionisation Cross Section
336 // *** see Benka, ADANDT 22, p220, f1
337
338 if (verboseLevel>0) G4cout << " at low velocity range, sigma PWBA L1 CS = " << sigmaPSSR_l1<< G4endl;
339
340 }
341
342 else
343
344 {
345
346 L1etaOverTheta2 = reducedEnergy/(tetal1*tetal1);
347 // Medium & high velocity
348 // *** 1) NO RELATIVISTIC CORRECTION
349 // *** 2) NO sigma_PSS_l1
350 // *** see Benka, ADANDT 22, p223
351
352 if ( (tetal1 >=0.2) && (tetal1 <=2.6670) && (L1etaOverTheta2>=0.1e-3) && (L1etaOverTheta2<=0.866e2) )
353
354 universalFunction_l1 = FunctionFL1(tetal1,L1etaOverTheta2);
355
356 if (verboseLevel>0) G4cout << "at medium and high velocity range, universalFunction_l1 =" << universalFunction_l1 << G4endl;
357
358 sigmaPSSR_l1 = (sigma0/tetal1)*universalFunction_l1;// Plane-wave Born -Aproximation L1-subshell ionisation Cross Section
359 // *** see Benka, ADANDT 22, p220, f1
360
361 if (verboseLevel>0) G4cout << " sigma PWBA L1 CS at medium and high velocity range = " << sigmaPSSR_l1<< G4endl;
362 }
363
364 G4double pssDeltal1 = (4./(systemMass *sigmaPSS_l1*tetal1))*(sigmaPSS_l1/velocityl1)*(sigmaPSS_l1/velocityl1);
365 // *** also called dzeta*delta
366 // *** see Brandt, Phys Rev A23, p1727, f B2
367
368 if (verboseLevel>0) G4cout << " pssDeltal1=" << pssDeltal1<< G4endl;
369
370 if (pssDeltal1>1) return 0.;
371
372 G4double energyLossl1 = std::pow(1-pssDeltal1,0.5);
373 // *** also called z
374 // *** see Brandt, Phys Rev A23, p1727, after f B2
375
376 if (verboseLevel>0) G4cout << " energyLossl1=" << energyLossl1<< G4endl;
377
378 G4double coulombDeflectionl1 =
379 (8.*pi*zIncident/systemMass)*std::pow(tetal1*sigmaPSS_l1,-2.)*std::pow(velocityl1/sigmaPSS_l1,-3.)*(zTarget/screenedzTarget);
380 // *** see Brandt, Phys Rev A20, v2s and f2 and B2
381 // *** with factor n2 compared to Brandt, Phys Rev A23, p1727, f B3
382
383 G4double cParameterl1 =2.* coulombDeflectionl1/(energyLossl1*(energyLossl1+1.));
384 // *** see Brandt, Phys Rev A23, p1727, f B4
385
386 G4double coulombDeflectionFunction_l1 = 9.*ExpIntFunction(10,cParameterl1); //Coulomb-deflection effect correction
387
388 if (verboseLevel>0) G4cout << " coulombDeflectionFunction_l1 =" << coulombDeflectionFunction_l1 << G4endl;
389
390 G4double crossSection_L1 = coulombDeflectionFunction_l1 * sigmaPSSR_l1;
391
392 //ECPSSR L1 -subshell cross section is estimated at perturbed-stationnairy-state(PSS)
393 //and reduced by the energy-loss(E),the Coulomb deflection(C),and the relativity(R) effects
394
395 if (verboseLevel>0) G4cout << " crossSection_L1 =" << crossSection_L1 << G4endl;
396
397 if (crossSection_L1 >= 0)
398 {
399 return crossSection_L1 * barn;
400 }
401
402 else {return 0;}
403}
404
405//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
406
408
409{
410 if (zTarget <=13 ) return 0.;
411
412 // this L2-CrossSection calculation method is done according to Werner Brandt and Grzegorz Lapicki, Phys.Rev.A20 N2 (1979),
413 // and using data tables of O. Benka et al. At.Data Nucl.Data Tables Vol.22 No.3 (1978).
414
415 G4NistManager* massManager = G4NistManager::Instance();
416
418
419 G4double zIncident = 0;
420
421 G4Proton* aProtone = G4Proton::Proton();
422 G4Alpha* aAlpha = G4Alpha::Alpha();
423
424 if (massIncident == aProtone->GetPDGMass() )
425
426 zIncident = (aProtone->GetPDGCharge())/eplus;
427
428 else
429 {
430 if (massIncident == aAlpha->GetPDGMass())
431
432 zIncident = (aAlpha->GetPDGCharge())/eplus;
433
434 else
435 {
436 G4cout << "*** WARNING in G4ecpssrBaseLixsModel::CalculateL2CrossSection : Proton or Alpha incident particles only. " << G4endl;
437 G4cout << massIncident << ", " << aAlpha->GetPDGMass() << " (alpha)" << aProtone->GetPDGMass() << " (proton)" << G4endl;
438 return 0;
439 }
440 }
441
442 G4double l2BindingEnergy = transitionManager->Shell(zTarget,2)->BindingEnergy(); //Observed binding energy of L2-subshell
443
444 G4double massTarget = (massManager->GetAtomicMassAmu(zTarget))*amu_c2;
445
446 G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2; //Mass of the system (projectile, target)
447
448 const G4double zlshell= 4.15;
449
450 G4double screenedzTarget = zTarget-zlshell; //Effective nuclear charge as seen by electrons in L2-subshell
451
452 const G4double rydbergMeV= 13.6056923e-6;
453
454 const G4double nl= 2.;
455
456 G4double tetal2 = (l2BindingEnergy*nl*nl)/((screenedzTarget*screenedzTarget)*rydbergMeV); //Screening parameter
457
458 if (verboseLevel>0) G4cout << " tetal2=" << tetal2<< G4endl;
459
460 G4double reducedEnergy = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
461
462 const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/barn ; //Bohr radius of hydrogen
463
464 G4double sigma0 = 8.*pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.);
465
466 G4double velocityl2 = CalculateVelocity(2, zTarget, massIncident, energyIncident); // Scaled velocity
467
468 if (verboseLevel>0) G4cout << " velocityl2=" << velocityl2<< G4endl;
469
470 const G4double l23AnalyticalApproximation= 1.25;
471
472 G4double x2 = (nl*l23AnalyticalApproximation)/velocityl2;
473
474 if (verboseLevel>0) G4cout << " x2=" << x2<< G4endl;
475
476 G4double electrIonizationEnergyl2=0.;
477
478 if ( x2<=0.035) electrIonizationEnergyl2= 0.75*pi*(std::log(1./(x2*x2))-1.);
479 else
480 {
481 if ( x2<=3.)
482 electrIonizationEnergyl2 =std::exp(-2.*x2)/(0.031+(0.213*std::pow(x2,0.5))+(0.005*x2)-(0.069*std::pow(x2,3./2.))+(0.324*x2*x2));
483 else
484 {if ( x2<=11.) electrIonizationEnergyl2 =2.*std::exp(-2.*x2)/std::pow(x2,1.6); }
485 }
486
487 G4double hFunctionl2 =(electrIonizationEnergyl2*2.*nl)/(tetal2*std::pow(velocityl2,3)); //takes into account the polarization effect
488
489 if (verboseLevel>0) G4cout << " hFunctionl2=" << hFunctionl2<< G4endl;
490
491 G4double gFunctionl2 = (1.+(10.*velocityl2)+(45.*velocityl2*velocityl2)+(102.*std::pow(velocityl2,3.))+(331.*std::pow(velocityl2,4.))+(6.7*std::pow(velocityl2,5.))+(58.*std::pow(velocityl2,6.))+(7.8*std::pow(velocityl2,7.))+ (0.888*std::pow(velocityl2,8.)) )/std::pow(1.+velocityl2,10.);
492 //takes into account the reduced binding effect
493
494 if (verboseLevel>0) G4cout << " gFunctionl2=" << gFunctionl2<< G4endl;
495
496 G4double sigmaPSS_l2 = 1.+(((2.*zIncident)/(screenedzTarget*tetal2))*(gFunctionl2-hFunctionl2)); //Binding-polarization factor
497
498 if (verboseLevel>0) G4cout << " sigmaPSS_l2=" << sigmaPSS_l2<< G4endl;
499
500 const G4double cNaturalUnit= 137.;
501
502 G4double yl2Formula=0.15*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(velocityl2/sigmaPSS_l2);
503
504 G4double l2relativityCorrection = std::pow((1.+(1.1*yl2Formula*yl2Formula)),0.5)+yl2Formula; // Relativistic correction parameter
505
506 G4double L2etaOverTheta2;
507
508 G4double universalFunction_l2 = 0.;
509
510 G4double sigmaPSSR_l2 ;
511
512 if ( velocityl2 < 20. )
513 {
514
515 L2etaOverTheta2 = (reducedEnergy*l2relativityCorrection)/((sigmaPSS_l2*tetal2)*(sigmaPSS_l2*tetal2));
516
517 if ( (tetal2*sigmaPSS_l2>=0.2) && (tetal2*sigmaPSS_l2<=2.6670) && (L2etaOverTheta2>=0.1e-3) && (L2etaOverTheta2<=0.866e2) )
518
519 universalFunction_l2 = FunctionFL2((tetal2*sigmaPSS_l2),L2etaOverTheta2);
520
521 sigmaPSSR_l2 = (sigma0/(tetal2*sigmaPSS_l2))*universalFunction_l2;
522
523 if (verboseLevel>0) G4cout << " sigma PWBA L2 CS at low velocity range = " << sigmaPSSR_l2<< G4endl;
524 }
525 else
526 {
527
528 L2etaOverTheta2 = reducedEnergy /(tetal2*tetal2);
529
530 if ( (tetal2>=0.2) && (tetal2<=2.6670) && (L2etaOverTheta2>=0.1e-3) && (L2etaOverTheta2<=0.866e2) )
531
532 universalFunction_l2 = FunctionFL2((tetal2),L2etaOverTheta2);
533
534 sigmaPSSR_l2 = (sigma0/tetal2)*universalFunction_l2;
535
536 if (verboseLevel>0) G4cout << " sigma PWBA L2 CS at medium and high velocity range = " << sigmaPSSR_l2<< G4endl;
537
538 }
539
540 G4double pssDeltal2 = (4./(systemMass*sigmaPSS_l2*tetal2))*(sigmaPSS_l2/velocityl2)*(sigmaPSS_l2/velocityl2);
541
542 if (pssDeltal2>1) return 0.;
543
544 G4double energyLossl2 = std::pow(1-pssDeltal2,0.5);
545
546 if (verboseLevel>0) G4cout << " energyLossl2=" << energyLossl2<< G4endl;
547
548 G4double coulombDeflectionl2
549 =(8.*pi*zIncident/systemMass)*std::pow(tetal2*sigmaPSS_l2,-2.)*std::pow(velocityl2/sigmaPSS_l2,-3.)*(zTarget/screenedzTarget);
550
551 G4double cParameterl2 = 2.*coulombDeflectionl2/(energyLossl2*(energyLossl2+1.));
552
553 G4double coulombDeflectionFunction_l2 = 11.*ExpIntFunction(12,cParameterl2); //Coulomb-deflection effect correction
554 // *** see Brandt, Phys Rev A10, p477, f25
555
556 if (verboseLevel>0) G4cout << " coulombDeflectionFunction_l2 =" << coulombDeflectionFunction_l2 << G4endl;
557
558 G4double crossSection_L2 = coulombDeflectionFunction_l2 * sigmaPSSR_l2;
559 //ECPSSR L2 -subshell cross section is estimated at perturbed-stationnairy-state(PSS)
560 //and reduced by the energy-loss(E),the Coulomb deflection(C),and the relativity(R) effects
561
562 if (verboseLevel>0) G4cout << " crossSection_L2 =" << crossSection_L2 << G4endl;
563
564 if (crossSection_L2 >= 0)
565 {
566 return crossSection_L2 * barn;
567 }
568 else {return 0;}
569}
570
571//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
572
573
575
576{
577 if (zTarget <=13) return 0.;
578
579 //this L3-CrossSection calculation method is done according to Werner Brandt and Grzegorz Lapicki, Phys.Rev.A20 N2 (1979),
580 //and using data tables of O. Benka et al. At.Data Nucl.Data Tables Vol.22 No.3 (1978).
581
582 G4NistManager* massManager = G4NistManager::Instance();
583
585
586 G4double zIncident = 0;
587
588 G4Proton* aProtone = G4Proton::Proton();
589 G4Alpha* aAlpha = G4Alpha::Alpha();
590
591 if (massIncident == aProtone->GetPDGMass() )
592
593 zIncident = (aProtone->GetPDGCharge())/eplus;
594
595 else
596 {
597 if (massIncident == aAlpha->GetPDGMass())
598
599 zIncident = (aAlpha->GetPDGCharge())/eplus;
600
601 else
602 {
603 G4cout << "*** WARNING in G4ecpssrBaseLixsModel::CalculateL3CrossSection : Proton or Alpha incident particles only. " << G4endl;
604 G4cout << massIncident << ", " << aAlpha->GetPDGMass() << " (alpha)" << aProtone->GetPDGMass() << " (proton)" << G4endl;
605 return 0;
606 }
607 }
608
609 G4double l3BindingEnergy = transitionManager->Shell(zTarget,3)->BindingEnergy();
610
611 G4double massTarget = (massManager->GetAtomicMassAmu(zTarget))*amu_c2;
612
613 G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2;//Mass of the system (projectile, target)
614
615 const G4double zlshell= 4.15;
616
617 G4double screenedzTarget = zTarget-zlshell;//Effective nuclear charge as seen by electrons in L3-subshell
618
619 const G4double rydbergMeV= 13.6056923e-6;
620
621 const G4double nl= 2.;
622
623 G4double tetal3 = (l3BindingEnergy*nl*nl)/((screenedzTarget*screenedzTarget)*rydbergMeV);//Screening parameter
624
625 if (verboseLevel>0) G4cout << " tetal3=" << tetal3<< G4endl;
626
627 G4double reducedEnergy = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
628
629 const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/barn ;//Bohr radius of hydrogen
630
631 G4double sigma0 = 8.*pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.);
632
633 G4double velocityl3 = CalculateVelocity(3, zTarget, massIncident, energyIncident);// Scaled velocity
634
635 if (verboseLevel>0) G4cout << " velocityl3=" << velocityl3<< G4endl;
636
637 const G4double l23AnalyticalApproximation= 1.25;
638
639 G4double x3 = (nl*l23AnalyticalApproximation)/velocityl3;
640
641 if (verboseLevel>0) G4cout << " x3=" << x3<< G4endl;
642
643 G4double electrIonizationEnergyl3=0.;
644
645 if ( x3<=0.035) electrIonizationEnergyl3= 0.75*pi*(std::log(1./(x3*x3))-1.);
646 else
647 {
648 if ( x3<=3.) electrIonizationEnergyl3 =std::exp(-2.*x3)/(0.031+(0.213*std::pow(x3,0.5))+(0.005*x3)-(0.069*std::pow(x3,3./2.))+(0.324*x3*x3));
649 else
650 {
651 if ( x3<=11.) electrIonizationEnergyl3 =2.*std::exp(-2.*x3)/std::pow(x3,1.6);}
652 }
653
654 G4double hFunctionl3 =(electrIonizationEnergyl3*2.*nl)/(tetal3*std::pow(velocityl3,3));//takes into account the polarization effect
655
656 if (verboseLevel>0) G4cout << " hFunctionl3=" << hFunctionl3<< G4endl;
657
658 G4double gFunctionl3 = (1.+(10.*velocityl3)+(45.*velocityl3*velocityl3)+(102.*std::pow(velocityl3,3.))+(331.*std::pow(velocityl3,4.))+(6.7*std::pow(velocityl3,5.))+(58.*std::pow(velocityl3,6.))+(7.8*std::pow(velocityl3,7.))+ (0.888*std::pow(velocityl3,8.)) )/std::pow(1.+velocityl3,10.);
659 //takes into account the reduced binding effect
660
661 if (verboseLevel>0) G4cout << " gFunctionl3=" << gFunctionl3<< G4endl;
662
663 G4double sigmaPSS_l3 = 1.+(((2.*zIncident)/(screenedzTarget*tetal3))*(gFunctionl3-hFunctionl3));//Binding-polarization factor
664
665 if (verboseLevel>0) G4cout << "sigmaPSS_l3 =" << sigmaPSS_l3<< G4endl;
666
667 const G4double cNaturalUnit= 137.;
668
669 G4double yl3Formula=0.15*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(velocityl3/sigmaPSS_l3);
670
671 G4double l3relativityCorrection = std::pow((1.+(1.1*yl3Formula*yl3Formula)),0.5)+yl3Formula; // Relativistic correction parameter
672
673 G4double L3etaOverTheta2;
674
675 G4double universalFunction_l3 = 0.;
676
677 G4double sigmaPSSR_l3;
678
679 if ( velocityl3 < 20. )
680 {
681
682 L3etaOverTheta2 = (reducedEnergy* l3relativityCorrection)/((sigmaPSS_l3*tetal3)*(sigmaPSS_l3*tetal3));
683
684 if ( (tetal3*sigmaPSS_l3>=0.2) && (tetal3*sigmaPSS_l3<=2.6670) && (L3etaOverTheta2>=0.1e-3) && (L3etaOverTheta2<=0.866e2) )
685
686 universalFunction_l3 = 2.*FunctionFL2((tetal3*sigmaPSS_l3), L3etaOverTheta2 );
687
688 sigmaPSSR_l3 = (sigma0/(tetal3*sigmaPSS_l3))*universalFunction_l3;
689
690 if (verboseLevel>0) G4cout << " sigma PWBA L3 CS at low velocity range = " << sigmaPSSR_l3<< G4endl;
691
692 }
693
694 else
695
696 {
697
698 L3etaOverTheta2 = reducedEnergy/(tetal3*tetal3);
699
700 if ( (tetal3>=0.2) && (tetal3<=2.6670) && (L3etaOverTheta2>=0.1e-3) && (L3etaOverTheta2<=0.866e2) )
701
702 universalFunction_l3 = 2.*FunctionFL2(tetal3, L3etaOverTheta2 );
703
704 sigmaPSSR_l3 = (sigma0/tetal3)*universalFunction_l3;
705
706 if (verboseLevel>0) G4cout << " sigma PWBA L3 CS at medium and high velocity range = " << sigmaPSSR_l3<< G4endl;
707 }
708
709 G4double pssDeltal3 = (4./(systemMass*sigmaPSS_l3*tetal3))*(sigmaPSS_l3/velocityl3)*(sigmaPSS_l3/velocityl3);
710
711 if (verboseLevel>0) G4cout << " pssDeltal3=" << pssDeltal3<< G4endl;
712
713 if (pssDeltal3>1) return 0.;
714
715 G4double energyLossl3 = std::pow(1-pssDeltal3,0.5);
716
717 if (verboseLevel>0) G4cout << " energyLossl3=" << energyLossl3<< G4endl;
718
719 G4double coulombDeflectionl3 =
720 (8.*pi*zIncident/systemMass)*std::pow(tetal3*sigmaPSS_l3,-2.)*std::pow(velocityl3/sigmaPSS_l3,-3.)*(zTarget/screenedzTarget);
721
722 G4double cParameterl3 = 2.*coulombDeflectionl3/(energyLossl3*(energyLossl3+1.));
723
724 G4double coulombDeflectionFunction_l3 = 11.*ExpIntFunction(12,cParameterl3);//Coulomb-deflection effect correction
725 // *** see Brandt, Phys Rev A10, p477, f25
726
727 if (verboseLevel>0) G4cout << " coulombDeflectionFunction_l3 =" << coulombDeflectionFunction_l3 << G4endl;
728
729 G4double crossSection_L3 = coulombDeflectionFunction_l3 * sigmaPSSR_l3;
730 //ECPSSR L3 -subshell cross section is estimated at perturbed-stationnairy-state(PSS)
731 //and reduced by the energy-loss(E),the Coulomb deflection(C),and the relativity(R) effects
732
733 if (verboseLevel>0) G4cout << " crossSection_L3 =" << crossSection_L3 << G4endl;
734
735 if (crossSection_L3 >= 0)
736 {
737 return crossSection_L3 * barn;
738 }
739 else {return 0;}
740}
741
742//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
743
744G4double G4ecpssrBaseLixsModel::CalculateVelocity(G4int subShell, G4int zTarget, G4double massIncident, G4double energyIncident)
745
746{
747
749
750 G4double liBindingEnergy = transitionManager->Shell(zTarget,subShell)->BindingEnergy();
751
752 G4Proton* aProtone = G4Proton::Proton();
753 G4Alpha* aAlpha = G4Alpha::Alpha();
754
755 if (!((massIncident == aProtone->GetPDGMass()) || (massIncident == aAlpha->GetPDGMass())))
756 {
757 G4cout << "*** WARNING in G4ecpssrBaseLixsModel::CalculateVelocity : Proton or Alpha incident particles only. " << G4endl;
758 G4cout << massIncident << ", " << aAlpha->GetPDGMass() << " (alpha)" << aProtone->GetPDGMass() << " (proton)" << G4endl;
759 return 0;
760 }
761
762 const G4double zlshell= 4.15;
763
764 G4double screenedzTarget = zTarget- zlshell;
765
766 const G4double rydbergMeV= 13.6056923e-6;
767
768 const G4double nl= 2.;
769
770 G4double tetali = (liBindingEnergy*nl*nl)/(screenedzTarget*screenedzTarget*rydbergMeV);
771
772 G4double reducedEnergy = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
773
774 G4double velocity = 2.*nl*std::pow(reducedEnergy,0.5)/tetali;
775 // *** see Brandt, Phys Rev A10, p10, f4
776
777 return velocity;
778}
779
780//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
781
782G4double G4ecpssrBaseLixsModel::FunctionFL1(G4double k, G4double theta)
783{
784
785 G4double sigma = 0.;
786 G4double valueT1 = 0;
787 G4double valueT2 = 0;
788 G4double valueE21 = 0;
789 G4double valueE22 = 0;
790 G4double valueE12 = 0;
791 G4double valueE11 = 0;
792 G4double xs11 = 0;
793 G4double xs12 = 0;
794 G4double xs21 = 0;
795 G4double xs22 = 0;
796
797 // PROTECTION TO ALLOW INTERPOLATION AT MINIMUM AND MAXIMUM Eta/Theta2 values
798
799 if (
800 theta==8.66e-4 ||
801 theta==8.66e-3 ||
802 theta==8.66e-2 ||
803 theta==8.66e-1 ||
804 theta==8.66e+00 ||
805 theta==8.66e+01
806 ) theta=theta-1e-12;
807
808 if (
809 theta==1.e-4 ||
810 theta==1.e-3 ||
811 theta==1.e-2 ||
812 theta==1.e-1 ||
813 theta==1.e+00 ||
814 theta==1.e+01
815 ) theta=theta+1e-12;
816
817 // END PROTECTION
818
819 std::vector<double>::iterator t2 = std::upper_bound(dummyVec1.begin(),dummyVec1.end(), k);
820 std::vector<double>::iterator t1 = t2-1;
821
822 std::vector<double>::iterator e12 = std::upper_bound(aVecMap1[(*t1)].begin(),aVecMap1[(*t1)].end(), theta);
823 std::vector<double>::iterator e11 = e12-1;
824
825 std::vector<double>::iterator e22 = std::upper_bound(aVecMap1[(*t2)].begin(),aVecMap1[(*t2)].end(), theta);
826 std::vector<double>::iterator e21 = e22-1;
827
828 valueT1 =*t1;
829 valueT2 =*t2;
830 valueE21 =*e21;
831 valueE22 =*e22;
832 valueE12 =*e12;
833 valueE11 =*e11;
834
835 xs11 = FL1Data[valueT1][valueE11];
836 xs12 = FL1Data[valueT1][valueE12];
837 xs21 = FL1Data[valueT2][valueE21];
838 xs22 = FL1Data[valueT2][valueE22];
839
840 if (verboseLevel>0)
841 G4cout
842 << valueT1 << " "
843 << valueT2 << " "
844 << valueE11 << " "
845 << valueE12 << " "
846 << valueE21 << " "
847 << valueE22 << " "
848 << xs11 << " "
849 << xs12 << " "
850 << xs21 << " "
851 << xs22 << " "
852 << G4endl;
853
854 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
855
856 if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) return (0.);
857
858 if (xsProduct != 0.)
859 {
860 sigma = QuadInterpolator( valueE11, valueE12,
861 valueE21, valueE22,
862 xs11, xs12,
863 xs21, xs22,
864 valueT1, valueT2,
865 k, theta );
866 }
867
868 return sigma;
869}
870
871//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
872
873G4double G4ecpssrBaseLixsModel::FunctionFL2(G4double k, G4double theta)
874{
875
876 G4double sigma = 0.;
877 G4double valueT1 = 0;
878 G4double valueT2 = 0;
879 G4double valueE21 = 0;
880 G4double valueE22 = 0;
881 G4double valueE12 = 0;
882 G4double valueE11 = 0;
883 G4double xs11 = 0;
884 G4double xs12 = 0;
885 G4double xs21 = 0;
886 G4double xs22 = 0;
887
888 // PROTECTION TO ALLOW INTERPOLATION AT MINIMUM AND MAXIMUM Eta/Theta2 values
889
890 if (
891 theta==8.66e-4 ||
892 theta==8.66e-3 ||
893 theta==8.66e-2 ||
894 theta==8.66e-1 ||
895 theta==8.66e+00 ||
896 theta==8.66e+01
897 ) theta=theta-1e-12;
898
899 if (
900 theta==1.e-4 ||
901 theta==1.e-3 ||
902 theta==1.e-2 ||
903 theta==1.e-1 ||
904 theta==1.e+00 ||
905 theta==1.e+01
906 ) theta=theta+1e-12;
907
908 // END PROTECTION
909
910 std::vector<double>::iterator t2 = std::upper_bound(dummyVec2.begin(),dummyVec2.end(), k);
911 std::vector<double>::iterator t1 = t2-1;
912
913 std::vector<double>::iterator e12 = std::upper_bound(aVecMap2[(*t1)].begin(),aVecMap2[(*t1)].end(), theta);
914 std::vector<double>::iterator e11 = e12-1;
915
916 std::vector<double>::iterator e22 = std::upper_bound(aVecMap2[(*t2)].begin(),aVecMap2[(*t2)].end(), theta);
917 std::vector<double>::iterator e21 = e22-1;
918
919 valueT1 =*t1;
920 valueT2 =*t2;
921 valueE21 =*e21;
922 valueE22 =*e22;
923 valueE12 =*e12;
924 valueE11 =*e11;
925
926 xs11 = FL2Data[valueT1][valueE11];
927 xs12 = FL2Data[valueT1][valueE12];
928 xs21 = FL2Data[valueT2][valueE21];
929 xs22 = FL2Data[valueT2][valueE22];
930
931 if (verboseLevel>0)
932 G4cout
933 << valueT1 << " "
934 << valueT2 << " "
935 << valueE11 << " "
936 << valueE12 << " "
937 << valueE21 << " "
938 << valueE22 << " "
939 << xs11 << " "
940 << xs12 << " "
941 << xs21 << " "
942 << xs22 << " "
943 << G4endl;
944
945 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
946
947 if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) return (0.);
948
949 if (xsProduct != 0.)
950 {
951 sigma = QuadInterpolator( valueE11, valueE12,
952 valueE21, valueE22,
953 xs11, xs12,
954 xs21, xs22,
955 valueT1, valueT2,
956 k, theta );
957 }
958
959 return sigma;
960}
961
962//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
963
964G4double G4ecpssrBaseLixsModel::LinLinInterpolate(G4double e1,
965 G4double e2,
966 G4double e,
967 G4double xs1,
968 G4double xs2)
969{
970 G4double value = xs1 + (xs2 - xs1)*(e - e1)/ (e2 - e1);
971 return value;
972}
973
974//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
975
976G4double G4ecpssrBaseLixsModel::LinLogInterpolate(G4double e1,
977 G4double e2,
978 G4double e,
979 G4double xs1,
980 G4double xs2)
981{
982 G4double d1 = std::log(xs1);
983 G4double d2 = std::log(xs2);
984 G4double value = std::exp(d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
985 return value;
986}
987
988//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
989
990G4double G4ecpssrBaseLixsModel::LogLogInterpolate(G4double e1,
991 G4double e2,
992 G4double e,
993 G4double xs1,
994 G4double xs2)
995{
996 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
997 G4double b = std::log10(xs2) - a*std::log10(e2);
998 G4double sigma = a*std::log10(e) + b;
999 G4double value = (std::pow(10.,sigma));
1000 return value;
1001}
1002
1003//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1004
1005G4double G4ecpssrBaseLixsModel::QuadInterpolator(G4double e11, G4double e12,
1006 G4double e21, G4double e22,
1007 G4double xs11, G4double xs12,
1008 G4double xs21, G4double xs22,
1009 G4double t1, G4double t2,
1010 G4double t, G4double e)
1011{
1012// Log-Log
1013 G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
1014 G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
1015 G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
1016
1017/*
1018// Lin-Log
1019 G4double interpolatedvalue1 = LinLogInterpolate(e11, e12, e, xs11, xs12);
1020 G4double interpolatedvalue2 = LinLogInterpolate(e21, e22, e, xs21, xs22);
1021 G4double value = LinLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
1022*/
1023
1024/*
1025// Lin-Lin
1026 G4double interpolatedvalue1 = LinLinInterpolate(e11, e12, e, xs11, xs12);
1027 G4double interpolatedvalue2 = LinLinInterpolate(e21, e22, e, xs21, xs22);
1028 G4double value = LinLinInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
1029*/
1030 return value;
1031
1032}
1033
@ FatalException
double G4double
Definition: G4Types.hh:64
int G4int
Definition: G4Types.hh:66
#define G4endl
Definition: G4ios.hh:52
G4DLLIMPORT std::ostream G4cout
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
G4double BindingEnergy() const
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const
static G4AtomicTransitionManager * Instance()
static G4NistManager * Instance()
G4double GetAtomicMassAmu(const G4String &symb) const
G4double GetPDGCharge() const
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4double CalculateL1CrossSection(G4int zTarget, G4double massIncident, G4double energyIncident)
G4double CalculateL2CrossSection(G4int zTarget, G4double massIncident, G4double energyIncident)
G4double ExpIntFunction(G4int n, G4double x)
G4double CalculateVelocity(G4int subShell, G4int zTarget, G4double massIncident, G4double energyIncident)
G4double CalculateL3CrossSection(G4int zTarget, G4double massIncident, G4double energyIncident)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41