Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions
G4ecpssrBaseLixsModel Class Reference

#include <G4ecpssrBaseLixsModel.hh>

Inheritance diagram for G4ecpssrBaseLixsModel:
G4VecpssrLiModel

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.

References FatalException, and G4Exception().

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

Definition at line 117 of file G4ecpssrBaseLixsModel.cc.

118 {}

Member Function Documentation

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

Implements G4VecpssrLiModel.

Definition at line 192 of file G4ecpssrBaseLixsModel.cc.

References G4Alpha::Alpha(), python.hepunit::amu_c2, python.hepunit::barn, G4AtomicShell::BindingEnergy(), python.hepunit::Bohr_radius, CalculateVelocity(), python.hepunit::electron_mass_c2, python.hepunit::eplus, ExpIntFunction(), G4cout, G4endl, G4NistManager::GetAtomicMassAmu(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), G4AtomicTransitionManager::Instance(), G4NistManager::Instance(), python.hepunit::pi, G4Proton::Proton(), and G4AtomicTransitionManager::Shell().

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 
202  G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
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 }
G4ecpssrBaseLixsModel::ExpIntFunction
G4double ExpIntFunction(G4int n, G4double x)
Definition: G4ecpssrBaseLixsModel.cc:122
G4AtomicShell::BindingEnergy
G4double BindingEnergy() const
Definition: G4AtomicShell.hh:79
G4ecpssrBaseLixsModel::CalculateVelocity
G4double CalculateVelocity(G4int subShell, G4int zTarget, G4double massIncident, G4double energyIncident)
Definition: G4ecpssrBaseLixsModel.cc:744
G4NistManager::Instance
static G4NistManager * Instance()
Definition: G4NistManager.cc:68
G4cout
G4GLOB_DLL std::ostream G4cout
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161
G4NistManager::GetAtomicMassAmu
G4double GetAtomicMassAmu(const G4String &symb) const
Definition: G4NistManager.hh:344
G4endl
#define G4endl
Definition: G4ios.hh:61
G4Alpha::Alpha
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
G4double
double G4double
Definition: G4Types.hh:76
python.hepunit.amu_c2
float amu_c2
Definition: hepunit.py:277
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:163
G4AtomicTransitionManager::Instance
static G4AtomicTransitionManager * Instance()
Definition: G4AtomicTransitionManager.cc:162
G4AtomicTransitionManager::Shell
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const
Definition: G4AtomicTransitionManager.cc:173
G4double G4ecpssrBaseLixsModel::CalculateL2CrossSection ( G4int  zTarget,
G4double  massIncident,
G4double  energyIncident 
)
virtual

Implements G4VecpssrLiModel.

Definition at line 407 of file G4ecpssrBaseLixsModel.cc.

References G4Alpha::Alpha(), python.hepunit::amu_c2, python.hepunit::barn, G4AtomicShell::BindingEnergy(), python.hepunit::Bohr_radius, CalculateVelocity(), python.hepunit::electron_mass_c2, python.hepunit::eplus, ExpIntFunction(), G4cout, G4endl, G4NistManager::GetAtomicMassAmu(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), G4AtomicTransitionManager::Instance(), G4NistManager::Instance(), python.hepunit::pi, G4Proton::Proton(), and G4AtomicTransitionManager::Shell().

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 
417  G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
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 }
G4ecpssrBaseLixsModel::ExpIntFunction
G4double ExpIntFunction(G4int n, G4double x)
Definition: G4ecpssrBaseLixsModel.cc:122
G4AtomicShell::BindingEnergy
G4double BindingEnergy() const
Definition: G4AtomicShell.hh:79
G4ecpssrBaseLixsModel::CalculateVelocity
G4double CalculateVelocity(G4int subShell, G4int zTarget, G4double massIncident, G4double energyIncident)
Definition: G4ecpssrBaseLixsModel.cc:744
G4NistManager::Instance
static G4NistManager * Instance()
Definition: G4NistManager.cc:68
G4cout
G4GLOB_DLL std::ostream G4cout
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161
G4NistManager::GetAtomicMassAmu
G4double GetAtomicMassAmu(const G4String &symb) const
Definition: G4NistManager.hh:344
G4endl
#define G4endl
Definition: G4ios.hh:61
G4Alpha::Alpha
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
G4double
double G4double
Definition: G4Types.hh:76
python.hepunit.amu_c2
float amu_c2
Definition: hepunit.py:277
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:163
G4AtomicTransitionManager::Instance
static G4AtomicTransitionManager * Instance()
Definition: G4AtomicTransitionManager.cc:162
G4AtomicTransitionManager::Shell
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const
Definition: G4AtomicTransitionManager.cc:173
G4double G4ecpssrBaseLixsModel::CalculateL3CrossSection ( G4int  zTarget,
G4double  massIncident,
G4double  energyIncident 
)
virtual

