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");
140 const G4int maxit= 100;
144 if (n<0 || x<0.0 || (x==0.0 && (n==0 || n==1)))
145 G4cout <<
"*** WARNING in G4ecpssrBaseLixsModel::ExpIntFunction: bad arguments in ExpIntFunction"
148 if (n==0) ans=std::exp(-x)/
x;
150 if (x==0.0) ans=1.0/nm1;
157 for (i=1;i<=maxit;i++) {
164 if (std::fabs(del-1.0) <
eps) {
170 ans = (nm1!=0 ? 1.0/nm1 : -std::log(x)-euler);
172 for (i=1;i<=maxit;i++) {
174 if (i !=nm1) del = -fact/(i-nm1);
177 for (ii=1;ii<=nm1;ii++) psi +=1.0/ii;
178 del=fact*(-std::log(x)+psi);
181 if (std::fabs(del) < std::fabs(ans)*
eps)
return ans;
195 if (zTarget <=4)
return 0.;
220 G4cout <<
"*** WARNING in G4ecpssrBaseLixsModel::CalculateL1CrossSection : Proton or Alpha incident particles only. " <<
G4endl;
230 G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2;
235 G4double screenedzTarget = zTarget-zlshell;
237 const G4double rydbergMeV= 13.6056923e-6;
242 G4double tetal1 = (l1BindingEnergy*nl*nl)/((screenedzTarget*screenedzTarget)*rydbergMeV);
247 G4double reducedEnergy = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
251 const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/
barn ;
253 G4double sigma0 = 8.*
pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.);
261 const G4double l1AnalyticalApproximation= 1.5;
262 G4double x1 =(nl*l1AnalyticalApproximation)/velocityl1;
268 G4double electrIonizationEnergyl1=0.;
273 if ( x1<=0.035) electrIonizationEnergyl1= 0.75*
pi*(std::log(1./(x1*x1))-1.);
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));
279 {
if ( x1<=11.) electrIonizationEnergyl1 =2.*std::exp(-2.*x1)/std::pow(x1,1.6);}
282 G4double hFunctionl1 =(electrIonizationEnergyl1*2.*nl)/(tetal1*std::pow(velocityl1,3));
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.);
292 G4double sigmaPSS_l1 = 1.+(((2.*zIncident)/(screenedzTarget*tetal1))*(gFunctionl1-hFunctionl1));
301 G4double yl1Formula=0.4*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(nl*velocityl1/sigmaPSS_l1);
306 G4double l1relativityCorrection = std::pow((1.+(1.1*yl1Formula*yl1Formula)),0.5)+yl1Formula;
320 if ( velocityl1 <20. )
323 L1etaOverTheta2 =(reducedEnergy* l1relativityCorrection)/((tetal1*sigmaPSS_l1)*(tetal1*sigmaPSS_l1));
329 if ( ((tetal1*sigmaPSS_l1) >=0.2) && ((tetal1*sigmaPSS_l1) <=2.6670) && (L1etaOverTheta2>=0.1e-3) && (L1etaOverTheta2<=0.866
e2) )
331 universalFunction_l1 =
FunctionFL1((tetal1*sigmaPSS_l1),L1etaOverTheta2);
335 sigmaPSSR_l1 = (sigma0/(tetal1*sigmaPSS_l1))*universalFunction_l1;
346 L1etaOverTheta2 = reducedEnergy/(tetal1*tetal1);
352 if ( (tetal1 >=0.2) && (tetal1 <=2.6670) && (L1etaOverTheta2>=0.1e-3) && (L1etaOverTheta2<=0.866
e2) )
354 universalFunction_l1 =
FunctionFL1(tetal1,L1etaOverTheta2);
356 if (
verboseLevel>0)
G4cout <<
"at medium and high velocity range, universalFunction_l1 =" << universalFunction_l1 <<
G4endl;
358 sigmaPSSR_l1 = (sigma0/tetal1)*universalFunction_l1;
364 G4double pssDeltal1 = (4./(systemMass *sigmaPSS_l1*tetal1))*(sigmaPSS_l1/velocityl1)*(sigmaPSS_l1/velocityl1);
370 if (pssDeltal1>1)
return 0.;
372 G4double energyLossl1 = std::pow(1-pssDeltal1,0.5);
379 (8.*
pi*zIncident/systemMass)*std::pow(tetal1*sigmaPSS_l1,-2.)*std::pow(velocityl1/sigmaPSS_l1,-3.)*(zTarget/screenedzTarget);
383 G4double cParameterl1 =2.* coulombDeflectionl1/(energyLossl1*(energyLossl1+1.));
390 G4double crossSection_L1 = coulombDeflectionFunction_l1 * sigmaPSSR_l1;
397 if (crossSection_L1 >= 0)
399 return crossSection_L1 *
barn;
410 if (zTarget <=13 )
return 0.;
436 G4cout <<
"*** WARNING in G4ecpssrBaseLixsModel::CalculateL2CrossSection : Proton or Alpha incident particles only. " <<
G4endl;
446 G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2;
450 G4double screenedzTarget = zTarget-zlshell;
452 const G4double rydbergMeV= 13.6056923e-6;
456 G4double tetal2 = (l2BindingEnergy*nl*nl)/((screenedzTarget*screenedzTarget)*rydbergMeV);
460 G4double reducedEnergy = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
462 const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/
barn ;
464 G4double sigma0 = 8.*
pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.);
470 const G4double l23AnalyticalApproximation= 1.25;
472 G4double x2 = (nl*l23AnalyticalApproximation)/velocityl2;
476 G4double electrIonizationEnergyl2=0.;
478 if ( x2<=0.035) electrIonizationEnergyl2= 0.75*
pi*(std::log(1./(x2*x2))-1.);
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));
484 {
if ( x2<=11.) electrIonizationEnergyl2 =2.*std::exp(-2.*x2)/std::pow(x2,1.6); }
487 G4double hFunctionl2 =(electrIonizationEnergyl2*2.*nl)/(tetal2*std::pow(velocityl2,3));
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.);
496 G4double sigmaPSS_l2 = 1.+(((2.*zIncident)/(screenedzTarget*tetal2))*(gFunctionl2-hFunctionl2));
502 G4double yl2Formula=0.15*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(velocityl2/sigmaPSS_l2);
504 G4double l2relativityCorrection = std::pow((1.+(1.1*yl2Formula*yl2Formula)),0.5)+yl2Formula;
512 if ( velocityl2 < 20. )
515 L2etaOverTheta2 = (reducedEnergy*l2relativityCorrection)/((sigmaPSS_l2*tetal2)*(sigmaPSS_l2*tetal2));
517 if ( (tetal2*sigmaPSS_l2>=0.2) && (tetal2*sigmaPSS_l2<=2.6670) && (L2etaOverTheta2>=0.1e-3) && (L2etaOverTheta2<=0.866
e2) )
519 universalFunction_l2 =
FunctionFL2((tetal2*sigmaPSS_l2),L2etaOverTheta2);
521 sigmaPSSR_l2 = (sigma0/(tetal2*sigmaPSS_l2))*universalFunction_l2;
528 L2etaOverTheta2 = reducedEnergy /(tetal2*tetal2);
530 if ( (tetal2>=0.2) && (tetal2<=2.6670) && (L2etaOverTheta2>=0.1e-3) && (L2etaOverTheta2<=0.866
e2) )
532 universalFunction_l2 =
FunctionFL2((tetal2),L2etaOverTheta2);
534 sigmaPSSR_l2 = (sigma0/tetal2)*universalFunction_l2;
540 G4double pssDeltal2 = (4./(systemMass*sigmaPSS_l2*tetal2))*(sigmaPSS_l2/velocityl2)*(sigmaPSS_l2/velocityl2);
542 if (pssDeltal2>1)
return 0.;
544 G4double energyLossl2 = std::pow(1-pssDeltal2,0.5);
549 =(8.*
pi*zIncident/systemMass)*std::pow(tetal2*sigmaPSS_l2,-2.)*std::pow(velocityl2/sigmaPSS_l2,-3.)*(zTarget/screenedzTarget);
551 G4double cParameterl2 = 2.*coulombDeflectionl2/(energyLossl2*(energyLossl2+1.));
558 G4double crossSection_L2 = coulombDeflectionFunction_l2 * sigmaPSSR_l2;
564 if (crossSection_L2 >= 0)
566 return crossSection_L2 *
barn;
577 if (zTarget <=13)
return 0.;
603 G4cout <<
"*** WARNING in G4ecpssrBaseLixsModel::CalculateL3CrossSection : Proton or Alpha incident particles only. " <<
G4endl;
613 G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2;
617 G4double screenedzTarget = zTarget-zlshell;
619 const G4double rydbergMeV= 13.6056923e-6;
623 G4double tetal3 = (l3BindingEnergy*nl*nl)/((screenedzTarget*screenedzTarget)*rydbergMeV);
627 G4double reducedEnergy = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
629 const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/
barn ;
631 G4double sigma0 = 8.*
pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.);
637 const G4double l23AnalyticalApproximation= 1.25;
639 G4double x3 = (nl*l23AnalyticalApproximation)/velocityl3;
643 G4double electrIonizationEnergyl3=0.;
645 if ( x3<=0.035) electrIonizationEnergyl3= 0.75*
pi*(std::log(1./(x3*x3))-1.);
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));
651 if ( x3<=11.) electrIonizationEnergyl3 =2.*std::exp(-2.*x3)/std::pow(x3,1.6);}
654 G4double hFunctionl3 =(electrIonizationEnergyl3*2.*nl)/(tetal3*std::pow(velocityl3,3));
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.);
663 G4double sigmaPSS_l3 = 1.+(((2.*zIncident)/(screenedzTarget*tetal3))*(gFunctionl3-hFunctionl3));
669 G4double yl3Formula=0.15*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(velocityl3/sigmaPSS_l3);
671 G4double l3relativityCorrection = std::pow((1.+(1.1*yl3Formula*yl3Formula)),0.5)+yl3Formula;
679 if ( velocityl3 < 20. )
682 L3etaOverTheta2 = (reducedEnergy* l3relativityCorrection)/((sigmaPSS_l3*tetal3)*(sigmaPSS_l3*tetal3));
684 if ( (tetal3*sigmaPSS_l3>=0.2) && (tetal3*sigmaPSS_l3<=2.6670) && (L3etaOverTheta2>=0.1e-3) && (L3etaOverTheta2<=0.866
e2) )
686 universalFunction_l3 = 2.*
FunctionFL2((tetal3*sigmaPSS_l3), L3etaOverTheta2 );
688 sigmaPSSR_l3 = (sigma0/(tetal3*sigmaPSS_l3))*universalFunction_l3;
698 L3etaOverTheta2 = reducedEnergy/(tetal3*tetal3);
700 if ( (tetal3>=0.2) && (tetal3<=2.6670) && (L3etaOverTheta2>=0.1e-3) && (L3etaOverTheta2<=0.866
e2) )
702 universalFunction_l3 = 2.*
FunctionFL2(tetal3, L3etaOverTheta2 );
704 sigmaPSSR_l3 = (sigma0/tetal3)*universalFunction_l3;
709 G4double pssDeltal3 = (4./(systemMass*sigmaPSS_l3*tetal3))*(sigmaPSS_l3/velocityl3)*(sigmaPSS_l3/velocityl3);
713 if (pssDeltal3>1)
return 0.;
715 G4double energyLossl3 = std::pow(1-pssDeltal3,0.5);
720 (8.*
pi*zIncident/systemMass)*std::pow(tetal3*sigmaPSS_l3,-2.)*std::pow(velocityl3/sigmaPSS_l3,-3.)*(zTarget/screenedzTarget);
722 G4double cParameterl3 = 2.*coulombDeflectionl3/(energyLossl3*(energyLossl3+1.));
729 G4double crossSection_L3 = coulombDeflectionFunction_l3 * sigmaPSSR_l3;
735 if (crossSection_L3 >= 0)
737 return crossSection_L3 *
barn;
757 G4cout <<
"*** WARNING in G4ecpssrBaseLixsModel::CalculateVelocity : Proton or Alpha incident particles only. " <<
G4endl;
764 G4double screenedzTarget = zTarget- zlshell;
766 const G4double rydbergMeV= 13.6056923e-6;
770 G4double tetali = (liBindingEnergy*nl*nl)/(screenedzTarget*screenedzTarget*rydbergMeV);
772 G4double reducedEnergy = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
774 G4double velocity = 2.*nl*std::pow(reducedEnergy,0.5)/tetali;
819 std::vector<double>::iterator t2 = std::upper_bound(
dummyVec1.begin(),
dummyVec1.end(), k);
820 std::vector<double>::iterator t1 = t2-1;
822 std::vector<double>::iterator e12 = std::upper_bound(
aVecMap1[(*t1)].begin(),
aVecMap1[(*t1)].end(), theta);
823 std::vector<double>::iterator e11 = e12-1;
825 std::vector<double>::iterator e22 = std::upper_bound(
aVecMap1[(*t2)].begin(),
aVecMap1[(*t2)].end(), theta);
826 std::vector<double>::iterator e21 = e22-1;
835 xs11 =
FL1Data[valueT1][valueE11];
836 xs12 =
FL1Data[valueT1][valueE12];
837 xs21 =
FL1Data[valueT2][valueE21];
838 xs22 =
FL1Data[valueT2][valueE22];
854 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
856 if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0)
return (0.);
910 std::vector<double>::iterator t2 = std::upper_bound(
dummyVec2.begin(),
dummyVec2.end(), k);
911 std::vector<double>::iterator t1 = t2-1;
913 std::vector<double>::iterator e12 = std::upper_bound(
aVecMap2[(*t1)].begin(),
aVecMap2[(*t1)].end(), theta);
914 std::vector<double>::iterator e11 = e12-1;
916 std::vector<double>::iterator e22 = std::upper_bound(
aVecMap2[(*t2)].begin(),
aVecMap2[(*t2)].end(), theta);
917 std::vector<double>::iterator e21 = e22-1;
926 xs11 =
FL2Data[valueT1][valueE11];
927 xs12 =
FL2Data[valueT1][valueE12];
928 xs21 =
FL2Data[valueT2][valueE21];
929 xs22 =
FL2Data[valueT2][valueE22];
945 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
947 if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0)
return (0.);
970 G4double value = xs1 + (xs2 - xs1)*(e - e1)/ (e2 -
e1);
984 G4double value = std::exp(d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
996 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
997 G4double b = std::log10(xs2) - a*std::log10(e2);
998 G4double sigma = a*std::log10(e) + b;
999 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 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