63 path = getenv(
"G4LEDATA");
66 G4Exception(
"G4ecpssrBaseKxsModel::G4ecpssrBaseKxsModel()",
"em0006",
FatalException,
"G4LEDATA environment variable not set" );
70 std::ostringstream fileName;
71 fileName << path <<
"/pixe/uf/FK.dat";
72 std::ifstream FK(fileName.str().c_str());
142 const G4int maxit= 100;
146 if (n<0 || x<0.0 || (x==0.0 && (n==0 || n==1))) {
147 G4cout <<
"*** WARNING in G4ecpssrBaseKxsModel::ExpIntFunction: bad arguments in ExpIntFunction" <<
G4endl;
151 if (n==0) ans=std::exp(-x)/x;
153 if (x==0.0) ans=1.0/nm1;
160 for (i=1;i<=maxit;i++) {
167 if (std::fabs(del-1.0) <
eps) {
173 ans = (nm1!=0 ? 1.0/nm1 : -std::log(x)-euler);
175 for (i=1;i<=maxit;i++) {
177 if (i !=nm1) del = -fact/(i-nm1);
180 for (ii=1;ii<=nm1;ii++) psi +=1.0/ii;
181 del=fact*(-std::log(x)+psi);
184 if (std::fabs(del) < std::fabs(ans)*
eps)
return ans;
222 G4cout <<
"*** WARNING in G4ecpssrBaseKxsModel::CalculateCrossSection : we can treat only Proton or Alpha incident particles " <<
G4endl;
237 G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2;
244 G4double screenedzTarget = zTarget-zkshell;
247 const G4double rydbergMeV= 13.6056923e-6;
249 G4double tetaK = kBindingEnergy/((screenedzTarget*screenedzTarget)*rydbergMeV);
254 G4double velocity =(2./(tetaK*screenedzTarget))*std::pow(((energyIncident*electron_mass_c2)/(massIncident*rydbergMeV)),0.5);
261 const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/
barn ;
265 G4double sigma0 = 8.*
pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.);
271 const G4double kAnalyticalApproximation= 1.5;
272 G4double x = kAnalyticalApproximation/velocity;
283 if ((0.< x) && (x <= 0.035))
285 electrIonizationEnergy= 0.75*
pi*(std::log(1./(x*x))-1.);
289 if ( (0.035 < x) && (x <=3.))
291 electrIonizationEnergy =std::exp(-2.*x)/(0.031+(0.213*std::pow(x,0.5))+(0.005*x)-(0.069*std::pow(x,3./2.))+(0.324*x*x));
296 if ( (3.< x) && (x<=11.))
298 electrIonizationEnergy =2.*std::exp(-2.*x)/std::pow(x,1.6);
301 else electrIonizationEnergy =0.;
307 G4double hFunction =(electrIonizationEnergy*2.)/(tetaK*std::pow(velocity,3));
312 G4double gFunction = (1.+(9.*velocity)+(31.*velocity*velocity)+(98.*std::pow(velocity,3.))+(12.*std::pow(velocity,4.))+(25.*std::pow(velocity,5.))
313 +(4.2*std::pow(velocity,6.))+(0.515*std::pow(velocity,7.)))/std::pow(1.+velocity,9.);
320 G4double sigmaPSS = 1.+(((2.*zIncident)/(screenedzTarget*tetaK))*(gFunction-hFunction));
331 const G4double cNaturalUnit= 1/fine_structure_const;
335 G4double ykFormula=0.4*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(velocity/sigmaPSS);
342 G4double relativityCorrection = std::pow((1.+(1.1*ykFormula*ykFormula)),0.5)+ykFormula;
348 G4double reducedVelocity = velocity*std::pow(relativityCorrection,0.5);
355 G4double etaOverTheta2 = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget)
356 /(sigmaPSS*tetaK)/(sigmaPSS*tetaK);
375 universalFunction = (std::pow(2.,9.)/45.)*std::pow(reducedVelocity/sigmaPSS,8.)*std::pow((1.+(1.72*(reducedVelocity/sigmaPSS)*(reducedVelocity/sigmaPSS))),-4.);
386 if ( etaOverTheta2 > 86.6 && (sigmaPSS*tetaK) > 0.4 && (sigmaPSS*tetaK) < 2.9996 )
402 G4double etaK = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
407 G4double etaT = (sigmaPSS*tetaK)*(sigmaPSS*tetaK)*(86.6);
418 "*** WARNING in G4ecpssrBaseKxsModel::CalculateCrossSection : unable to interpolate FK function in high velocity region ! ***" <<
G4endl;
426 G4double GK = C2/(4*etaK) + C3/(32*etaK*etaK);
430 G4double GT = C2/(4*etaT) + C3/(32*etaT*etaT);
438 G4double fKK = C1*std::log(etaK) + DT - GK;
442 G4double universalFunction3= fKK/(etaK/tetaK);
447 universalFunction=universalFunction3;
451 else if ( etaOverTheta2 >= 1.e-3 && etaOverTheta2 <= 86.6 && (sigmaPSS*tetaK) >= 0.4 && (sigmaPSS*tetaK) <= 2.9996 )
460 if (universalFunction2<=0)
463 "*** WARNING : G4ecpssrBaseKxsModel::CalculateCrossSection is unable to interpolate FK function in medium velocity region ! ***" <<
G4endl;
467 if (
verboseLevel>0)
G4cout <<
" universalFunction2=" << universalFunction2 <<
" for theta=" << sigmaPSS*tetaK <<
" and etaOverTheta2=" << etaOverTheta2 <<
G4endl;
469 universalFunction=universalFunction2;
476 G4double sigmaPSSR = (sigma0/(sigmaPSS*tetaK))*universalFunction;
483 G4double pssDeltaK = (4./(systemMass*sigmaPSS*tetaK))*(sigmaPSS/velocity)*(sigmaPSS/velocity);
489 if (pssDeltaK>1)
return 0.;
491 G4double energyLoss = std::pow(1-pssDeltaK,0.5);
497 G4double energyLossFunction = (std::pow(2.,-9)/8.)*((((9.*energyLoss)-1.)*std::pow(1.+energyLoss,9.))+(((9.*energyLoss)+1.)*std::pow(1.-energyLoss,9.)));
505 G4double coulombDeflection = (4.*
pi*zIncident/systemMass)*std::pow(tetaK*sigmaPSS,-2.)*std::pow(velocity/sigmaPSS,-3.)*(zTarget/screenedzTarget);
510 G4double cParameter = 2.*coulombDeflection/(energyLoss*(energyLoss+1.));
526 crossSection = energyLossFunction* coulombDeflectionFunction*sigmaPSSR;
532 if (crossSection >= 0) {
533 return crossSection *
barn;
577 std::vector<double>::iterator t2 = std::upper_bound(
dummyVec.begin(),
dummyVec.end(), k);
578 std::vector<double>::iterator t1 = t2-1;
580 std::vector<double>::iterator e12 = std::upper_bound(
aVecMap[(*t1)].begin(),
aVecMap[(*t1)].end(), theta);
581 std::vector<double>::iterator e11 = e12-1;
583 std::vector<double>::iterator e22 = std::upper_bound(
aVecMap[(*t2)].begin(),
aVecMap[(*t2)].end(), theta);
584 std::vector<double>::iterator e21 = e22-1;
593 xs11 =
FKData[valueT1][valueE11];
594 xs12 =
FKData[valueT1][valueE12];
595 xs21 =
FKData[valueT2][valueE21];
596 xs22 =
FKData[valueT2][valueE22];
626 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
628 if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0)
return (0.);
653 G4double value = std::exp(d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
665 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
666 G4double b = std::log10(xs2) - a*std::log10(e2);
667 G4double sigma = a*std::log10(e) + b;
668 G4double value = (std::pow(10.,sigma));
G4double FunctionFK(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)
G4double LogLogInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
static const G4double eps
G4double BindingEnergy() const
G4double LinLogInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
std::vector< double > dummyVec
static G4NistManager * Instance()
G4CrossSectionDataSet * tableC3
G4CrossSectionDataSet * tableC1
G4GLOB_DLL std::ostream G4cout
virtual G4double FindValue(G4double e, G4int componentId=0) const
static G4Proton * Proton()
G4double CalculateCrossSection(G4int, G4double, G4double)
G4double ExpIntFunction(G4int n, G4double x)
void print(G4double elem)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
G4double GetPDGMass() const
G4double GetAtomicMassAmu(const G4String &symb) const
G4CrossSectionDataSet * tableC2
static const double eplus
G4double GetPDGCharge() const
static G4AtomicTransitionManager * Instance()
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const
virtual G4bool LoadData(const G4String &argFileName)