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)

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().

00044 {
00045     verboseLevel=0;
00046 
00047     // Storing FLi data needed for 0.2 to 3.0  velocities region
00048 
00049     char *path = getenv("G4LEDATA");
00050 
00051     if (!path) {      
00052       G4Exception("G4ecpssrLCrossSection::G4ecpssrBaseLixsModel()","em0006", FatalException ,"G4LEDATA environment variable not set");
00053       return;
00054     }
00055     std::ostringstream fileName1;
00056     std::ostringstream fileName2;
00057     
00058     fileName1 << path << "/pixe/uf/FL1.dat";
00059     fileName2 << path << "/pixe/uf/FL2.dat";
00060 
00061     // Reading of FL1.dat
00062 
00063     std::ifstream FL1(fileName1.str().c_str());
00064     if (!FL1) G4Exception("G4ecpssrLCrossSection::G4ecpssrBaseLixsModel()","em0003",FatalException, "error opening FL1 data file");
00065   
00066     dummyVec1.push_back(0.);
00067 
00068     while(!FL1.eof())
00069     {
00070         double x1;
00071         double y1;
00072 
00073         FL1>>x1>>y1;
00074 
00075         //  Mandatory vector initialization
00076         if (x1 != dummyVec1.back())
00077         {
00078           dummyVec1.push_back(x1);
00079           aVecMap1[x1].push_back(-1.);
00080         }
00081 
00082         FL1>>FL1Data[x1][y1];
00083 
00084         if (y1 != aVecMap1[x1].back()) aVecMap1[x1].push_back(y1);
00085     }
00086 
00087     // Reading of FL2.dat
00088     
00089     std::ifstream FL2(fileName2.str().c_str());
00090     if (!FL2) G4Exception("G4ecpssrLCrossSection::G4ecpssrBaseLixsModel()","em0003", FatalException," error opening FL2 data file");
00091     
00092     dummyVec2.push_back(0.);
00093 
00094     while(!FL2.eof())
00095     {
00096         double x2;
00097         double y2;
00098 
00099         FL2>>x2>>y2;
00100 
00101         //  Mandatory vector initialization
00102         if (x2 != dummyVec2.back())
00103         {
00104           dummyVec2.push_back(x2);
00105           aVecMap2[x2].push_back(-1.);
00106         }
00107 
00108         FL2>>FL2Data[x2][y2];
00109 
00110         if (y2 != aVecMap2[x2].back()) aVecMap2[x2].push_back(y2);
00111     }
00112 
00113 }

G4ecpssrBaseLixsModel::~G4ecpssrBaseLixsModel (  ) 

Definition at line 117 of file G4ecpssrBaseLixsModel.cc.

00118 {}


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(), G4AtomicShell::BindingEnergy(), CalculateVelocity(), ExpIntFunction(), G4cout, G4endl, G4NistManager::GetAtomicMassAmu(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), G4AtomicTransitionManager::Instance(), G4NistManager::Instance(), G4INCL::Math::pi, G4Proton::Proton(), and G4AtomicTransitionManager::Shell().

00193 {
00194 
00195   if (zTarget <=4) return 0.;
00196 
00197  //this L1-CrossSection calculation method is done according to Werner Brandt and Grzegorz Lapicki, Phys.Rev.A20 N2 (1979),
00198  //and using data tables of O. Benka et al. At.Data Nucl.Data Tables Vol.22 No.3 (1978).
00199 
00200   G4NistManager* massManager = G4NistManager::Instance();
00201 
00202   G4AtomicTransitionManager*  transitionManager =  G4AtomicTransitionManager::Instance();
00203 
00204   G4double  zIncident = 0;
00205   G4Proton* aProtone = G4Proton::Proton();
00206   G4Alpha* aAlpha = G4Alpha::Alpha();
00207 
00208   if (massIncident == aProtone->GetPDGMass() )
00209 
00210    zIncident = (aProtone->GetPDGCharge())/eplus;
00211 
00212   else
00213     {
00214       if (massIncident == aAlpha->GetPDGMass())
00215 
00216           zIncident  = (aAlpha->GetPDGCharge())/eplus;
00217 
00218       else
00219         {
00220           G4cout << "*** WARNING in G4ecpssrBaseLixsModel::CalculateL1CrossSection : Proton or Alpha incident particles only. " << G4endl;
00221           G4cout << massIncident << ", " << aAlpha->GetPDGMass() << " (alpha)" << aProtone->GetPDGMass() << " (proton)" << G4endl;
00222           return 0;
00223         }
00224     }
00225 
00226   G4double l1BindingEnergy = transitionManager->Shell(zTarget,1)->BindingEnergy(); //Observed binding energy of L1-subshell
00227 
00228   G4double massTarget = (massManager->GetAtomicMassAmu(zTarget))*amu_c2;
00229 
00230   G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2; //Mass of the system (projectile, target)
00231 
00232   const G4double zlshell= 4.15;
00233   // *** see Benka, ADANDT 22, p 223
00234 
00235   G4double screenedzTarget = zTarget-zlshell; //Effective nuclear charge as seen by electrons in L1-sub shell
00236 
00237   const G4double rydbergMeV= 13.6056923e-6;
00238 
00239   const G4double nl= 2.;
00240   // *** see Benka, ADANDT 22, p 220, f3
00241 
00242   G4double tetal1 = (l1BindingEnergy*nl*nl)/((screenedzTarget*screenedzTarget)*rydbergMeV); //Screening parameter
00243   // *** see Benka, ADANDT 22, p 220, f3
00244 
00245     if (verboseLevel>0) G4cout << "  tetal1=" <<  tetal1<< G4endl;
00246 
00247   G4double reducedEnergy = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
00248   // *** also called etaS
00249   // *** see Benka, ADANDT 22, p 220, f3
00250 
00251   const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/barn ; //Bohr radius of hydrogen
00252 
00253   G4double sigma0 = 8.*pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.);
00254   // *** see Benka, ADANDT 22, p 220, f2, for protons
00255   // *** see Basbas, Phys Rev A7, p 1000
00256 
00257   G4double velocityl1 = CalculateVelocity(1, zTarget, massIncident, energyIncident); // Scaled velocity
00258 
00259     if (verboseLevel>0) G4cout << "  velocityl1=" << velocityl1<< G4endl;
00260 
00261   const G4double l1AnalyticalApproximation= 1.5;
00262   G4double x1 =(nl*l1AnalyticalApproximation)/velocityl1;
00263   // *** 1.5 is cK = cL1 (it is 1.25 for L2 & L3)
00264   // *** see Brandt, Phys Rev A20, p 469, f16 in expression of h
00265   
00266    if (verboseLevel>0) G4cout << "  x1=" << x1<< G4endl;
00267 
00268   G4double electrIonizationEnergyl1=0.;
00269   // *** see Basbas, Phys Rev A17, p1665, f27
00270   // *** see Brandt, Phys Rev A20, p469
00271   // *** see Liu, Comp Phys Comm 97, p325, f A5
00272 
00273   if ( x1<=0.035)  electrIonizationEnergyl1= 0.75*pi*(std::log(1./(x1*x1))-1.);
00274   else
00275     {
00276       if ( x1<=3.)
00277         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));
00278       else
00279         {if ( x1<=11.) electrIonizationEnergyl1 =2.*std::exp(-2.*x1)/std::pow(x1,1.6);}
00280     }
00281 
00282   G4double hFunctionl1 =(electrIonizationEnergyl1*2.*nl)/(tetal1*std::pow(velocityl1,3)); //takes into account the polarization effect
00283   // *** see Brandt, Phys Rev A20, p 469, f16
00284 
00285     if (verboseLevel>0) G4cout << "  hFunctionl1=" << hFunctionl1<< G4endl;
00286 
00287   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
00288   // *** see Brandt, Phys Rev A20, p 469, f19
00289 
00290     if (verboseLevel>0) G4cout << "  gFunctionl1=" << gFunctionl1<< G4endl;
00291 
00292   G4double sigmaPSS_l1 = 1.+(((2.*zIncident)/(screenedzTarget*tetal1))*(gFunctionl1-hFunctionl1)); //Binding-polarization factor
00293   // *** also called dzeta
00294   // *** also called epsilon
00295   // *** see Basbas, Phys Rev A17, p1667, f45
00296 
00297     if (verboseLevel>0) G4cout << "sigmaPSS_l1 =" << sigmaPSS_l1<< G4endl;
00298 
00299   const G4double cNaturalUnit= 137.;
00300   
00301   G4double yl1Formula=0.4*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(nl*velocityl1/sigmaPSS_l1);
00302   // *** also called yS
00303   // *** see Brandt, Phys Rev A20, p467, f6
00304   // *** see Brandt, Phys Rev A23, p1728
00305       
00306   G4double l1relativityCorrection = std::pow((1.+(1.1*yl1Formula*yl1Formula)),0.5)+yl1Formula; // Relativistic correction parameter
00307   // *** also called mRS
00308   // *** see Brandt, Phys Rev A20, p467, f6
00309   
00310   //G4double reducedVelocity_l1 = velocityl1*std::pow(l1relativityCorrection,0.5); //Reduced velocity parameter
00311 
00312   G4double L1etaOverTheta2;
00313 
00314   G4double  universalFunction_l1 = 0.;
00315 
00316   G4double sigmaPSSR_l1;
00317 
00318   // low velocity formula
00319   // *****************  
00320   if ( velocityl1 <20.  )
00321   {
00322 
00323     L1etaOverTheta2 =(reducedEnergy* l1relativityCorrection)/((tetal1*sigmaPSS_l1)*(tetal1*sigmaPSS_l1));
00324     // *** 1) RELATIVISTIC CORRECTION ADDED
00325     // *** 2) sigma_PSS_l1 ADDED
00326     // *** reducedEnergy is etaS, l1relativityCorrection is mRS
00327     // *** see Phys Rev A20, p468, top
00328     
00329     if ( ((tetal1*sigmaPSS_l1) >=0.2) && ((tetal1*sigmaPSS_l1) <=2.6670) && (L1etaOverTheta2>=0.1e-3) && (L1etaOverTheta2<=0.866e2) )
00330   
00331       universalFunction_l1 = FunctionFL1((tetal1*sigmaPSS_l1),L1etaOverTheta2);
00332 
00333       if (verboseLevel>0) G4cout << "at low velocity range, universalFunction_l1  =" << universalFunction_l1 << G4endl;
00334 
00335     sigmaPSSR_l1 = (sigma0/(tetal1*sigmaPSS_l1))*universalFunction_l1;// Plane-wave Born -Aproximation L1-subshell ionisation Cross Section
00336   // *** see Benka, ADANDT 22, p220, f1
00337 
00338       if (verboseLevel>0) G4cout << "  at low velocity range, sigma PWBA L1 CS  = " << sigmaPSSR_l1<< G4endl;
00339 
00340    }
00341 
00342   else
00343 
00344   {
00345 
00346     L1etaOverTheta2 = reducedEnergy/(tetal1*tetal1);
00347     // Medium & high velocity
00348     // *** 1) NO RELATIVISTIC CORRECTION
00349     // *** 2) NO sigma_PSS_l1
00350     // *** see Benka, ADANDT 22, p223
00351 
00352     if ( (tetal1 >=0.2) && (tetal1 <=2.6670) && (L1etaOverTheta2>=0.1e-3) && (L1etaOverTheta2<=0.866e2) ) 
00353     
00354       universalFunction_l1 = FunctionFL1(tetal1,L1etaOverTheta2);
00355 
00356     if (verboseLevel>0) G4cout << "at medium and high velocity range, universalFunction_l1  =" << universalFunction_l1 << G4endl;
00357 
00358     sigmaPSSR_l1 = (sigma0/tetal1)*universalFunction_l1;// Plane-wave Born -Aproximation L1-subshell ionisation Cross Section
00359   // *** see Benka, ADANDT 22, p220, f1
00360 
00361     if (verboseLevel>0) G4cout << "  sigma PWBA L1 CS at medium and high velocity range = " << sigmaPSSR_l1<< G4endl;    
00362   }
00363  
00364   G4double pssDeltal1 = (4./(systemMass *sigmaPSS_l1*tetal1))*(sigmaPSS_l1/velocityl1)*(sigmaPSS_l1/velocityl1);
00365   // *** also called dzeta*delta
00366   // *** see Brandt, Phys Rev A23, p1727, f B2
00367 
00368     if (verboseLevel>0) G4cout << "  pssDeltal1=" << pssDeltal1<< G4endl;
00369 
00370   if (pssDeltal1>1) return 0.;
00371 
00372   G4double energyLossl1 = std::pow(1-pssDeltal1,0.5);
00373   // *** also called z
00374   // *** see Brandt, Phys Rev A23, p1727, after f B2
00375 
00376     if (verboseLevel>0) G4cout << "  energyLossl1=" << energyLossl1<< G4endl;
00377 
00378   G4double coulombDeflectionl1 = 
00379    (8.*pi*zIncident/systemMass)*std::pow(tetal1*sigmaPSS_l1,-2.)*std::pow(velocityl1/sigmaPSS_l1,-3.)*(zTarget/screenedzTarget);
00380   // *** see Brandt, Phys Rev A20, v2s and f2 and B2 
00381   // *** with factor n2 compared to Brandt, Phys Rev A23, p1727, f B3
00382 
00383   G4double cParameterl1 =2.* coulombDeflectionl1/(energyLossl1*(energyLossl1+1.));
00384   // *** see Brandt, Phys Rev A23, p1727, f B4
00385 
00386   G4double coulombDeflectionFunction_l1 = 9.*ExpIntFunction(10,cParameterl1); //Coulomb-deflection effect correction
00387 
00388     if (verboseLevel>0) G4cout << "  coulombDeflectionFunction_l1 =" << coulombDeflectionFunction_l1 << G4endl;
00389 
00390   G4double crossSection_L1 = coulombDeflectionFunction_l1 * sigmaPSSR_l1;
00391 
00392   //ECPSSR L1 -subshell cross section is estimated at perturbed-stationnairy-state(PSS)
00393   //and reduced by the energy-loss(E),the Coulomb deflection(C),and the relativity(R) effects
00394 
00395     if (verboseLevel>0) G4cout << "  crossSection_L1 =" << crossSection_L1 << G4endl;
00396 
00397   if (crossSection_L1 >= 0) 
00398   {
00399     return crossSection_L1 * barn;
00400   }
00401   
00402   else {return 0;}
00403 }

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

