Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ecpssrBaseLixsModel Class Reference

#include <G4ecpssrBaseLixsModel.hh>

Inheritance diagram for G4ecpssrBaseLixsModel:
Collaboration diagram for G4ecpssrBaseLixsModel:

Public Member Functions

 G4ecpssrBaseLixsModel ()
 
 ~G4ecpssrBaseLixsModel ()
 
G4double CalculateL1CrossSection (G4int zTarget, G4double massIncident, G4double energyIncident)
 
G4double CalculateL2CrossSection (G4int zTarget, G4double massIncident, G4double energyIncident)
 
G4double CalculateL3CrossSection (G4int zTarget, G4double massIncident, G4double energyIncident)
 
G4double CalculateVelocity (G4int subShell, G4int zTarget, G4double massIncident, G4double energyIncident)
 
G4double ExpIntFunction (G4int n, G4double x)
 
- Public Member Functions inherited from G4VecpssrLiModel
 G4VecpssrLiModel ()
 
virtual ~G4VecpssrLiModel ()
 

Detailed Description

Definition at line 55 of file G4ecpssrBaseLixsModel.hh.

Constructor & Destructor Documentation

G4ecpssrBaseLixsModel::G4ecpssrBaseLixsModel ( )

Definition at line 43 of file G4ecpssrBaseLixsModel.cc.

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 }
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41

Here is the call graph for this function:

G4ecpssrBaseLixsModel::~G4ecpssrBaseLixsModel ( )

Definition at line 116 of file G4ecpssrBaseLixsModel.cc.

117 {}

Member Function Documentation

G4double G4ecpssrBaseLixsModel::CalculateL1CrossSection ( G4int  zTarget,
G4double  massIncident,
G4double  energyIncident 
)
virtual

Implements G4VecpssrLiModel.

Definition at line 191 of file G4ecpssrBaseLixsModel.cc.

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 }
G4double ExpIntFunction(G4int n, G4double x)
G4double BindingEnergy() const
G4double CalculateVelocity(G4int subShell, G4int zTarget, G4double massIncident, G4double energyIncident)
static G4NistManager * Instance()
G4GLOB_DLL std::ostream G4cout
static constexpr double eplus
Definition: G4SIunits.hh:199
static G4Proton * Proton()
Definition: G4Proton.cc:93
float electron_mass_c2
Definition: hepunit.py:274
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double GetPDGMass() const
G4double GetAtomicMassAmu(const G4String &symb) const
#define G4endl
Definition: G4ios.hh:61
static constexpr double pi
Definition: G4SIunits.hh:75
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
double G4double
Definition: G4Types.hh:76
float amu_c2
Definition: hepunit.py:277
static constexpr double barn
Definition: G4SIunits.hh:105
G4double GetPDGCharge() const
static G4AtomicTransitionManager * Instance()
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const

Here is the call graph for this function:

G4double G4ecpssrBaseLixsModel::CalculateL2CrossSection ( G4int  zTarget,
G4double  massIncident,
G4double  energyIncident 
)
virtual

Implements G4VecpssrLiModel.

Definition at line 394 of file G4ecpssrBaseLixsModel.cc.

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 }
G4double ExpIntFunction(G4int n, G4double x)
G4double BindingEnergy() const
G4double CalculateVelocity(G4int subShell, G4int zTarget, G4double massIncident, G4double energyIncident)
static G4NistManager * Instance()
G4GLOB_DLL std::ostream G4cout
static constexpr double eplus
Definition: G4SIunits.hh:199
static G4Proton * Proton()
Definition: G4Proton.cc:93
float electron_mass_c2
Definition: hepunit.py:274
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double GetPDGMass() const
G4double GetAtomicMassAmu(const G4String &symb) const
#define G4endl
Definition: G4ios.hh:61
static constexpr double pi
Definition: G4SIunits.hh:75
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
double G4double
Definition: G4Types.hh:76
float amu_c2
Definition: hepunit.py:277
static constexpr double barn
Definition: G4SIunits.hh:105
G4double GetPDGCharge() const
static G4AtomicTransitionManager * Instance()
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const

Here is the call graph for this function:

G4double G4ecpssrBaseLixsModel::CalculateL3CrossSection ( G4int  zTarget,
G4double  massIncident,
G4double  energyIncident 
)
virtual

Implements G4VecpssrLiModel.

Definition at line 557 of file G4ecpssrBaseLixsModel.cc.

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 }
G4double ExpIntFunction(G4int n, G4double x)
G4double BindingEnergy() const
G4double CalculateVelocity(G4int subShell, G4int zTarget, G4double massIncident, G4double energyIncident)
static G4NistManager * Instance()
G4GLOB_DLL std::ostream G4cout
static constexpr double eplus
Definition: G4SIunits.hh:199
static G4Proton * Proton()
Definition: G4Proton.cc:93
float electron_mass_c2
Definition: hepunit.py:274
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double GetPDGMass() const
G4double GetAtomicMassAmu(const G4String &symb) const
#define G4endl
Definition: G4ios.hh:61
static constexpr double pi
Definition: G4SIunits.hh:75
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
double G4double
Definition: G4Types.hh:76
float amu_c2
Definition: hepunit.py:277
static constexpr double barn
Definition: G4SIunits.hh:105
G4double GetPDGCharge() const
static G4AtomicTransitionManager * Instance()
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const

Here is the call graph for this function:

G4double G4ecpssrBaseLixsModel::CalculateVelocity ( G4int  subShell,
G4int  zTarget,
G4double  massIncident,
G4double  energyIncident 
)

Definition at line 727 of file G4ecpssrBaseLixsModel.cc.

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 }
G4double BindingEnergy() const
G4GLOB_DLL std::ostream G4cout
static G4Proton * Proton()
Definition: G4Proton.cc:93
float electron_mass_c2
Definition: hepunit.py:274
G4double GetPDGMass() const
#define G4endl
Definition: G4ios.hh:61
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
double G4double
Definition: G4Types.hh:76
static G4AtomicTransitionManager * Instance()
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4ecpssrBaseLixsModel::ExpIntFunction ( G4int  n,
G4double  x 
)

Definition at line 121 of file G4ecpssrBaseLixsModel.cc.

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 }
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
static const G4double eps
tuple x
Definition: test.py:50
int G4int
Definition: G4Types.hh:78
tuple b
Definition: test.py:12
G4GLOB_DLL std::ostream G4cout
const G4int n
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
tuple c
Definition: test.py:13

Here is the call graph for this function:

Here is the caller graph for this function:


The documentation for this class was generated from the following files: