223 G4cout <<
"*** WARNING in G4ecpssrBaseKxsModel::CalculateCrossSection : we can treat only Proton or Alpha incident particles " <<
G4endl;
228 if (verboseLevel>0)
G4cout <<
" massIncident=" << massIncident<<
G4endl;
232 if (verboseLevel>0)
G4cout <<
" kBindingEnergy=" << kBindingEnergy/
eV<<
G4endl;
236 if (verboseLevel>0)
G4cout <<
" massTarget=" << massTarget<<
G4endl;
240 if (verboseLevel>0)
G4cout <<
" systemMass=" << systemMass<<
G4endl;
245 G4double screenedzTarget = zTarget-zkshell;
248 const G4double rydbergMeV= 13.6056923e-6;
250 G4double tetaK = kBindingEnergy/((screenedzTarget*screenedzTarget)*rydbergMeV);
253 if (verboseLevel>0)
G4cout <<
" tetaK=" << tetaK<<
G4endl;
255 G4double velocity =(2./(tetaK*screenedzTarget))*std::pow(((energyIncident*
electron_mass_c2)/(massIncident*rydbergMeV)),0.5);
260 if (verboseLevel>0)
G4cout <<
" velocity=" << velocity<<
G4endl;
264 if (verboseLevel>0)
G4cout <<
" bohrPow2Barn=" << bohrPow2Barn<<
G4endl;
266 G4double sigma0 = 8.*
pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.);
270 if (verboseLevel>0)
G4cout <<
" sigma0=" << sigma0<<
G4endl;
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.;
306 if (verboseLevel>0)
G4cout <<
" electrIonizationEnergy=" << electrIonizationEnergy<<
G4endl;
308 G4double hFunction =(electrIonizationEnergy*2.)/(tetaK*std::pow(velocity,3));
311 if (verboseLevel>0)
G4cout <<
" hFunction=" << hFunction<<
G4endl;
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.);
317 if (verboseLevel>0)
G4cout <<
" gFunction=" << gFunction<<
G4endl;
321 G4double sigmaPSS = 1.+(((2.*zIncident)/(screenedzTarget*tetaK))*(gFunction-hFunction));
326 if (verboseLevel>0)
G4cout <<
" sigmaPSS=" << sigmaPSS<<
G4endl;
328 if (verboseLevel>0)
G4cout <<
" sigmaPSS*tetaK=" << sigmaPSS*tetaK<<
G4endl;
334 if (verboseLevel>0)
G4cout <<
" cNaturalUnit=" << cNaturalUnit<<
G4endl;
336 G4double ykFormula=0.4*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(velocity/sigmaPSS);
341 if (verboseLevel>0)
G4cout <<
" ykFormula=" << ykFormula<<
G4endl;
343 G4double relativityCorrection = std::pow((1.+(1.1*ykFormula*ykFormula)),0.5)+ykFormula;
347 if (verboseLevel>0)
G4cout <<
" relativityCorrection=" << relativityCorrection<<
G4endl;
349 G4double reducedVelocity = velocity*std::pow(relativityCorrection,0.5);
354 if (verboseLevel>0)
G4cout <<
" reducedVelocity=" << reducedVelocity<<
G4endl;
356 G4double etaOverTheta2 = (energyIncident*
electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget)
357 /(sigmaPSS*tetaK)/(sigmaPSS*tetaK);
361 if (verboseLevel>0)
G4cout <<
" etaOverTheta2=" << etaOverTheta2<<
G4endl;
374 if (verboseLevel>0)
G4cout <<
" Notice : FK is computed from low velocity formula" <<
G4endl;
376 universalFunction = (std::pow(2.,9.)/45.)*std::pow(reducedVelocity/sigmaPSS,8.)*std::pow((1.+(1.72*(reducedVelocity/sigmaPSS)*(reducedVelocity/sigmaPSS))),-4.);
379 if (verboseLevel>0)
G4cout <<
" universalFunction by Brandt 1981 =" << universalFunction<<
G4endl;
387 if ( etaOverTheta2 > 86.6 && (sigmaPSS*tetaK) > 0.4 && (sigmaPSS*tetaK) < 2.9996 )
391 if (verboseLevel>0)
G4cout <<
" Notice : FK is computed from high velocity formula" <<
G4endl;
393 if (verboseLevel>0)
G4cout <<
" sigmaPSS*tetaK=" << sigmaPSS*tetaK <<
G4endl;
406 if (verboseLevel>0)
G4cout <<
" etaK=" << etaK <<
G4endl;
408 G4double etaT = (sigmaPSS*tetaK)*(sigmaPSS*tetaK)*(86.6);
411 if (verboseLevel>0)
G4cout <<
" etaT=" << etaT <<
G4endl;
413 G4double fKT = FunctionFK((sigmaPSS*tetaK),86.6)*(etaT/(sigmaPSS*tetaK));
416 if (FunctionFK((sigmaPSS*tetaK),86.6)<=0.)
419 "*** WARNING in G4ecpssrBaseKxsModel::CalculateCrossSection : unable to interpolate FK function in high velocity region ! ***" <<
G4endl;
423 if (verboseLevel>0)
G4cout <<
" FunctionFK=" << FunctionFK((sigmaPSS*tetaK),86.6) <<
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);
446 if (verboseLevel>0)
G4cout <<
" universalFunction3=" << universalFunction3 <<
G4endl;
448 universalFunction=universalFunction3;
452 else if ( etaOverTheta2 >= 1.e-3 && etaOverTheta2 <= 86.6 && (sigmaPSS*tetaK) >= 0.4 && (sigmaPSS*tetaK) <= 2.9996 )
457 if (verboseLevel>0)
G4cout <<
" Notice : FK is computed from INTERPOLATED data" <<
G4endl;
459 G4double universalFunction2 = FunctionFK((sigmaPSS*tetaK),etaOverTheta2);
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;
480 if (verboseLevel>0)
G4cout <<
" sigmaPSSR=" << sigmaPSSR<<
G4endl;
484 G4double pssDeltaK = (4./(systemMass*sigmaPSS*tetaK))*(sigmaPSS/velocity)*(sigmaPSS/velocity);
488 if (verboseLevel>0)
G4cout <<
" pssDeltaK=" << pssDeltaK<<
G4endl;
490 if (pssDeltaK>1)
return 0.;
492 G4double energyLoss = std::pow(1-pssDeltaK,0.5);
496 if (verboseLevel>0)
G4cout <<
" energyLoss=" << energyLoss<<
G4endl;
498 G4double energyLossFunction = (std::pow(2.,-9)/8.)*((((9.*energyLoss)-1.)*std::pow(1.+energyLoss,9.))+(((9.*energyLoss)+1.)*std::pow(1.-energyLoss,9.)));
502 if (verboseLevel>0)
G4cout <<
" energyLossFunction=" << energyLossFunction<<
G4endl;
506 G4double coulombDeflection = (4.*
pi*zIncident/systemMass)*std::pow(tetaK*sigmaPSS,-2.)*std::pow(velocity/sigmaPSS,-3.)*(zTarget/screenedzTarget);
509 if (verboseLevel>0)
G4cout <<
" cParameter-short=" << coulombDeflection<<
G4endl;
511 G4double cParameter = 2.*coulombDeflection/(energyLoss*(energyLoss+1.));
514 if (verboseLevel>0)
G4cout <<
" cParameter-full=" << cParameter<<
G4endl;
521 if (verboseLevel>0)
G4cout <<
" coulombDeflectionFunction =" << coulombDeflectionFunction <<
G4endl;
527 crossSection = energyLossFunction* coulombDeflectionFunction*sigmaPSSR;
533 if (crossSection >= 0) {
534 return crossSection *
barn;
G4double BindingEnergy() const
static G4NistManager * Instance()
G4GLOB_DLL std::ostream G4cout
static constexpr double eplus
virtual G4double FindValue(G4double e, G4int componentId=0) const
static G4Proton * Proton()
static constexpr double eV
G4double ExpIntFunction(G4int n, G4double x)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4double GetPDGMass() const
G4double GetAtomicMassAmu(const G4String &symb) const
static constexpr double pi
static constexpr double barn
G4double GetPDGCharge() const
static G4AtomicTransitionManager * Instance()
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const