Implements G4VecpssrLiModel.

Definition at line 574 of file G4ecpssrBaseLixsModel.cc.

References G4Alpha::Alpha(), python.hepunit::amu_c2, python.hepunit::barn, G4AtomicShell::BindingEnergy(), python.hepunit::Bohr_radius, CalculateVelocity(), python.hepunit::electron_mass_c2, python.hepunit::eplus, ExpIntFunction(), G4cout, G4endl, G4NistManager::GetAtomicMassAmu(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), G4AtomicTransitionManager::Instance(), G4NistManager::Instance(), python.hepunit::pi, G4Proton::Proton(), and G4AtomicTransitionManager::Shell().

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 
584  G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
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 }
G4ecpssrBaseLixsModel::ExpIntFunction
G4double ExpIntFunction(G4int n, G4double x)
Definition: G4ecpssrBaseLixsModel.cc:122
G4AtomicShell::BindingEnergy
G4double BindingEnergy() const
Definition: G4AtomicShell.hh:79
G4ecpssrBaseLixsModel::CalculateVelocity
G4double CalculateVelocity(G4int subShell, G4int zTarget, G4double massIncident, G4double energyIncident)
Definition: G4ecpssrBaseLixsModel.cc:744
G4NistManager::Instance
static G4NistManager * Instance()
Definition: G4NistManager.cc:68
G4cout
G4GLOB_DLL std::ostream G4cout
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161
G4NistManager::GetAtomicMassAmu
G4double GetAtomicMassAmu(const G4String &symb) const
Definition: G4NistManager.hh:344
G4endl
#define G4endl
Definition: G4ios.hh:61
G4Alpha::Alpha
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
G4double
double G4double
Definition: G4Types.hh:76
python.hepunit.amu_c2
float amu_c2
Definition: hepunit.py:277
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:163
G4AtomicTransitionManager::Instance
static G4AtomicTransitionManager * Instance()
Definition: G4AtomicTransitionManager.cc:162
G4AtomicTransitionManager::Shell
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const
Definition: G4AtomicTransitionManager.cc:173
G4double G4ecpssrBaseLixsModel::CalculateVelocity ( G4int  subShell,
G4int  zTarget,
G4double  massIncident,
G4double  energyIncident 
)

Definition at line 744 of file G4ecpssrBaseLixsModel.cc.

References G4Alpha::Alpha(), G4AtomicShell::BindingEnergy(), python.hepunit::electron_mass_c2, G4cout, G4endl, G4ParticleDefinition::GetPDGMass(), G4AtomicTransitionManager::Instance(), G4Proton::Proton(), and G4AtomicTransitionManager::Shell().

Referenced by CalculateL1CrossSection(), CalculateL2CrossSection(), and CalculateL3CrossSection().

746 {
747 
748  G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
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 }
G4AtomicShell::BindingEnergy
G4double BindingEnergy() const
Definition: G4AtomicShell.hh:79
G4cout
G4GLOB_DLL std::ostream G4cout
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161
G4endl
#define G4endl
Definition: G4ios.hh:61
G4Alpha::Alpha
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
G4double
double G4double
Definition: G4Types.hh:76
G4AtomicTransitionManager::Instance
static G4AtomicTransitionManager * Instance()
Definition: G4AtomicTransitionManager.cc:162
G4AtomicTransitionManager::Shell
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const
Definition: G4AtomicTransitionManager.cc:173
G4double G4ecpssrBaseLixsModel::ExpIntFunction ( G4int  n,
G4double  x 
)

Definition at line 122 of file G4ecpssrBaseLixsModel.cc.

References test::a, test::b, test::c, G4cout, G4endl, n, and test::x.

Referenced by CalculateL1CrossSection(), CalculateL2CrossSection(), and CalculateL3CrossSection().

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 }
G4int
int G4int
Definition: G4Types.hh:78
G4cout
G4GLOB_DLL std::ostream G4cout
n
const G4int n
Definition: G4UrQMD1_3Interface.hh:144
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
The documentation for this class was generated from the following files:
Generated on Wed Apr 30 2014 15:56:21 for Geant4.10 by   doxygen 1.8.7