46 char *path = std::getenv(
"G4LEDATA");
49 G4Exception(
"G4ecpssrLCrossSection::G4ecpssrBaseLixsModel()",
"em0006",
FatalException ,
"G4LEDATA environment variable not set");
52 std::ostringstream fileName1;
53 std::ostringstream fileName2;
55 fileName1 << path <<
"/pixe/uf/FL1.dat";
56 fileName2 << path <<
"/pixe/uf/FL2.dat";
59 std::ifstream FL1(fileName1.str().c_str());
60 if (!FL1)
G4Exception(
"G4ecpssrLCrossSection::G4ecpssrBaseLixsModel()",
"em0003",
FatalException,
"error opening FL1 data file");
85 std::ifstream FL2(fileName2.str().c_str());
86 if (!FL2)
G4Exception(
"G4ecpssrLCrossSection::G4ecpssrBaseLixsModel()",
"em0003",
FatalException,
" error opening FL2 data file");
133 static const G4double euler= 0.5772156649;
134 static const G4int maxit= 100;
135 static const G4double fpmin = 1.0e-30;
138 if (
n<0 || x<0.0 || (x==0.0 && (
n==0 ||
n==1)))
139 G4cout <<
"*** WARNING in G4ecpssrBaseLixsModel::ExpIntFunction: bad arguments in ExpIntFunction"
142 if (
n==0) ans=
G4Exp(-x)/x;
144 if (x==0.0) ans=1.0/nm1;
151 for (i=1;i<=maxit;i++) {
158 if (std::fabs(del-1.0) <
eps) {
164 ans = (nm1!=0 ? 1.0/nm1 : -std::log(x)-euler);
166 for (i=1;i<=maxit;i++) {
168 if (i !=nm1) del = -fact/(i-nm1);
171 for (ii=1;ii<=nm1;ii++) psi +=1.0/ii;
172 del=fact*(-std::log(x)+psi);
175 if (std::fabs(del) < std::fabs(ans)*
eps)
return ans;
189 if (zTarget <=4)
return 0.;
214 G4cout <<
"*** WARNING in G4ecpssrBaseLixsModel::CalculateL1CrossSection : Proton or Alpha incident particles only. " <<
G4endl;
223 static const G4double zlshell= 4.15;
225 G4double screenedzTarget = zTarget-zlshell;
226 static const G4double rydbergMeV= 13.6056923e-6;
230 G4double tetal1 = (l1BindingEnergy*nl*nl)/((screenedzTarget*screenedzTarget)*rydbergMeV);
235 G4double reducedEnergy = (energyIncident*
electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
241 G4double sigma0 = 8.*
pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.);
249 static const G4double l1AnalyticalApproximation= 1.5;
250 G4double x1 =(nl*l1AnalyticalApproximation)/velocityl1;
256 G4double electrIonizationEnergyl1=0.;
261 if ( x1<=0.035) electrIonizationEnergyl1= 0.75*
pi*(std::log(1./(x1*x1))-1.);
265 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));
267 {
if ( x1<=11.) electrIonizationEnergyl1 =2.*
G4Exp(-2.*x1)/std::pow(x1,1.6);}
270 G4double hFunctionl1 =(electrIonizationEnergyl1*2.*nl)/(tetal1*std::pow(velocityl1,3));
275 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.);
280 G4double sigmaPSS_l1 = 1.+(((2.*zIncident)/(screenedzTarget*tetal1))*(gFunctionl1-hFunctionl1));
289 G4double yl1Formula=0.4*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(nl*velocityl1/sigmaPSS_l1);
294 G4double l1relativityCorrection = std::pow((1.+(1.1*yl1Formula*yl1Formula)),0.5)+yl1Formula;
308 if ( velocityl1 <20. )
311 L1etaOverTheta2 =(reducedEnergy* l1relativityCorrection)/((tetal1*sigmaPSS_l1)*(tetal1*sigmaPSS_l1));
317 if ( ((tetal1*sigmaPSS_l1) >=0.2) && ((tetal1*sigmaPSS_l1) <=2.6670) && (L1etaOverTheta2>=0.1e-3) && (L1etaOverTheta2<=0.866e2) )
318 universalFunction_l1 =
FunctionFL1((tetal1*sigmaPSS_l1),L1etaOverTheta2);
322 sigmaPSSR_l1 = (sigma0/(tetal1*sigmaPSS_l1))*universalFunction_l1;
329 L1etaOverTheta2 = reducedEnergy/(tetal1*tetal1);
335 if ( (tetal1 >=0.2) && (tetal1 <=2.6670) && (L1etaOverTheta2>=0.1e-3) && (L1etaOverTheta2<=0.866e2) )
336 universalFunction_l1 =
FunctionFL1(tetal1,L1etaOverTheta2);
338 if (
verboseLevel>0)
G4cout <<
"at medium and high velocity range, universalFunction_l1 =" << universalFunction_l1 <<
G4endl;
340 sigmaPSSR_l1 = (sigma0/tetal1)*universalFunction_l1;
346 G4double pssDeltal1 = (4./(systemMass *sigmaPSS_l1*tetal1))*(sigmaPSS_l1/velocityl1)*(sigmaPSS_l1/velocityl1);
352 if (pssDeltal1>1)
return 0.;
354 G4double energyLossl1 = std::pow(1-pssDeltal1,0.5);
361 (8.*
pi*zIncident/systemMass)*std::pow(tetal1*sigmaPSS_l1,-2.)*std::pow(velocityl1/sigmaPSS_l1,-3.)*(zTarget/screenedzTarget);
365 G4double cParameterl1 =2.* coulombDeflectionl1/(energyLossl1*(energyLossl1+1.));
372 G4double crossSection_L1 = coulombDeflectionFunction_l1 * sigmaPSSR_l1;
379 if (crossSection_L1 >= 0)
381 return crossSection_L1 *
barn;
392 if (zTarget <=13 )
return 0.;
416 G4cout <<
"*** WARNING in G4ecpssrBaseLixsModel::CalculateL2CrossSection : Proton or Alpha incident particles only. " <<
G4endl;
430 G4double screenedzTarget = zTarget-zlshell;
432 const G4double rydbergMeV= 13.6056923e-6;
436 G4double tetal2 = (l2BindingEnergy*nl*nl)/((screenedzTarget*screenedzTarget)*rydbergMeV);
440 G4double reducedEnergy = (energyIncident*
electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
444 G4double sigma0 = 8.*
pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.);
450 const G4double l23AnalyticalApproximation= 1.25;
452 G4double x2 = (nl*l23AnalyticalApproximation)/velocityl2;
456 G4double electrIonizationEnergyl2=0.;
458 if ( x2<=0.035) electrIonizationEnergyl2= 0.75*
pi*(std::log(1./(x2*x2))-1.);
462 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));
464 {
if ( x2<=11.) electrIonizationEnergyl2 =2.*
G4Exp(-2.*x2)/std::pow(x2,1.6); }
467 G4double hFunctionl2 =(electrIonizationEnergyl2*2.*nl)/(tetal2*std::pow(velocityl2,3));
471 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.);
476 G4double sigmaPSS_l2 = 1.+(((2.*zIncident)/(screenedzTarget*tetal2))*(gFunctionl2-hFunctionl2));
482 G4double yl2Formula=0.15*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(velocityl2/sigmaPSS_l2);
484 G4double l2relativityCorrection = std::pow((1.+(1.1*yl2Formula*yl2Formula)),0.5)+yl2Formula;
492 if ( velocityl2 < 20. )
495 L2etaOverTheta2 = (reducedEnergy*l2relativityCorrection)/((sigmaPSS_l2*tetal2)*(sigmaPSS_l2*tetal2));
497 if ( (tetal2*sigmaPSS_l2>=0.2) && (tetal2*sigmaPSS_l2<=2.6670) && (L2etaOverTheta2>=0.1e-3) && (L2etaOverTheta2<=0.866e2) )
498 universalFunction_l2 =
FunctionFL2((tetal2*sigmaPSS_l2),L2etaOverTheta2);
500 sigmaPSSR_l2 = (sigma0/(tetal2*sigmaPSS_l2))*universalFunction_l2;
507 L2etaOverTheta2 = reducedEnergy /(tetal2*tetal2);
509 if ( (tetal2>=0.2) && (tetal2<=2.6670) && (L2etaOverTheta2>=0.1e-3) && (L2etaOverTheta2<=0.866e2) )
510 universalFunction_l2 =
FunctionFL2((tetal2),L2etaOverTheta2);
512 sigmaPSSR_l2 = (sigma0/tetal2)*universalFunction_l2;
518 G4double pssDeltal2 = (4./(systemMass*sigmaPSS_l2*tetal2))*(sigmaPSS_l2/velocityl2)*(sigmaPSS_l2/velocityl2);
520 if (pssDeltal2>1)
return 0.;
522 G4double energyLossl2 = std::pow(1-pssDeltal2,0.5);
527 =(8.*
pi*zIncident/systemMass)*std::pow(tetal2*sigmaPSS_l2,-2.)*std::pow(velocityl2/sigmaPSS_l2,-3.)*(zTarget/screenedzTarget);
529 G4double cParameterl2 = 2.*coulombDeflectionl2/(energyLossl2*(energyLossl2+1.));
536 G4double crossSection_L2 = coulombDeflectionFunction_l2 * sigmaPSSR_l2;
542 if (crossSection_L2 >= 0)
544 return crossSection_L2 *
barn;
555 if (zTarget <=13)
return 0.;
581 G4cout <<
"*** WARNING in G4ecpssrBaseLixsModel::CalculateL3CrossSection : Proton or Alpha incident particles only. " <<
G4endl;
595 G4double screenedzTarget = zTarget-zlshell;
597 const G4double rydbergMeV= 13.6056923e-6;
601 G4double tetal3 = (l3BindingEnergy*nl*nl)/((screenedzTarget*screenedzTarget)*rydbergMeV);
605 G4double reducedEnergy = (energyIncident*
electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
609 G4double sigma0 = 8.*
pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.);
615 const G4double l23AnalyticalApproximation= 1.25;
617 G4double x3 = (nl*l23AnalyticalApproximation)/velocityl3;
621 G4double electrIonizationEnergyl3=0.;
623 if ( x3<=0.035) electrIonizationEnergyl3= 0.75*
pi*(std::log(1./(x3*x3))-1.);
626 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));
629 if ( x3<=11.) electrIonizationEnergyl3 =2.*
G4Exp(-2.*x3)/std::pow(x3,1.6);}
632 G4double hFunctionl3 =(electrIonizationEnergyl3*2.*nl)/(tetal3*std::pow(velocityl3,3));
636 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.);
641 G4double sigmaPSS_l3 = 1.+(((2.*zIncident)/(screenedzTarget*tetal3))*(gFunctionl3-hFunctionl3));
647 G4double yl3Formula=0.15*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(velocityl3/sigmaPSS_l3);
649 G4double l3relativityCorrection = std::pow((1.+(1.1*yl3Formula*yl3Formula)),0.5)+yl3Formula;
657 if ( velocityl3 < 20. )
660 L3etaOverTheta2 = (reducedEnergy* l3relativityCorrection)/((sigmaPSS_l3*tetal3)*(sigmaPSS_l3*tetal3));
662 if ( (tetal3*sigmaPSS_l3>=0.2) && (tetal3*sigmaPSS_l3<=2.6670) && (L3etaOverTheta2>=0.1e-3) && (L3etaOverTheta2<=0.866e2) )
664 universalFunction_l3 = 2.*
FunctionFL2((tetal3*sigmaPSS_l3), L3etaOverTheta2 );
666 sigmaPSSR_l3 = (sigma0/(tetal3*sigmaPSS_l3))*universalFunction_l3;
676 L3etaOverTheta2 = reducedEnergy/(tetal3*tetal3);
678 if ( (tetal3>=0.2) && (tetal3<=2.6670) && (L3etaOverTheta2>=0.1e-3) && (L3etaOverTheta2<=0.866e2) )
680 universalFunction_l3 = 2.*
FunctionFL2(tetal3, L3etaOverTheta2 );
682 sigmaPSSR_l3 = (sigma0/tetal3)*universalFunction_l3;
687 G4double pssDeltal3 = (4./(systemMass*sigmaPSS_l3*tetal3))*(sigmaPSS_l3/velocityl3)*(sigmaPSS_l3/velocityl3);
691 if (pssDeltal3>1)
return 0.;
693 G4double energyLossl3 = std::pow(1-pssDeltal3,0.5);
698 (8.*
pi*zIncident/systemMass)*std::pow(tetal3*sigmaPSS_l3,-2.)*std::pow(velocityl3/sigmaPSS_l3,-3.)*(zTarget/screenedzTarget);
700 G4double cParameterl3 = 2.*coulombDeflectionl3/(energyLossl3*(energyLossl3+1.));
707 G4double crossSection_L3 = coulombDeflectionFunction_l3 * sigmaPSSR_l3;
713 if (crossSection_L3 >= 0)
715 return crossSection_L3 *
barn;
735 G4cout <<
"*** WARNING in G4ecpssrBaseLixsModel::CalculateVelocity : Proton or Alpha incident particles only. " <<
G4endl;
742 G4double screenedzTarget = zTarget- zlshell;
744 constexpr G4double rydbergMeV= 13.6056923e-6;
748 G4double tetali = (liBindingEnergy*nl*nl)/(screenedzTarget*screenedzTarget*rydbergMeV);
750 G4double reducedEnergy = (energyIncident*
electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
752 G4double velocity = 2.*nl*std::pow(reducedEnergy,0.5)/tetali;
799 auto e12 = std::upper_bound(
aVecMap1[(*t1)].begin(),
aVecMap1[(*t1)].end(), theta);
802 auto e22 = std::upper_bound(
aVecMap1[(*t2)].begin(),
aVecMap1[(*t2)].end(), theta);
812 xs11 =
FL1Data[valueT1][valueE11];
813 xs12 =
FL1Data[valueT1][valueE12];
814 xs21 =
FL1Data[valueT2][valueE21];
815 xs22 =
FL1Data[valueT2][valueE22];
831 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
833 if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0)
return (0.);
889 auto e12 = std::upper_bound(
aVecMap2[(*t1)].begin(),
aVecMap2[(*t1)].end(), theta);
891 auto e22 = std::upper_bound(
aVecMap2[(*t2)].begin(),
aVecMap2[(*t2)].end(), theta);
901 xs11 =
FL2Data[valueT1][valueE11];
902 xs12 =
FL2Data[valueT1][valueE12];
903 xs21 =
FL2Data[valueT2][valueE21];
904 xs22 =
FL2Data[valueT2][valueE22];
920 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
922 if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0)
return (0.);
971 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(
e2)-std::log10(
e1));
972 G4double b = std::log10(xs2) - a*std::log10(
e2);
973 G4double sigma = a*std::log10(e) + b;
974 G4double value = (std::pow(10.,sigma));
static const G4double e1[44]
static const G4double e2[44]
static const G4double eps
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
static constexpr double barn
static constexpr double eplus
static constexpr double pi
G4GLOB_DLL std::ostream G4cout
G4double BindingEnergy() const
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const
static G4AtomicTransitionManager * Instance()
static G4NistManager * Instance()
G4double GetAtomicMassAmu(const G4String &symb) const
G4double GetPDGMass() const
G4double GetPDGCharge() const
static G4Proton * Proton()
G4double CalculateL1CrossSection(G4int zTarget, G4double massIncident, G4double energyIncident) override
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)
G4double FunctionFL2(G4double k, G4double theta)
G4double LinLogInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
G4double LinLinInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
G4double LogLogInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
G4double ExpIntFunction(G4int n, G4double x)
std::vector< G4double > dummyVec1
std::vector< G4double > dummyVec2
G4double CalculateVelocity(G4int subShell, G4int zTarget, G4double massIncident, G4double energyIncident)
G4double FunctionFL1(G4double k, G4double theta)
G4double CalculateL3CrossSection(G4int zTarget, G4double massIncident, G4double energyIncident) override
G4double CalculateL2CrossSection(G4int zTarget, G4double massIncident, G4double energyIncident) override