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