Implements G4VecpssrLiModel.

Definition at line 407 of file G4ecpssrBaseLixsModel.cc.

References G4Alpha::Alpha(), G4AtomicShell::BindingEnergy(), CalculateVelocity(), ExpIntFunction(), G4cout, G4endl, G4NistManager::GetAtomicMassAmu(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), G4AtomicTransitionManager::Instance(), G4NistManager::Instance(), G4INCL::Math::pi, G4Proton::Proton(), and G4AtomicTransitionManager::Shell().

00409 {
00410   if (zTarget <=13 ) return 0.;
00411 
00412   // this L2-CrossSection calculation method is done according to Werner Brandt and Grzegorz Lapicki, Phys.Rev.A20 N2 (1979),
00413   // and using data tables of O. Benka et al. At.Data Nucl.Data Tables Vol.22 No.3 (1978).
00414 
00415   G4NistManager* massManager = G4NistManager::Instance();
00416 
00417   G4AtomicTransitionManager*  transitionManager =  G4AtomicTransitionManager::Instance();
00418 
00419   G4double  zIncident = 0;
00420   
00421   G4Proton* aProtone = G4Proton::Proton();
00422   G4Alpha* aAlpha = G4Alpha::Alpha();
00423 
00424   if (massIncident == aProtone->GetPDGMass() )
00425 
00426    zIncident = (aProtone->GetPDGCharge())/eplus;
00427 
00428   else
00429     {
00430       if (massIncident == aAlpha->GetPDGMass())
00431 
00432           zIncident  = (aAlpha->GetPDGCharge())/eplus;
00433 
00434       else
00435         {
00436           G4cout << "*** WARNING in G4ecpssrBaseLixsModel::CalculateL2CrossSection : Proton or Alpha incident particles only. " << G4endl;
00437           G4cout << massIncident << ", " << aAlpha->GetPDGMass() << " (alpha)" << aProtone->GetPDGMass() << " (proton)" << G4endl;
00438           return 0;
00439         }
00440     }
00441 
00442   G4double l2BindingEnergy = transitionManager->Shell(zTarget,2)->BindingEnergy(); //Observed binding energy of L2-subshell
00443 
00444   G4double massTarget = (massManager->GetAtomicMassAmu(zTarget))*amu_c2;
00445 
00446   G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2; //Mass of the system (projectile, target)
00447 
00448   const G4double zlshell= 4.15;
00449 
00450   G4double screenedzTarget = zTarget-zlshell; //Effective nuclear charge as seen by electrons in L2-subshell
00451 
00452   const G4double rydbergMeV= 13.6056923e-6;
00453 
00454   const G4double nl= 2.;
00455 
00456   G4double tetal2 = (l2BindingEnergy*nl*nl)/((screenedzTarget*screenedzTarget)*rydbergMeV); //Screening parameter
00457 
00458     if (verboseLevel>0) G4cout << "  tetal2=" <<  tetal2<< G4endl;
00459 
00460   G4double reducedEnergy = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget); 
00461 
00462   const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/barn ; //Bohr radius of hydrogen
00463 
00464   G4double sigma0 = 8.*pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.);
00465 
00466   G4double velocityl2 = CalculateVelocity(2,  zTarget, massIncident, energyIncident); // Scaled velocity
00467 
00468     if (verboseLevel>0) G4cout << "  velocityl2=" << velocityl2<< G4endl;
00469 
00470   const G4double l23AnalyticalApproximation= 1.25;
00471 
00472   G4double x2 = (nl*l23AnalyticalApproximation)/velocityl2;
00473     
00474     if (verboseLevel>0) G4cout << "  x2=" << x2<< G4endl;
00475 
00476   G4double electrIonizationEnergyl2=0.;
00477 
00478   if ( x2<=0.035)  electrIonizationEnergyl2= 0.75*pi*(std::log(1./(x2*x2))-1.);
00479   else
00480     {
00481       if ( x2<=3.)
00482         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));
00483       else
00484         {if ( x2<=11.) electrIonizationEnergyl2 =2.*std::exp(-2.*x2)/std::pow(x2,1.6);  }
00485     }
00486 
00487   G4double hFunctionl2 =(electrIonizationEnergyl2*2.*nl)/(tetal2*std::pow(velocityl2,3)); //takes into account  the polarization effect
00488 
00489    if (verboseLevel>0) G4cout << "  hFunctionl2=" << hFunctionl2<< G4endl;
00490 
00491   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.);
00492   //takes into account the reduced binding effect
00493 
00494     if (verboseLevel>0) G4cout << "  gFunctionl2=" << gFunctionl2<< G4endl;
00495 
00496   G4double sigmaPSS_l2 = 1.+(((2.*zIncident)/(screenedzTarget*tetal2))*(gFunctionl2-hFunctionl2)); //Binding-polarization factor
00497 
00498     if (verboseLevel>0) G4cout << "  sigmaPSS_l2=" << sigmaPSS_l2<< G4endl;
00499 
00500   const G4double cNaturalUnit= 137.;
00501      
00502   G4double yl2Formula=0.15*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(velocityl2/sigmaPSS_l2);
00503   
00504   G4double l2relativityCorrection = std::pow((1.+(1.1*yl2Formula*yl2Formula)),0.5)+yl2Formula; // Relativistic correction parameter
00505     
00506   G4double L2etaOverTheta2;
00507 
00508   G4double  universalFunction_l2 = 0.;
00509 
00510   G4double sigmaPSSR_l2 ;
00511  
00512   if ( velocityl2 < 20. )
00513   {
00514 
00515     L2etaOverTheta2 = (reducedEnergy*l2relativityCorrection)/((sigmaPSS_l2*tetal2)*(sigmaPSS_l2*tetal2));
00516 
00517     if ( (tetal2*sigmaPSS_l2>=0.2) && (tetal2*sigmaPSS_l2<=2.6670) && (L2etaOverTheta2>=0.1e-3) && (L2etaOverTheta2<=0.866e2) )
00518 
00519       universalFunction_l2 = FunctionFL2((tetal2*sigmaPSS_l2),L2etaOverTheta2);
00520 
00521     sigmaPSSR_l2 = (sigma0/(tetal2*sigmaPSS_l2))*universalFunction_l2;
00522 
00523     if (verboseLevel>0) G4cout << "  sigma PWBA L2 CS at low velocity range = " << sigmaPSSR_l2<< G4endl;
00524   }    
00525   else
00526   { 
00527     
00528     L2etaOverTheta2 = reducedEnergy /(tetal2*tetal2);
00529 
00530     if ( (tetal2>=0.2) && (tetal2<=2.6670) && (L2etaOverTheta2>=0.1e-3) && (L2etaOverTheta2<=0.866e2) )
00531 
00532       universalFunction_l2 = FunctionFL2((tetal2),L2etaOverTheta2);
00533    
00534     sigmaPSSR_l2 = (sigma0/tetal2)*universalFunction_l2;
00535 
00536       if (verboseLevel>0) G4cout << "  sigma PWBA L2 CS at medium and high velocity range = " << sigmaPSSR_l2<< G4endl;
00537 
00538   }
00539  
00540   G4double pssDeltal2 = (4./(systemMass*sigmaPSS_l2*tetal2))*(sigmaPSS_l2/velocityl2)*(sigmaPSS_l2/velocityl2);
00541 
00542   if (pssDeltal2>1) return 0.;
00543 
00544   G4double energyLossl2 = std::pow(1-pssDeltal2,0.5);
00545  
00546     if (verboseLevel>0) G4cout << "  energyLossl2=" << energyLossl2<< G4endl;
00547 
00548   G4double coulombDeflectionl2 
00549     =(8.*pi*zIncident/systemMass)*std::pow(tetal2*sigmaPSS_l2,-2.)*std::pow(velocityl2/sigmaPSS_l2,-3.)*(zTarget/screenedzTarget);
00550 
00551   G4double cParameterl2 = 2.*coulombDeflectionl2/(energyLossl2*(energyLossl2+1.));
00552 
00553   G4double coulombDeflectionFunction_l2 = 11.*ExpIntFunction(12,cParameterl2); //Coulomb-deflection effect correction
00554   // *** see Brandt, Phys Rev A10, p477, f25
00555   
00556     if (verboseLevel>0) G4cout << "  coulombDeflectionFunction_l2 =" << coulombDeflectionFunction_l2 << G4endl;
00557 
00558   G4double crossSection_L2 = coulombDeflectionFunction_l2 * sigmaPSSR_l2;
00559   //ECPSSR L2 -subshell cross section is estimated at perturbed-stationnairy-state(PSS)
00560   //and reduced by the energy-loss(E),the Coulomb deflection(C),and the relativity(R) effects
00561 
00562     if (verboseLevel>0) G4cout << "  crossSection_L2 =" << crossSection_L2 << G4endl;
00563 
00564   if (crossSection_L2 >= 0) 
00565   {
00566      return crossSection_L2 * barn;
00567   }
00568   else {return 0;}
00569 }

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

