49   char *path = getenv(
"G4LEDATA");
 
   52     G4Exception(
"G4ecpssrLCrossSection::G4ecpssrBaseLixsModel()",
"em0006", 
FatalException ,
"G4LEDATA environment variable not set");
 
   55   std::ostringstream fileName1;
 
   56   std::ostringstream fileName2;
 
   58   fileName1 << path << 
"/pixe/uf/FL1.dat";
 
   59   fileName2 << path << 
"/pixe/uf/FL2.dat";
 
   63   std::ifstream FL1(fileName1.str().c_str());
 
   64   if (!FL1) 
G4Exception(
"G4ecpssrLCrossSection::G4ecpssrBaseLixsModel()",
"em0003",
FatalException, 
"error opening FL1 data file");
 
   89   std::ifstream FL2(fileName2.str().c_str());
 
   90   if (!FL2) 
G4Exception(
"G4ecpssrLCrossSection::G4ecpssrBaseLixsModel()",
"em0003", 
FatalException,
" error opening FL2 data file");
 
  138   static const G4double euler= 0.5772156649;
 
  139   static const G4int maxit= 100;
 
  140   static const G4double fpmin = 1.0e-30;
 
  143   if (n<0 || x<0.0 || (x==0.0 && (n==0 || n==1)))
 
  144   G4cout << 
"*** WARNING in G4ecpssrBaseLixsModel::ExpIntFunction: bad arguments in ExpIntFunction" 
  147     if (n==0) ans=
G4Exp(-x)/
x;
 
  149       if (x==0.0) ans=1.0/nm1;
 
  156           for (i=1;i<=maxit;i++) {
 
  163             if (std::fabs(del-1.0) < 
eps) {
 
  169           ans = (nm1!=0 ? 1.0/nm1 : -std::log(x)-euler);
 
  171           for (i=1;i<=maxit;i++) {
 
  173             if (i !=nm1) del = -fact/(i-nm1);
 
  176               for (ii=1;ii<=nm1;ii++) psi +=1.0/ii;
 
  177               del=fact*(-std::log(x)+psi);
 
  180             if (std::fabs(del) < std::fabs(ans)*
eps) 
return ans;
 
  194   if (zTarget <=4) 
return 0.;
 
  219           G4cout << 
"*** WARNING in G4ecpssrBaseLixsModel::CalculateL1CrossSection : Proton or Alpha incident particles only. " << 
G4endl;
 
  227   G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2; 
 
  228   static const G4double zlshell= 4.15;
 
  230   G4double screenedzTarget = zTarget-zlshell; 
 
  231   static const G4double rydbergMeV= 13.6056923e-6;
 
  235   G4double tetal1 = (l1BindingEnergy*nl*nl)/((screenedzTarget*screenedzTarget)*rydbergMeV); 
 
  240   G4double reducedEnergy = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
 
  244   static const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/
barn ; 
 
  246   G4double sigma0 = 8.*
pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.);
 
  254   static const G4double l1AnalyticalApproximation= 1.5;
 
  255   G4double x1 =(nl*l1AnalyticalApproximation)/velocityl1;
 
  261   G4double electrIonizationEnergyl1=0.;
 
  266   if ( x1<=0.035)  electrIonizationEnergyl1= 0.75*
pi*(std::log(1./(x1*x1))-1.);
 
  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));
 
  272         {
if ( x1<=11.) electrIonizationEnergyl1 =2.*
G4Exp(-2.*x1)/std::pow(x1,1.6);}
 
  275   G4double hFunctionl1 =(electrIonizationEnergyl1*2.*nl)/(tetal1*std::pow(velocityl1,3)); 
 
  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.);
 
  285   G4double sigmaPSS_l1 = 1.+(((2.*zIncident)/(screenedzTarget*tetal1))*(gFunctionl1-hFunctionl1)); 
 
  294   G4double yl1Formula=0.4*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(nl*velocityl1/sigmaPSS_l1);
 
  299   G4double l1relativityCorrection = std::pow((1.+(1.1*yl1Formula*yl1Formula)),0.5)+yl1Formula; 
 
  313   if ( velocityl1 <20.  )
 
  316     L1etaOverTheta2 =(reducedEnergy* l1relativityCorrection)/((tetal1*sigmaPSS_l1)*(tetal1*sigmaPSS_l1));
 
  322     if ( ((tetal1*sigmaPSS_l1) >=0.2) && ((tetal1*sigmaPSS_l1) <=2.6670) && (L1etaOverTheta2>=0.1e-3) && (L1etaOverTheta2<=0.866
e2) )
 
  323       universalFunction_l1 = 
