64 path = getenv(
"G4LEDATA");
67 G4Exception(
"G4ecpssrBaseKxsModel::G4ecpssrBaseKxsModel()",
"em0006",
FatalException,
"G4LEDATA environment variable not set" );
71 std::ostringstream fileName;
72 fileName << path <<
"/pixe/uf/FK.dat";
73 std::ifstream FK(fileName.str().c_str());
143 const G4int maxit= 100;
147 if (n<0 || x<0.0 || (x==0.0 && (n==0 || n==1))) {
148 G4cout <<
"*** WARNING in G4ecpssrBaseKxsModel::ExpIntFunction: bad arguments in ExpIntFunction" <<
G4endl;
152 if (n==0) ans=
G4Exp(-x)/
x;
154 if (x==0.0) ans=1.0/nm1;
161 for (i=1;i<=maxit;i++) {
168 if (std::fabs(del-1.0) <
eps) {
174 ans = (nm1!=0 ? 1.0/nm1 : -std::log(x)-euler);
176 for (i=1;i<=maxit;i++) {
178 if (i !=nm1) del = -fact/(i-nm1);
181 for (ii=1;ii<=nm1;ii++) psi +=1.0/ii;
182 del=fact*(-std::log(x)+psi);
185 if (std::fabs(del) < std::fabs(ans)*
eps)
return ans;
223 G4cout <<
"*** WARNING in G4ecpssrBaseKxsModel::CalculateCrossSection : we can treat only Proton or Alpha incident particles " <<
G4endl;
245 G4double screenedzTarget = zTarget-zkshell;
248 const G4double rydbergMeV= 13.6056923e-6;
250 G4double tetaK = kBindingEnergy/((screenedzTarget*screenedzTarget)*rydbergMeV);
255 G4double velocity =(2./(tetaK*screenedzTarget))*std::pow(((energyIncident*
electron_mass_c2)/(massIncident*rydbergMeV)),0.5);
266 G4double sigma0 = 8.*
pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.);
272 const G4double kAnalyticalApproximation= 1.5;
273 G4double x = kAnalyticalApproximation/velocity;
284 if ((0.< x) && (x <= 0.035))
286 electrIonizationEnergy= 0.75*
pi*(std::log(1./(x*x))-1.);
290 if ( (0.035 < x) && (x <=3.))
292 electrIonizationEnergy =
G4Exp(-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));
297 if ( (3.< x) && (x<=11.))
299 electrIonizationEnergy =2.*
G4Exp(-2.*x)/std::pow(x,1.6);
302 else electrIonizationEnergy =0.;
308 G4double hFunction =(electrIonizationEnergy*2.)/(tetaK*std::pow(velocity,3));
313 G4double gFunction = (1.+(9.*velocity)+(31.*velocity*velocity)+(98.*std::pow(velocity,3.))+(12.*std::pow(velocity,4.))+(25.*std::pow(velocity,5.))
314 +(4.2*std::pow(velocity,6.))+(0.515*std::pow(velocity,7.)))/std::pow(1.+velocity,9.);
321 G4double sigmaPSS = 1.+(((2.*zIncident)/(screenedzTarget*tetaK))*(gFunction-hFunction));
336 G4double ykFormula=0.4*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(velocity/sigmaPSS);
343 G4double relativityCorrection = std::pow((1.+(1.1*ykFormula*ykFormula)),0.5)+ykFormula;
349 G4double reducedVelocity = velocity*std::pow(relativityCorrection,0.5);
356 G4double etaOverTheta2 = (energyIncident*
electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget)
357 /(sigmaPSS*tetaK)/(sigmaPSS*tetaK);
376 universalFunction = (std::pow(2.,9.)/45.)*std::pow(reducedVelocity/sigmaPSS,8.)*std::pow((1.+(1.72*(reducedVelocity/sigmaPSS)*(reducedVelocity/sigmaPSS))),-4.);
387 if ( etaOverTheta2 > 86.6 && (sigmaPSS*tetaK) > 0.4 && (sigmaPSS*tetaK) < 2.9996 )
408 G4double etaT = (sigmaPSS*tetaK)*(sigmaPSS*tetaK)*(86.6);
419 "*** WARNING in G4ecpssrBaseKxsModel::CalculateCrossSection : unable to interpolate FK function in high velocity region ! ***" <<
G4endl;
427 G4double GK = C2/(4*etaK) + C3/(32*etaK*etaK);
431 G4double GT = C2/(4*etaT) + C3/(32*etaT*etaT);
439 G4double fKK = C1*std::log(etaK) + DT - GK;
443 G4double universalFunction3= fKK/(etaK/tetaK);
448 universalFunction=universalFunction3;
452 else if ( etaOverTheta2 >= 1.
e-3 && etaOverTheta2 <= 86.6 && (sigmaPSS*tetaK) >= 0.4 && (sigmaPSS*tetaK) <= 2.9996 )
461 if (universalFunction2<=0)
464 "*** WARNING : G4ecpssrBaseKxsModel::CalculateCrossSection is unable to interpolate FK function in medium velocity region ! ***" <<
G4endl;
468 if (
verboseLevel>0)
G4cout <<
" universalFunction2=" << universalFunction2 <<
" for theta=" << sigmaPSS*tetaK <<
" and etaOverTheta2=" << etaOverTheta2 <<
G4endl;
470 universalFunction=universalFunction2;
477 G4double sigmaPSSR = (sigma0/(sigmaPSS*tetaK))*universalFunction;
484 G4double pssDeltaK = (4./(systemMass*sigmaPSS*tetaK))*(sigmaPSS/velocity)*(sigmaPSS/velocity);
490 if (pssDeltaK>1)
return 0.;
492 G4double energyLoss = std::pow(1-pssDeltaK,0.5);
498 G4double energyLossFunction = (std::pow(2.,-9)/8.)*((((9.*energyLoss)-1.)*std::pow(1.+energyLoss,9.))+(((9.*energyLoss)+1.)*std::pow(1.-energyLoss,9.)));
506 G4double coulombDeflection = (4.*
pi*zIncident/systemMass)*std::pow(tetaK*sigmaPSS,-2.)*std::pow(velocity/sigmaPSS,-3.)*(zTarget/screenedzTarget);
511 G4double cParameter = 2.*coulombDeflection/(energyLoss*(energyLoss+1.));
527 crossSection = energyLossFunction* coulombDeflectionFunction*sigmaPSSR;
533 if (crossSection >= 0) {
534 return crossSection *
barn;
578 std::vector<double>::iterator
t2 = std::upper_bound(
dummyVec.begin(),
dummyVec.end(), k);
579 std::vector<double>::iterator
t1 = t2-1;
581 std::vector<double>::iterator e12 = std::upper_bound(
aVecMap[(*t1)].begin(),
aVecMap[(*t1)].end(), theta);
582 std::vector<double>::iterator e11 = e12-1;
584 std::vector<double>::iterator e22 = std::upper_bound(
aVecMap[(*t2)].begin(),
aVecMap[(*t2)].end(), theta);
585 std::vector<double>::iterator e21 = e22-1;
594 xs11 =
FKData[valueT1][valueE11];
595 xs12 =
FKData[valueT1][valueE12];
596 xs21 =
FKData[valueT2][valueE21];
597 xs22 =
FKData[valueT2][valueE22];
627 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
629 if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0)
return (0.);
666 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
667 G4double b = std::log10(xs2) - a*std::log10(e2);
668 G4double sigma = a*std::log10(e) + b;
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
void print(G4double elem)
std::vector< ExP01TrackerHit * > a
G4double FunctionFK(G4double k, G4double theta)
G4double GetPDGCharge() const
static G4Proton * Proton()
G4double GetAtomicMassAmu(const G4String &symb) const
G4double GetPDGMass() const
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const
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 constexpr double amu_c2
virtual G4double FindValue(G4double e, G4int componentId=0) const
const XML_Char int const XML_Char * value
G4CrossSectionDataSet * tableC3
static constexpr double electron_mass_c2
G4CrossSectionDataSet * tableC2
virtual G4bool LoadData(const G4String &argFileName)
static constexpr double eV
G4double ExpIntFunction(G4int n, G4double x)
std::vector< double > dummyVec
static constexpr double eplus
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
static constexpr double barn
G4double LinLogInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
G4GLOB_DLL std::ostream G4cout
static constexpr double pi
static const G4double eps
static G4AtomicTransitionManager * Instance()
static constexpr double fine_structure_const
static constexpr double Bohr_radius
G4double CalculateCrossSection(G4int, G4double, G4double)
G4CrossSectionDataSet * tableC1
static G4NistManager * Instance()
G4double BindingEnergy() const