Implements G4VecpssrLiModel.

Definition at line 574 of file G4ecpssrBaseLixsModel.cc.

References G4Alpha::Alpha(), G4AtomicShell::BindingEnergy(), CalculateVelocity(), ExpIntFunction(), G4cout, G4endl, G4NistManager::GetAtomicMassAmu(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), G4AtomicTransitionManager::Instance(), G4NistManager::Instance(), G4INCL::Math::pi, G4Proton::Proton(), and G4AtomicTransitionManager::Shell().

00576 {
00577   if (zTarget <=13) return 0.;
00578 
00579   //this L3-CrossSection calculation method is done according to Werner Brandt and Grzegorz Lapicki, Phys.Rev.A20 N2 (1979),
00580   //and using data tables of O. Benka et al. At.Data Nucl.Data Tables Vol.22 No.3 (1978).
00581 
00582   G4NistManager* massManager = G4NistManager::Instance();
00583 
00584   G4AtomicTransitionManager*  transitionManager =  G4AtomicTransitionManager::Instance();
00585 
00586   G4double  zIncident = 0;
00587 
00588   G4Proton* aProtone = G4Proton::Proton();
00589   G4Alpha* aAlpha = G4Alpha::Alpha();
00590 
00591   if (massIncident == aProtone->GetPDGMass() )
00592 
00593    zIncident = (aProtone->GetPDGCharge())/eplus;
00594 
00595   else
00596     {
00597       if (massIncident == aAlpha->GetPDGMass())
00598 
00599           zIncident  = (aAlpha->GetPDGCharge())/eplus;
00600 
00601       else
00602         {
00603           G4cout << "*** WARNING in G4ecpssrBaseLixsModel::CalculateL3CrossSection : Proton or Alpha incident particles only. " << G4endl;
00604           G4cout << massIncident << ", " << aAlpha->GetPDGMass() << " (alpha)" << aProtone->GetPDGMass() << " (proton)" << G4endl;
00605           return 0;
00606         }
00607     }
00608 
00609   G4double l3BindingEnergy = transitionManager->Shell(zTarget,3)->BindingEnergy();
00610 
00611   G4double massTarget = (massManager->GetAtomicMassAmu(zTarget))*amu_c2;
00612 
00613   G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2;//Mass of the system (projectile, target)
00614 
00615   const G4double zlshell= 4.15;
00616 
00617   G4double screenedzTarget = zTarget-zlshell;//Effective nuclear charge as seen by electrons in L3-subshell
00618 
00619   const G4double rydbergMeV= 13.6056923e-6;
00620 
00621   const G4double nl= 2.;
00622 
00623   G4double tetal3 = (l3BindingEnergy*nl*nl)/((screenedzTarget*screenedzTarget)*rydbergMeV);//Screening parameter
00624 
00625     if (verboseLevel>0) G4cout << "  tetal3=" <<  tetal3<< G4endl;
00626 
00627   G4double reducedEnergy = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
00628 
00629   const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/barn ;//Bohr radius of hydrogen
00630 
00631   G4double sigma0 = 8.*pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.);
00632 
00633   G4double velocityl3 = CalculateVelocity(3, zTarget, massIncident, energyIncident);// Scaled velocity
00634 
00635     if (verboseLevel>0) G4cout << "  velocityl3=" << velocityl3<< G4endl;
00636 
00637   const G4double l23AnalyticalApproximation= 1.25;
00638 
00639   G4double x3 = (nl*l23AnalyticalApproximation)/velocityl3;
00640 
00641     if (verboseLevel>0) G4cout << "  x3=" << x3<< G4endl;
00642 
00643   G4double electrIonizationEnergyl3=0.;
00644 
00645   if ( x3<=0.035)  electrIonizationEnergyl3= 0.75*pi*(std::log(1./(x3*x3))-1.);
00646     else
00647     {
00648       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));
00649       else
00650         {
00651           if ( x3<=11.) electrIonizationEnergyl3 =2.*std::exp(-2.*x3)/std::pow(x3,1.6);}
00652     }
00653 
00654   G4double hFunctionl3 =(electrIonizationEnergyl3*2.*nl)/(tetal3*std::pow(velocityl3,3));//takes into account the polarization effect
00655 
00656     if (verboseLevel>0) G4cout << "  hFunctionl3=" << hFunctionl3<< G4endl;
00657 
00658   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.);
00659   //takes into account the reduced binding effect
00660 
00661     if (verboseLevel>0) G4cout << "  gFunctionl3=" << gFunctionl3<< G4endl;
00662 
00663   G4double sigmaPSS_l3 = 1.+(((2.*zIncident)/(screenedzTarget*tetal3))*(gFunctionl3-hFunctionl3));//Binding-polarization factor
00664 
00665     if (verboseLevel>0) G4cout << "sigmaPSS_l3 =" << sigmaPSS_l3<< G4endl;
00666 
00667   const G4double cNaturalUnit= 137.;
00668         
00669   G4double yl3Formula=0.15*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(velocityl3/sigmaPSS_l3);
00670         
00671   G4double l3relativityCorrection = std::pow((1.+(1.1*yl3Formula*yl3Formula)),0.5)+yl3Formula; // Relativistic correction parameter
00672   
00673   G4double L3etaOverTheta2;
00674  
00675   G4double  universalFunction_l3 = 0.;
00676  
00677   G4double sigmaPSSR_l3;
00678 
00679   if ( velocityl3 < 20. )
00680   {
00681 
00682     L3etaOverTheta2 = (reducedEnergy* l3relativityCorrection)/((sigmaPSS_l3*tetal3)*(sigmaPSS_l3*tetal3));
00683 
00684     if ( (tetal3*sigmaPSS_l3>=0.2) && (tetal3*sigmaPSS_l3<=2.6670) && (L3etaOverTheta2>=0.1e-3) && (L3etaOverTheta2<=0.866e2) )
00685 
00686       universalFunction_l3 = 2.*FunctionFL2((tetal3*sigmaPSS_l3), L3etaOverTheta2  );
00687 
00688     sigmaPSSR_l3 = (sigma0/(tetal3*sigmaPSS_l3))*universalFunction_l3;
00689 
00690     if (verboseLevel>0) G4cout << "  sigma PWBA L3 CS at low velocity range = " << sigmaPSSR_l3<< G4endl;
00691 
00692   }
00693 
00694   else 
00695 
00696   {
00697 
00698     L3etaOverTheta2 = reducedEnergy/(tetal3*tetal3);
00699 
00700     if ( (tetal3>=0.2) && (tetal3<=2.6670) && (L3etaOverTheta2>=0.1e-3) && (L3etaOverTheta2<=0.866e2) ) 
00701      
00702       universalFunction_l3 = 2.*FunctionFL2(tetal3, L3etaOverTheta2  );
00703 
00704     sigmaPSSR_l3 = (sigma0/tetal3)*universalFunction_l3;
00705 
00706       if (verboseLevel>0) G4cout << "  sigma PWBA L3 CS at medium and high velocity range = " << sigmaPSSR_l3<< G4endl;
00707   }
00708   
00709   G4double pssDeltal3 = (4./(systemMass*sigmaPSS_l3*tetal3))*(sigmaPSS_l3/velocityl3)*(sigmaPSS_l3/velocityl3);
00710 
00711     if (verboseLevel>0) G4cout << "  pssDeltal3=" << pssDeltal3<< G4endl;
00712 
00713   if (pssDeltal3>1) return 0.;
00714 
00715   G4double energyLossl3 = std::pow(1-pssDeltal3,0.5);
00716 
00717   if (verboseLevel>0) G4cout << "  energyLossl3=" << energyLossl3<< G4endl;
00718 
00719   G4double coulombDeflectionl3 = 
00720     (8.*pi*zIncident/systemMass)*std::pow(tetal3*sigmaPSS_l3,-2.)*std::pow(velocityl3/sigmaPSS_l3,-3.)*(zTarget/screenedzTarget);
00721 
00722   G4double cParameterl3 = 2.*coulombDeflectionl3/(energyLossl3*(energyLossl3+1.));
00723 
00724   G4double coulombDeflectionFunction_l3 = 11.*ExpIntFunction(12,cParameterl3);//Coulomb-deflection effect correction
00725   // *** see Brandt, Phys Rev A10, p477, f25
00726 
00727     if (verboseLevel>0) G4cout << "  coulombDeflectionFunction_l3 =" << coulombDeflectionFunction_l3 << G4endl;
00728 
00729   G4double crossSection_L3 =  coulombDeflectionFunction_l3 * sigmaPSSR_l3;
00730   //ECPSSR L3 -subshell cross section is estimated at perturbed-stationnairy-state(PSS)
00731   //and reduced by the energy-loss(E),the Coulomb deflection(C),and the relativity(R) effects
00732 
00733     if (verboseLevel>0) G4cout << "  crossSection_L3 =" << crossSection_L3 << G4endl;
00734  
00735   if (crossSection_L3 >= 0) 
00736   {
00737     return crossSection_L3 * barn;
00738   }
00739   else {return 0;}
00740 }

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