FunctionFL1((tetal1*sigmaPSS_l1),L1etaOverTheta2);
 
  327     sigmaPSSR_l1 = (sigma0/(tetal1*sigmaPSS_l1))*universalFunction_l1;
 
  334     L1etaOverTheta2 = reducedEnergy/(tetal1*tetal1);
 
  340     if ( (tetal1 >=0.2) && (tetal1 <=2.6670) && (L1etaOverTheta2>=0.1e-3) && (L1etaOverTheta2<=0.866
e2) )
 
  341       universalFunction_l1 = 
FunctionFL1(tetal1,L1etaOverTheta2);
 
  343     if (
verboseLevel>0) 
G4cout << 
"at medium and high velocity range, universalFunction_l1  =" << universalFunction_l1 << 
G4endl;
 
  345     sigmaPSSR_l1 = (sigma0/tetal1)*universalFunction_l1;
 
  351   G4double pssDeltal1 = (4./(systemMass *sigmaPSS_l1*tetal1))*(sigmaPSS_l1/velocityl1)*(sigmaPSS_l1/velocityl1);
 
  357   if (pssDeltal1>1) 
return 0.;
 
  359   G4double energyLossl1 = std::pow(1-pssDeltal1,0.5);
 
  366    (8.*
pi*zIncident/systemMass)*std::pow(tetal1*sigmaPSS_l1,-2.)*std::pow(velocityl1/sigmaPSS_l1,-3.)*(zTarget/screenedzTarget);
 
  370   G4double cParameterl1 =2.* coulombDeflectionl1/(energyLossl1*(energyLossl1+1.));
 
  377   G4double crossSection_L1 = coulombDeflectionFunction_l1 * sigmaPSSR_l1;
 
  384   if (crossSection_L1 >= 0)
 
  386     return crossSection_L1 * 
barn;
 
  397   if (zTarget <=13 ) 
return 0.;
 
  421           G4cout << 
"*** WARNING in G4ecpssrBaseLixsModel::CalculateL2CrossSection : Proton or Alpha incident particles only. " << 
G4endl;
 
  431   G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2; 
 
  435   G4double screenedzTarget = zTarget-zlshell; 
 
  437   const G4double rydbergMeV= 13.6056923e-6;
 
  441   G4double tetal2 = (l2BindingEnergy*nl*nl)/((screenedzTarget*screenedzTarget)*rydbergMeV); 
 
  445   G4double reducedEnergy = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
 
  447   const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/
barn ; 
 
  449   G4double sigma0 = 8.*
pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.);
 
  455   const G4double l23AnalyticalApproximation= 1.25;
 
  457   G4double x2 = (nl*l23AnalyticalApproximation)/velocityl2;
 
  461   G4double electrIonizationEnergyl2=0.;
 
  463   if ( x2<=0.035)  electrIonizationEnergyl2= 0.75*
pi*(std::log(1./(x2*x2))-1.);
 
  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));
 
  469         {
if ( x2<=11.) electrIonizationEnergyl2 =2.*
G4Exp(-2.*x2)/std::pow(x2,1.6);  }
 
  472   G4double hFunctionl2 =(electrIonizationEnergyl2*2.*nl)/(tetal2*std::pow(velocityl2,3)); 
 
  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.);
 
  481   G4double sigmaPSS_l2 = 1.+(((2.*zIncident)/(screenedzTarget*tetal2))*(gFunctionl2-hFunctionl2)); 
 
  487   G4double yl2Formula=0.15*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(velocityl2/sigmaPSS_l2);
 
  489   G4double l2relativityCorrection = std::pow((1.+(1.1*yl2Formula*yl2Formula)),0.5)+yl2Formula; 
 
  497   if ( velocityl2 < 20. )
 
  500     L2etaOverTheta2 = (reducedEnergy*l2relativityCorrection)/((sigmaPSS_l2*tetal2)*(sigmaPSS_l2*tetal2));
 
  502     if ( (tetal2*sigmaPSS_l2>=0.2) && (tetal2*sigmaPSS_l2<=2.6670) && (L2etaOverTheta2>=0.1e-3) && (L2etaOverTheta2<=0.866
e2) )
 
  503       universalFunction_l2 = 
FunctionFL2((tetal2*sigmaPSS_l2),L2etaOverTheta2);
 
  505     sigmaPSSR_l2 = (sigma0/(tetal2*sigmaPSS_l2))*universalFunction_l2;
 
  512     L2etaOverTheta2 = reducedEnergy /(tetal2*tetal2);
 
  514     if ( (tetal2>=0.2) && (tetal2<=2.6670) && (L2etaOverTheta2>=0.1e-3) && (L2etaOverTheta2<=0.866
e2) )
 
  515       universalFunction_l2 = 
FunctionFL2((tetal2),L2etaOverTheta2);
 
  517     sigmaPSSR_l2 = (sigma0/tetal2)*universalFunction_l2;
 
  523   G4double pssDeltal2 = (4./(systemMass*sigmaPSS_l2*tetal2))*(sigmaPSS_l2/velocityl2)*(sigmaPSS_l2/velocityl2);
 
  525   if (pssDeltal2>1) 
