ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ecpssrBaseLixsModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ecpssrBaseLixsModel.cc
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 
29 #include "G4ecpssrBaseLixsModel.hh"
30 
31 #include "globals.hh"
32 #include "G4PhysicalConstants.hh"
33 #include "G4SystemOfUnits.hh"
34 #include "G4AtomicTransitionManager.hh"
35 #include "G4NistManager.hh"
36 #include "G4Proton.hh"
37 #include "G4Alpha.hh"
38 #include "G4LinLogInterpolation.hh"
39 #include "G4Exp.hh"
40 
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
42 
43 G4ecpssrBaseLixsModel::G4ecpssrBaseLixsModel()
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 
116 G4ecpssrBaseLixsModel::~G4ecpssrBaseLixsModel()
117 {}
118 
119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
120 
121 G4double G4ecpssrBaseLixsModel::ExpIntFunction(G4int n,G4double x)
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 
191 G4double G4ecpssrBaseLixsModel::CalculateL1CrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
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 
201  G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
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 
394 G4double G4ecpssrBaseLixsModel::CalculateL2CrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
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 
404  G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
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 
557 G4double G4ecpssrBaseLixsModel::CalculateL3CrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
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 
567  G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
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 
727 G4double G4ecpssrBaseLixsModel::CalculateVelocity(G4int subShell, G4int zTarget, G4double massIncident, G4double energyIncident)
728 
729 {
730 
731  G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
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 
765 G4double 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 
856 G4double 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 
947 G4double 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 
959 G4double 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 
973 G4double 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 
988 G4double 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 }