Definition at line 744 of file G4ecpssrBaseLixsModel.cc.

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

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

00746 {
00747 
00748   G4AtomicTransitionManager*  transitionManager =  G4AtomicTransitionManager::Instance();
00749 
00750   G4double liBindingEnergy = transitionManager->Shell(zTarget,subShell)->BindingEnergy();
00751 
00752   G4Proton* aProtone = G4Proton::Proton();
00753   G4Alpha* aAlpha = G4Alpha::Alpha();
00754 
00755   if (!((massIncident == aProtone->GetPDGMass()) || (massIncident == aAlpha->GetPDGMass())))
00756     {
00757       G4cout << "*** WARNING in G4ecpssrBaseLixsModel::CalculateVelocity : Proton or Alpha incident particles only. " << G4endl;
00758       G4cout << massIncident << ", " << aAlpha->GetPDGMass() << " (alpha)" << aProtone->GetPDGMass() << " (proton)" << G4endl;
00759       return 0;
00760     }
00761 
00762   const G4double zlshell= 4.15;
00763 
00764   G4double screenedzTarget = zTarget- zlshell;
00765 
00766   const G4double rydbergMeV= 13.6056923e-6;
00767 
00768   const G4double nl= 2.;
00769 
00770   G4double tetali = (liBindingEnergy*nl*nl)/(screenedzTarget*screenedzTarget*rydbergMeV);
00771 
00772   G4double reducedEnergy = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
00773 
00774   G4double velocity = 2.*nl*std::pow(reducedEnergy,0.5)/tetali;
00775   // *** see Brandt, Phys Rev A10, p10, f4
00776   
00777   return velocity;
00778 }

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

