#include <G4ecpssrBaseKxsModel.hh>
Inheritance diagram for G4ecpssrBaseKxsModel:
Public Member Functions | |
G4ecpssrBaseKxsModel () | |
~G4ecpssrBaseKxsModel () | |
G4double | CalculateCrossSection (G4int, G4double, G4double) |
G4double | ExpIntFunction (G4int n, G4double x) |
Definition at line 38 of file G4ecpssrBaseKxsModel.hh.
G4ecpssrBaseKxsModel::G4ecpssrBaseKxsModel | ( | ) |
Definition at line 45 of file G4ecpssrBaseKxsModel.cc.
References FatalException, G4Exception(), and G4CrossSectionDataSet::LoadData().
00046 { 00047 verboseLevel=0; 00048 00049 // Storing C coefficients for high velocity formula 00050 00051 G4String fileC1("pixe/uf/c1"); 00052 tableC1 = new G4CrossSectionDataSet(new G4SemiLogInterpolation, 1.,1.); 00053 00054 G4String fileC2("pixe/uf/c2"); 00055 tableC2 = new G4CrossSectionDataSet(new G4SemiLogInterpolation, 1.,1.); 00056 00057 G4String fileC3("pixe/uf/c3"); 00058 tableC3 = new G4CrossSectionDataSet(new G4SemiLogInterpolation, 1.,1.); 00059 00060 // Storing FK data needed for medium velocities region 00061 char *path = 0; 00062 00063 path = getenv("G4LEDATA"); 00064 00065 if (!path) { 00066 G4Exception("G4ecpssrBaseKxsModel::G4ecpssrBaseKxsModel()", "em0006", FatalException,"G4LEDATA environment variable not set" ); 00067 return; 00068 } 00069 00070 std::ostringstream fileName; 00071 fileName << path << "/pixe/uf/FK.dat"; 00072 std::ifstream FK(fileName.str().c_str()); 00073 00074 if (!FK) 00075 G4Exception("G4ecpssrBaseKxsModel::G4ecpssrBaseKxsModel()", "em0003", FatalException,"error opening FK data file" ); 00076 00077 dummyVec.push_back(0.); 00078 00079 while(!FK.eof()) 00080 { 00081 double x; 00082 double y; 00083 00084 FK>>x>>y; 00085 00086 // Mandatory vector initialization 00087 if (x != dummyVec.back()) 00088 { 00089 dummyVec.push_back(x); 00090 aVecMap[x].push_back(-1.); 00091 } 00092 00093 FK>>FKData[x][y]; 00094 00095 if (y != aVecMap[x].back()) aVecMap[x].push_back(y); 00096 00097 } 00098 00099 tableC1->LoadData(fileC1); 00100 tableC2->LoadData(fileC2); 00101 tableC3->LoadData(fileC3); 00102 00103 }
G4ecpssrBaseKxsModel::~G4ecpssrBaseKxsModel | ( | ) |
Implements G4VecpssrKModel.
Definition at line 196 of file G4ecpssrBaseKxsModel.cc.
References G4Alpha::Alpha(), G4AtomicShell::BindingEnergy(), C1, C3, ExpIntFunction(), G4CrossSectionDataSet::FindValue(), G4cout, G4endl, G4NistManager::GetAtomicMassAmu(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), GT, G4AtomicTransitionManager::Instance(), G4NistManager::Instance(), G4INCL::Math::pi, G4Proton::Proton(), and G4AtomicTransitionManager::Shell().
00198 { 00199 00200 // this K-CrossSection calculation method is done according to W.Brandt and G.Lapicki, Phys.Rev.A23(1981)// 00201 00202 G4NistManager* massManager = G4NistManager::Instance(); 00203 00204 G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance(); 00205 00206 G4double zIncident = 0; 00207 G4Proton* aProtone = G4Proton::Proton(); 00208 G4Alpha* aAlpha = G4Alpha::Alpha(); 00209 00210 if (massIncident == aProtone->GetPDGMass() ) 00211 { 00212 zIncident = (aProtone->GetPDGCharge())/eplus; 00213 } 00214 else 00215 { 00216 if (massIncident == aAlpha->GetPDGMass()) 00217 { 00218 zIncident = (aAlpha->GetPDGCharge())/eplus; 00219 } 00220 else 00221 { 00222 G4cout << "*** WARNING in G4ecpssrBaseKxsModel::CalculateCrossSection : we can treat only Proton or Alpha incident particles " << G4endl; 00223 return 0; 00224 } 00225 } 00226 00227 if (verboseLevel>0) G4cout << " massIncident=" << massIncident<< G4endl; 00228 00229 G4double kBindingEnergy = transitionManager->Shell(zTarget,0)->BindingEnergy(); 00230 00231 if (verboseLevel>0) G4cout << " kBindingEnergy=" << kBindingEnergy/eV<< G4endl; 00232 00233 G4double massTarget = (massManager->GetAtomicMassAmu(zTarget))*amu_c2; 00234 00235 if (verboseLevel>0) G4cout << " massTarget=" << massTarget<< G4endl; 00236 00237 G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2; //the mass of the system (projectile, target) 00238 00239 if (verboseLevel>0) G4cout << " systemMass=" << systemMass<< G4endl; 00240 00241 const G4double zkshell= 0.3; 00242 // *** see Brandt, Phys Rev A23, p 1727 00243 00244 G4double screenedzTarget = zTarget-zkshell; // screenedzTarget is the screened nuclear charge of the target 00245 // *** see Brandt, Phys Rev A23, p 1727 00246 00247 const G4double rydbergMeV= 13.6056923e-6; 00248 00249 G4double tetaK = kBindingEnergy/((screenedzTarget*screenedzTarget)*rydbergMeV); //tetaK denotes the reduced binding energy of the electron 00250 // *** see Rice, ADANDT 20, p 504, f 2 00251 00252 if (verboseLevel>0) G4cout << " tetaK=" << tetaK<< G4endl; 00253 00254 G4double velocity =(2./(tetaK*screenedzTarget))*std::pow(((energyIncident*electron_mass_c2)/(massIncident*rydbergMeV)),0.5); 00255 // *** also called xiK 00256 // *** see Brandt, Phys Rev A23, p 1727 00257 // *** see Basbas, Phys Rev A17, p 1656, f4 00258 00259 if (verboseLevel>0) G4cout << " velocity=" << velocity<< G4endl; 00260 00261 const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/barn ; 00262 00263 if (verboseLevel>0) G4cout << " bohrPow2Barn=" << bohrPow2Barn<< G4endl; 00264 00265 G4double sigma0 = 8.*pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.); //sigma0 is the initial cross section of K shell at stable state 00266 // *** see Benka, ADANDT 22, p 220, f2, for protons 00267 // *** see Basbas, Phys Rev A7, p 1000 00268 00269 if (verboseLevel>0) G4cout << " sigma0=" << sigma0<< G4endl; 00270 00271 const G4double kAnalyticalApproximation= 1.5; 00272 G4double x = kAnalyticalApproximation/velocity; 00273 // *** see Brandt, Phys Rev A23, p 1727 00274 // *** see Brandt, Phys Rev A20, p 469, f16 in expression of h 00275 00276 if (verboseLevel>0) G4cout << " x=" << x<< G4endl; 00277 00278 G4double electrIonizationEnergy; 00279 // *** see Basbas, Phys Rev A17, p1665, f27 00280 // *** see Brandt, Phys Rev A20, p469 00281 // *** see Liu, Comp Phys Comm 97, p325, f A5 00282 00283 if ((0.< x) && (x <= 0.035)) 00284 { 00285 electrIonizationEnergy= 0.75*pi*(std::log(1./(x*x))-1.); 00286 } 00287 else 00288 { 00289 if ( (0.035 < x) && (x <=3.)) 00290 { 00291 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)); 00292 } 00293 00294 else 00295 { 00296 if ( (3.< x) && (x<=11.)) 00297 { 00298 electrIonizationEnergy =2.*std::exp(-2.*x)/std::pow(x,1.6); 00299 } 00300 00301 else electrIonizationEnergy =0.; 00302 } 00303 } 00304 00305 if (verboseLevel>0) G4cout << " electrIonizationEnergy=" << electrIonizationEnergy<< G4endl; 00306 00307 G4double hFunction =(electrIonizationEnergy*2.)/(tetaK*std::pow(velocity,3)); //hFunction represents the correction for polarization effet 00308 // *** see Brandt, Phys Rev A20, p 469, f16 00309 00310 if (verboseLevel>0) G4cout << " hFunction=" << hFunction<< G4endl; 00311 00312 G4double gFunction = (1.+(9.*velocity)+(31.*velocity*velocity)+(98.*std::pow(velocity,3.))+(12.*std::pow(velocity,4.))+(25.*std::pow(velocity,5.)) 00313 +(4.2*std::pow(velocity,6.))+(0.515*std::pow(velocity,7.)))/std::pow(1.+velocity,9.); //gFunction represents the correction for binding effet 00314 // *** see Brandt, Phys Rev A20, p 469, f19 00315 00316 if (verboseLevel>0) G4cout << " gFunction=" << gFunction<< G4endl; 00317 00318 //----------------------------------------------------------------------------------------------------------------------------- 00319 00320 G4double sigmaPSS = 1.+(((2.*zIncident)/(screenedzTarget*tetaK))*(gFunction-hFunction)); //describes the perturbed stationnairy state of the affected atomic electon 00321 // *** also called dzeta 00322 // *** also called epsilon 00323 // *** see Basbas, Phys Rev A17, p1667, f45 00324 00325 if (verboseLevel>0) G4cout << " sigmaPSS=" << sigmaPSS<< G4endl; 00326 00327 if (verboseLevel>0) G4cout << " sigmaPSS*tetaK=" << sigmaPSS*tetaK<< G4endl; 00328 00329 //---------------------------------------------------------------------------------------------------------------------------- 00330 00331 const G4double cNaturalUnit= 1/fine_structure_const; // it's the speed of light according to Atomic-Unit-System 00332 00333 if (verboseLevel>0) G4cout << " cNaturalUnit=" << cNaturalUnit<< G4endl; 00334 00335 G4double ykFormula=0.4*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(velocity/sigmaPSS); 00336 // *** also called yS 00337 // *** see Brandt, Phys Rev A20, p467, f6 00338 // *** see Brandt, Phys Rev A23, p1728 00339 00340 if (verboseLevel>0) G4cout << " ykFormula=" << ykFormula<< G4endl; 00341 00342 G4double relativityCorrection = std::pow((1.+(1.1*ykFormula*ykFormula)),0.5)+ykFormula;// the relativistic correction parameter 00343 // *** also called mRS 00344 // *** see Brandt, Phys Rev A20, p467, f6 00345 00346 if (verboseLevel>0) G4cout << " relativityCorrection=" << relativityCorrection<< G4endl; 00347 00348 G4double reducedVelocity = velocity*std::pow(relativityCorrection,0.5); // presents the reduced collision velocity parameter 00349 // *** also called xiR 00350 // *** see Brandt, Phys Rev A20, p468, f7 00351 // *** see Brandt, Phys Rev A23, p1728 00352 00353 if (verboseLevel>0) G4cout << " reducedVelocity=" << reducedVelocity<< G4endl; 00354 00355 G4double etaOverTheta2 = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget) 00356 /(sigmaPSS*tetaK)/(sigmaPSS*tetaK); 00357 // *** see Benka, ADANDT 22, p220, f4 for eta 00358 // then we use sigmaPSS*tetaK == epsilon*tetaK 00359 00360 if (verboseLevel>0) G4cout << " etaOverTheta2=" << etaOverTheta2<< G4endl; 00361 00362 G4double universalFunction = 0; 00363 00364 // low velocity formula 00365 // ***************** 00366 if ( velocity < 1. ) 00367 // OR 00368 //if ( reducedVelocity/sigmaPSS < 1.) 00369 // *** see Brandt, Phys Rev A23, p1727 00370 // *** reducedVelocity/sigmaPSS is also called xiR/dzeta 00371 // ***************** 00372 { 00373 if (verboseLevel>0) G4cout << " Notice : FK is computed from low velocity formula" << G4endl; 00374 00375 universalFunction = (std::pow(2.,9.)/45.)*std::pow(reducedVelocity/sigmaPSS,8.)*std::pow((1.+(1.72*(reducedVelocity/sigmaPSS)*(reducedVelocity/sigmaPSS))),-4.);// is the reduced universal cross section 00376 // *** see Brandt, Phys Rev A23, p1728 00377 00378 if (verboseLevel>0) G4cout << " universalFunction by Brandt 1981 =" << universalFunction<< G4endl; 00379 00380 } 00381 00382 else 00383 00384 { 00385 00386 if ( etaOverTheta2 > 86.6 && (sigmaPSS*tetaK) > 0.4 && (sigmaPSS*tetaK) < 2.9996 ) 00387 { 00388 // High and medium energies. Method from Rice ADANDT 20, p506, 1977 on tables from Benka 1978 00389 00390 if (verboseLevel>0) G4cout << " Notice : FK is computed from high velocity formula" << G4endl; 00391 00392 if (verboseLevel>0) G4cout << " sigmaPSS*tetaK=" << sigmaPSS*tetaK << G4endl; 00393 00394 G4double C1= tableC1->FindValue(sigmaPSS*tetaK); 00395 G4double C2= tableC2->FindValue(sigmaPSS*tetaK); 00396 G4double C3= tableC3->FindValue(sigmaPSS*tetaK); 00397 00398 if (verboseLevel>0) G4cout << " C1=" << C1 << G4endl; 00399 if (verboseLevel>0) G4cout << " C2=" << C2 << G4endl; 00400 if (verboseLevel>0) G4cout << " C3=" << C3 << G4endl; 00401 00402 G4double etaK = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget); 00403 // *** see Benka, ADANDT 22, p220, f4 for eta 00404 00405 if (verboseLevel>0) G4cout << " etaK=" << etaK << G4endl; 00406 00407 G4double etaT = (sigmaPSS*tetaK)*(sigmaPSS*tetaK)*(86.6); // at any theta, the largest tabulated etaOverTheta2 is 86.6 00408 // *** see Rice, ADANDT 20, p506 00409 00410 if (verboseLevel>0) G4cout << " etaT=" << etaT << G4endl; 00411 00412 G4double fKT = FunctionFK((sigmaPSS*tetaK),86.6)*(etaT/(sigmaPSS*tetaK)); 00413 // *** see Rice, ADANDT 20, p506 00414 00415 if (FunctionFK((sigmaPSS*tetaK),86.6)<=0.) 00416 { 00417 G4cout << 00418 "*** WARNING in G4ecpssrBaseKxsModel::CalculateCrossSection : unable to interpolate FK function in high velocity region ! ***" << G4endl; 00419 return 0; 00420 } 00421 00422 if (verboseLevel>0) G4cout << " FunctionFK=" << FunctionFK((sigmaPSS*tetaK),86.6) << G4endl; 00423 00424 if (verboseLevel>0) G4cout << " fKT=" << fKT << G4endl; 00425 00426 G4double GK = C2/(4*etaK) + C3/(32*etaK*etaK); 00427 00428 if (verboseLevel>0) G4cout << " GK=" << GK << G4endl; 00429 00430 G4double GT = C2/(4*etaT) + C3/(32*etaT*etaT); 00431 00432 if (verboseLevel>0) G4cout << " GT=" << GT << G4endl; 00433 00434 G4double DT = fKT - C1*std::log(etaT) + GT; 00435 00436 if (verboseLevel>0) G4cout << " DT=" << DT << G4endl; 00437 00438 G4double fKK = C1*std::log(etaK) + DT - GK; 00439 00440 if (verboseLevel>0) G4cout << " fKK=" << fKK << G4endl; 00441 00442 G4double universalFunction3= fKK/(etaK/tetaK); 00443 // *** see Rice, ADANDT 20, p505, f7 00444 00445 if (verboseLevel>0) G4cout << " universalFunction3=" << universalFunction3 << G4endl; 00446 00447 universalFunction=universalFunction3; 00448 00449 } 00450 00451 else if ( etaOverTheta2 >= 1.e-3 && etaOverTheta2 <= 86.6 && (sigmaPSS*tetaK) >= 0.4 && (sigmaPSS*tetaK) <= 2.9996 ) 00452 00453 { 00454 // From Benka 1978 00455 00456 if (verboseLevel>0) G4cout << " Notice : FK is computed from INTERPOLATED data" << G4endl; 00457 00458 G4double universalFunction2 = FunctionFK((sigmaPSS*tetaK),etaOverTheta2); 00459 00460 if (universalFunction2<=0) 00461 { 00462 G4cout << 00463 "*** WARNING : G4ecpssrBaseKxsModel::CalculateCrossSection is unable to interpolate FK function in medium velocity region ! ***" << G4endl; 00464 return 0; 00465 } 00466 00467 if (verboseLevel>0) G4cout << " universalFunction2=" << universalFunction2 << " for theta=" << sigmaPSS*tetaK << " and etaOverTheta2=" << etaOverTheta2 << G4endl; 00468 00469 universalFunction=universalFunction2; 00470 } 00471 00472 } 00473 00474 //---------------------------------------------------------------------------------------------------------------------- 00475 00476 G4double sigmaPSSR = (sigma0/(sigmaPSS*tetaK))*universalFunction; //sigmaPSSR is the straight-line K-shell ionization cross section 00477 // *** see Benka, ADANDT 22, p220, f1 00478 00479 if (verboseLevel>0) G4cout << " sigmaPSSR=" << sigmaPSSR<< G4endl; 00480 00481 //----------------------------------------------------------------------------------------------------------------------- 00482 00483 G4double pssDeltaK = (4./(systemMass*sigmaPSS*tetaK))*(sigmaPSS/velocity)*(sigmaPSS/velocity); 00484 // *** also called dzetaK*deltaK 00485 // *** see Brandt, Phys Rev A23, p1727, f B2 00486 00487 if (verboseLevel>0) G4cout << " pssDeltaK=" << pssDeltaK<< G4endl; 00488 00489 if (pssDeltaK>1) return 0.; 00490 00491 G4double energyLoss = std::pow(1-pssDeltaK,0.5); //energyLoss incorporates the straight-line energy-loss 00492 // *** also called zK 00493 // *** see Brandt, Phys Rev A23, p1727, after f B2 00494 00495 if (verboseLevel>0) G4cout << " energyLoss=" << energyLoss<< G4endl; 00496 00497 G4double energyLossFunction = (std::pow(2.,-9)/8.)*((((9.*energyLoss)-1.)*std::pow(1.+energyLoss,9.))+(((9.*energyLoss)+1.)*std::pow(1.-energyLoss,9.)));//energy loss function 00498 // *** also called fs 00499 // *** see Brandt, Phys Rev A23, p1718, f7 00500 00501 if (verboseLevel>0) G4cout << " energyLossFunction=" << energyLossFunction<< G4endl; 00502 00503 //---------------------------------------------------------------------------------------------------------------------------------------------- 00504 00505 G4double coulombDeflection = (4.*pi*zIncident/systemMass)*std::pow(tetaK*sigmaPSS,-2.)*std::pow(velocity/sigmaPSS,-3.)*(zTarget/screenedzTarget); //incorporates Coulomb deflection parameter 00506 // *** see Brandt, Phys Rev A23, p1727, f B3 00507 00508 if (verboseLevel>0) G4cout << " cParameter-short=" << coulombDeflection<< G4endl; 00509 00510 G4double cParameter = 2.*coulombDeflection/(energyLoss*(energyLoss+1.)); 00511 // *** see Brandt, Phys Rev A23, p1727, f B4 00512 00513 if (verboseLevel>0) G4cout << " cParameter-full=" << cParameter<< G4endl; 00514 00515 G4double coulombDeflectionFunction = 9.*ExpIntFunction(10,cParameter); //this function describes Coulomb-deflection effect 00516 // *** see Brandt, Phys Rev A23, p1727 00517 00518 if (verboseLevel>0) G4cout << " ExpIntFunction(10,cParameter) =" << ExpIntFunction(10,cParameter) << G4endl; 00519 00520 if (verboseLevel>0) G4cout << " coulombDeflectionFunction =" << coulombDeflectionFunction << G4endl; 00521 00522 //-------------------------------------------------------------------------------------------------------------------------------------------------- 00523 00524 G4double crossSection = 0; 00525 00526 crossSection = energyLossFunction* coulombDeflectionFunction*sigmaPSSR; //this ECPSSR cross section is estimated at perturbed-stationnairy-state(PSS) 00527 //and it's reduced by the energy-loss(E),the Coulomb deflection(C), 00528 //and the relativity(R) effects 00529 00530 //-------------------------------------------------------------------------------------------------------------------------------------------------- 00531 00532 if (crossSection >= 0) { 00533 return crossSection * barn; 00534 } 00535 else {return 0;} 00536 00537 }
Definition at line 124 of file G4ecpssrBaseKxsModel.cc.
References G4cout, and G4endl.
Referenced by CalculateCrossSection().
00126 { 00127 // this "ExpIntFunction" function allows fast evaluation of the n order exponential integral function En(x) 00128 00129 G4int i; 00130 G4int ii; 00131 G4int nm1; 00132 G4double a; 00133 G4double b; 00134 G4double c; 00135 G4double d; 00136 G4double del; 00137 G4double fact; 00138 G4double h; 00139 G4double psi; 00140 G4double ans = 0; 00141 const G4double euler= 0.5772156649; 00142 const G4int maxit= 100; 00143 const G4double fpmin = 1.0e-30; 00144 const G4double eps = 1.0e-7; 00145 nm1=n-1; 00146 if (n<0 || x<0.0 || (x==0.0 && (n==0 || n==1))) { 00147 G4cout << "*** WARNING in G4ecpssrBaseKxsModel::ExpIntFunction: bad arguments in ExpIntFunction" << G4endl; 00148 G4cout << n << ", " << x << G4endl; 00149 } 00150 else { 00151 if (n==0) ans=std::exp(-x)/x; 00152 else { 00153 if (x==0.0) ans=1.0/nm1; 00154 else { 00155 if (x > 1.0) { 00156 b=x+n; 00157 c=1.0/fpmin; 00158 d=1.0/b; 00159 h=d; 00160 for (i=1;i<=maxit;i++) { 00161 a=-i*(nm1+i); 00162 b +=2.0; 00163 d=1.0/(a*d+b); 00164 c=b+a/c; 00165 del=c*d; 00166 h *=del; 00167 if (std::fabs(del-1.0) < eps) { 00168 ans=h*std::exp(-x); 00169 return ans; 00170 } 00171 } 00172 } else { 00173 ans = (nm1!=0 ? 1.0/nm1 : -std::log(x)-euler); 00174 fact=1.0; 00175 for (i=1;i<=maxit;i++) { 00176 fact *=-x/i; 00177 if (i !=nm1) del = -fact/(i-nm1); 00178 else { 00179 psi = -euler; 00180 for (ii=1;ii<=nm1;ii++) psi +=1.0/ii; 00181 del=fact*(-std::log(x)+psi); 00182 } 00183 ans += del; 00184 if (std::fabs(del) < std::fabs(ans)*eps) return ans; 00185 } 00186 } 00187 } 00188 } 00189 } 00190 return ans; 00191 }