return 0.;
 
  527   G4double energyLossl2 = std::pow(1-pssDeltal2,0.5);
 
  532     =(8.*
pi*zIncident/systemMass)*std::pow(tetal2*sigmaPSS_l2,-2.)*std::pow(velocityl2/sigmaPSS_l2,-3.)*(zTarget/screenedzTarget);
 
  534   G4double cParameterl2 = 2.*coulombDeflectionl2/(energyLossl2*(energyLossl2+1.));
 
  541   G4double crossSection_L2 = coulombDeflectionFunction_l2 * sigmaPSSR_l2;
 
  547   if (crossSection_L2 >= 0)
 
  549      return crossSection_L2 * 
barn;
 
  560   if (zTarget <=13) 
return 0.;
 
  586           G4cout << 
"*** WARNING in G4ecpssrBaseLixsModel::CalculateL3CrossSection : Proton or Alpha incident particles only. " << 
G4endl;
 
  596   G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2;
 
  600   G4double screenedzTarget = zTarget-zlshell;
 
  602   const G4double rydbergMeV= 13.6056923e-6;
 
  606   G4double tetal3 = (l3BindingEnergy*nl*nl)/((screenedzTarget*screenedzTarget)*rydbergMeV);
 
  610   G4double reducedEnergy = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
 
  612   const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/
barn ;
 
  614   G4double sigma0 = 8.*
pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.);
 
  620   const G4double l23AnalyticalApproximation= 1.25;
 
  622   G4double x3 = (nl*l23AnalyticalApproximation)/velocityl3;
 
  626   G4double electrIonizationEnergyl3=0.;
 
  628   if ( x3<=0.035)  electrIonizationEnergyl3= 0.75*
pi*(std::log(1./(x3*x3))-1.);
 
  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));
 
  634           if ( x3<=11.) electrIonizationEnergyl3 =2.*
G4Exp(-2.*x3)/std::pow(x3,1.6);}
 
  637   G4double hFunctionl3 =(electrIonizationEnergyl3*2.*nl)/(tetal3*std::pow(velocityl3,3));
 
  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.);
 
  646   G4double sigmaPSS_l3 = 1.+(((2.*zIncident)/(screenedzTarget*tetal3))*(gFunctionl3-hFunctionl3));
 
  652   G4double yl3Formula=0.15*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(velocityl3/sigmaPSS_l3);
 
  654   G4double l3relativityCorrection = std::pow((1.+(1.1*yl3Formula*yl3Formula)),0.5)+yl3Formula; 
 
  662   if ( velocityl3 < 20. )
 
  665     L3etaOverTheta2 = (reducedEnergy* l3relativityCorrection)/((sigmaPSS_l3*tetal3)*(sigmaPSS_l3*tetal3));
 
  667     if ( (tetal3*sigmaPSS_l3>=0.2) && (tetal3*sigmaPSS_l3<=2.6670) && (L3etaOverTheta2>=0.1e-3) && (L3etaOverTheta2<=0.866
e2) )
 
  669       universalFunction_l3 = 2.*
FunctionFL2((tetal3*sigmaPSS_l3), L3etaOverTheta2  );
 
  671     sigmaPSSR_l3 = (sigma0/(tetal3*sigmaPSS_l3))*universalFunction_l3;
 
  681     L3etaOverTheta2 = reducedEnergy/(tetal3*tetal3);
 
  683     if ( (tetal3>=0.2) && (tetal3<=2.6670) && (L3etaOverTheta2>=0.1e-3) && (L3etaOverTheta2<=0.866
e2) )
 
  685       universalFunction_l3 = 2.*
FunctionFL2(tetal3, L3etaOverTheta2  );
 
  687     sigmaPSSR_l3 = (sigma0/tetal3)*universalFunction_l3;
 
  692   G4double pssDeltal3 = (4./(systemMass*sigmaPSS_l3*tetal3))*(sigmaPSS_l3/velocityl3)*(sigmaPSS_l3/velocityl3);
 
  696   if (pssDeltal3>1) 
return 0.;
 
  698   G4double energyLossl3 = std::pow(1-pssDeltal3,0.5);
 
  703     (8.*
pi*zIncident/systemMass)*std::pow(tetal3*sigmaPSS_l3,-2.)*std::pow(velocityl3/sigmaPSS_l3,-3.)*(zTarget/screenedzTarget);
 
  705   G4double cParameterl3 = 2.*coulombDeflectionl3/(energyLossl3*(energyLossl3+1.));
 
  712   G4double crossSection_L3 =  coulombDeflectionFunction_l3 * sigmaPSSR_l3;
 
  718   if (crossSection_L3 >= 0)
 
  720     return crossSection_L3 * 
