Geant4  10.01.p02
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 
30 #include "G4ecpssrBaseLixsModel.hh"
31 
32 #include "globals.hh"
33 #include "G4PhysicalConstants.hh"
34 #include "G4SystemOfUnits.hh"
36 #include "G4NistManager.hh"
37 #include "G4Proton.hh"
38 #include "G4Alpha.hh"
39 #include "G4LinLogInterpolation.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  }
187 return 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 
744 G4double 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 
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 
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 
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 
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 
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 
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 
G4double LinLinInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
static const G4double d1
G4double ExpIntFunction(G4int n, G4double x)
G4double FunctionFL1(G4double k, G4double theta)
const G4double pi
std::vector< double > dummyVec1
G4double FunctionFL2(G4double k, G4double theta)
static const G4double e2
G4double QuadInterpolator(G4double e11, G4double e12, G4double e21, G4double e22, G4double x11, G4double x12, G4double x21, G4double x22, G4double t1, G4double t2, G4double t, G4double e)
static const G4double eps
G4double a
Definition: TRTMaterials.hh:39
G4double BindingEnergy() const
int G4int
Definition: G4Types.hh:78
G4double CalculateVelocity(G4int subShell, G4int zTarget, G4double massIncident, G4double energyIncident)
static G4NistManager * Instance()
G4double CalculateL2CrossSection(G4int zTarget, G4double massIncident, G4double energyIncident)
G4double CalculateL1CrossSection(G4int zTarget, G4double massIncident, G4double energyIncident)
G4GLOB_DLL std::ostream G4cout
G4double CalculateL3CrossSection(G4int zTarget, G4double massIncident, G4double energyIncident)
static G4Proton * Proton()
Definition: G4Proton.cc:93
const G4int n
static const G4double e1
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
std::vector< double > dummyVec2
G4double LogLogInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
G4double GetPDGMass() const
G4double GetAtomicMassAmu(const G4String &symb) const
G4double LinLogInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
#define G4endl
Definition: G4ios.hh:61
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
static const double barn
Definition: G4SIunits.hh:95
double G4double
Definition: G4Types.hh:76
static const double eplus
Definition: G4SIunits.hh:178
G4double GetPDGCharge() const
static const G4double d2
static G4AtomicTransitionManager * Instance()
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const