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