barn;
 
  740       G4cout << 
"*** WARNING in G4ecpssrBaseLixsModel::CalculateVelocity : Proton or Alpha incident particles only. " << 
G4endl;
 
  747   G4double screenedzTarget = zTarget- zlshell;
 
  749   const G4double rydbergMeV= 13.6056923e-6;
 
  753   G4double tetali = (liBindingEnergy*nl*nl)/(screenedzTarget*screenedzTarget*rydbergMeV);
 
  755   G4double reducedEnergy = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
 
  757   G4double velocity = 2.*nl*std::pow(reducedEnergy,0.5)/tetali;
 
  802     std::vector<double>::iterator t2 = std::upper_bound(
dummyVec1.begin(),
dummyVec1.end(), k);
 
  803     std::vector<double>::iterator t1 = t2-1;
 
  805     std::vector<double>::iterator e12 = std::upper_bound(
aVecMap1[(*t1)].begin(),
aVecMap1[(*t1)].end(), theta);
 
  806     std::vector<double>::iterator e11 = e12-1;
 
  808     std::vector<double>::iterator e22 = std::upper_bound(
aVecMap1[(*t2)].begin(),
aVecMap1[(*t2)].end(), theta);
 
  809     std::vector<double>::iterator e21 = e22-1;
 
  818     xs11 = 
FL1Data[valueT1][valueE11];
 
  819     xs12 = 
FL1Data[valueT1][valueE12];
 
  820     xs21 = 
FL1Data[valueT2][valueE21];
 
  821     xs22 = 
FL1Data[valueT2][valueE22];
 
  837   G4double xsProduct = xs11 * xs12 * xs21 * xs22;
 
  839   if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) 
return (0.);
 
  893     std::vector<double>::iterator t2 = std::upper_bound(
dummyVec2.begin(),
dummyVec2.end(), k);
 
  894     std::vector<double>::iterator t1 = t2-1;
 
  896     std::vector<double>::iterator e12 = std::upper_bound(
aVecMap2[(*t1)].begin(),
aVecMap2[(*t1)].end(), theta);
 
  897     std::vector<double>::iterator e11 = e12-1;
 
  899     std::vector<double>::iterator e22 = std::upper_bound(
aVecMap2[(*t2)].begin(),
aVecMap2[(*t2)].end(), theta);
 
  900     std::vector<double>::iterator e21 = e22-1;
 
  909     xs11 = 
FL2Data[valueT1][valueE11];
 
  910     xs12 = 
FL2Data[valueT1][valueE12];
 
  911     xs21 = 
FL2Data[valueT2][valueE21];
 
  912     xs22 = 
FL2Data[valueT2][valueE22];
 
  928   G4double xsProduct = xs11 * xs12 * xs21 * xs22;
 
  930   if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) 
return (0.);
 
  953   G4double value = xs1 + (xs2 - xs1)*(e - e1)/ (e2 - 
e1);
 
  979   G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
 
  980   G4double b = std::log10(xs2) - a*std::log10(e2);
 
  981   G4double sigma = a*std::log10(e) + b;
 
  982   G4double value = (std::pow(10.,sigma));
 
G4double LinLinInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
 
G4double ExpIntFunction(G4int n, G4double x)
 
G4double FunctionFL1(G4double k, G4double theta)
 
std::vector< double > dummyVec1
 
G4double FunctionFL2(G4double k, G4double theta)
 
G4double QuadInterpolator(G4double e11, G4double e12, G4double e21, G4double e22, G4double x11, G4double x12, G4double x21, G4double x22, G4double t1, G4double t2, G4double t, G4double e)
 
static const G4double eps
 
G4double BindingEnergy() const 
 
G4double CalculateVelocity(G4int subShell, G4int zTarget, G4double massIncident, G4double energyIncident)
 
static G4NistManager * Instance()
 
G4double CalculateL2CrossSection(G4int zTarget, G4double massIncident, G4double energyIncident)
 
G4double CalculateL1CrossSection(G4int zTarget, G4double massIncident, G4double energyIncident)
 
G4GLOB_DLL std::ostream G4cout
 
G4double CalculateL3CrossSection(G4int zTarget, G4double massIncident, G4double energyIncident)
 
static G4Proton * Proton()
 
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
 
std::vector< double > dummyVec2
 
G4double G4Exp(G4double initial_x)
Exponential Function double precision. 
 
G4double LogLogInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
 
G4double GetPDGMass() const 
 
const G4double x[NPOINTSGL]
 
G4double GetAtomicMassAmu(const G4String &symb) const 
 
G4double LinLogInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
 
static const double eplus
 
G4double GetPDGCharge() const 
 
static G4AtomicTransitionManager * Instance()
 
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const