Definition at line 122 of file G4ecpssrBaseLixsModel.cc.

References G4cout, and G4endl.

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

00124 {
00125 // this function allows fast evaluation of the n order exponential integral function En(x)
00126 
00127   G4int i;
00128   G4int ii;
00129   G4int nm1;
00130   G4double a;
00131   G4double b;
00132   G4double c;
00133   G4double d;
00134   G4double del;
00135   G4double fact;
00136   G4double h;
00137   G4double psi;
00138   G4double ans = 0;
00139   const G4double euler= 0.5772156649;
00140   const G4int maxit= 100;
00141   const G4double fpmin = 1.0e-30;
00142   const G4double eps = 1.0e-7;
00143   nm1=n-1;
00144   if (n<0 || x<0.0 || (x==0.0 && (n==0 || n==1)))
00145   G4cout << "*** WARNING in G4ecpssrBaseLixsModel::ExpIntFunction: bad arguments in ExpIntFunction" 
00146          << G4endl;
00147   else {
00148        if (n==0) ans=std::exp(-x)/x;
00149         else {
00150            if (x==0.0) ans=1.0/nm1;
00151               else {
00152                    if (x > 1.0) {
00153                                 b=x+n;
00154                                 c=1.0/fpmin;
00155                                 d=1.0/b;
00156                                 h=d;
00157                                 for (i=1;i<=maxit;i++) {
00158                                   a=-i*(nm1+i);
00159                                   b +=2.0;
00160                                   d=1.0/(a*d+b);
00161                                   c=b+a/c;
00162                                   del=c*d;
00163                                   h *=del;
00164                                       if (std::fabs(del-1.0) < eps) {
00165                                         ans=h*std::exp(-x);
00166                                         return ans;
00167                                       }
00168                                 }
00169                    } else {
00170                      ans = (nm1!=0 ? 1.0/nm1 : -std::log(x)-euler);
00171                      fact=1.0;
00172                      for (i=1;i<=maxit;i++) {
00173                        fact *=-x/i;
00174                        if (i !=nm1) del = -fact/(i-nm1);
00175                        else {
00176                          psi = -euler;
00177                          for (ii=1;ii<=nm1;ii++) psi +=1.0/ii;
00178                          del=fact*(-std::log(x)+psi);
00179                        }
00180                        ans += del;
00181                        if (std::fabs(del) < std::fabs(ans)*eps) return ans;
00182                      }
00183                    }
00184               }
00185         }
00186   }
00187 return ans;
00188 }


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:51:49 2013 for Geant4 by  doxygen 1.4.7