00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include <cmath>
00030 #include <iostream>
00031
00032 #include "G4ecpssrBaseKxsModel.hh"
00033
00034 #include "globals.hh"
00035 #include "G4PhysicalConstants.hh"
00036 #include "G4SystemOfUnits.hh"
00037 #include "G4AtomicTransitionManager.hh"
00038 #include "G4NistManager.hh"
00039 #include "G4Proton.hh"
00040 #include "G4Alpha.hh"
00041 #include "G4SemiLogInterpolation.hh"
00042
00043
00044
00045 G4ecpssrBaseKxsModel::G4ecpssrBaseKxsModel()
00046 {
00047 verboseLevel=0;
00048
00049
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
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
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 }
00104
00105
00106
00107 void print (G4double elem)
00108 {
00109 G4cout << elem << " ";
00110 }
00111
00112
00113 G4ecpssrBaseKxsModel::~G4ecpssrBaseKxsModel()
00114 {
00115
00116 delete tableC1;
00117 delete tableC2;
00118 delete tableC3;
00119
00120 }
00121
00122
00123
00124 G4double G4ecpssrBaseKxsModel::ExpIntFunction(G4int n,G4double x)
00125
00126 {
00127
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 }
00192
00193
00194
00195
00196 G4double G4ecpssrBaseKxsModel::CalculateCrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
00197
00198 {
00199
00200
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;
00238
00239 if (verboseLevel>0) G4cout << " systemMass=" << systemMass<< G4endl;
00240
00241 const G4double zkshell= 0.3;
00242
00243
00244 G4double screenedzTarget = zTarget-zkshell;
00245
00246
00247 const G4double rydbergMeV= 13.6056923e-6;
00248
00249 G4double tetaK = kBindingEnergy/((screenedzTarget*screenedzTarget)*rydbergMeV);
00250
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
00256
00257
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.);
00266
00267
00268
00269 if (verboseLevel>0) G4cout << " sigma0=" << sigma0<< G4endl;
00270
00271 const G4double kAnalyticalApproximation= 1.5;
00272 G4double x = kAnalyticalApproximation/velocity;
00273
00274
00275
00276 if (verboseLevel>0) G4cout << " x=" << x<< G4endl;
00277
00278 G4double electrIonizationEnergy;
00279
00280
00281
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));
00308
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.);
00314
00315
00316 if (verboseLevel>0) G4cout << " gFunction=" << gFunction<< G4endl;
00317
00318
00319
00320 G4double sigmaPSS = 1.+(((2.*zIncident)/(screenedzTarget*tetaK))*(gFunction-hFunction));
00321
00322
00323
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;
00332
00333 if (verboseLevel>0) G4cout << " cNaturalUnit=" << cNaturalUnit<< G4endl;
00334
00335 G4double ykFormula=0.4*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(velocity/sigmaPSS);
00336
00337
00338
00339
00340 if (verboseLevel>0) G4cout << " ykFormula=" << ykFormula<< G4endl;
00341
00342 G4double relativityCorrection = std::pow((1.+(1.1*ykFormula*ykFormula)),0.5)+ykFormula;
00343
00344
00345
00346 if (verboseLevel>0) G4cout << " relativityCorrection=" << relativityCorrection<< G4endl;
00347
00348 G4double reducedVelocity = velocity*std::pow(relativityCorrection,0.5);
00349
00350
00351
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
00358
00359
00360 if (verboseLevel>0) G4cout << " etaOverTheta2=" << etaOverTheta2<< G4endl;
00361
00362 G4double universalFunction = 0;
00363
00364
00365
00366 if ( velocity < 1. )
00367
00368
00369
00370
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.);
00376
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
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
00404
00405 if (verboseLevel>0) G4cout << " etaK=" << etaK << G4endl;
00406
00407 G4double etaT = (sigmaPSS*tetaK)*(sigmaPSS*tetaK)*(86.6);
00408
00409
00410 if (verboseLevel>0) G4cout << " etaT=" << etaT << G4endl;
00411
00412 G4double fKT = FunctionFK((sigmaPSS*tetaK),86.6)*(etaT/(sigmaPSS*tetaK));
00413
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
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
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;
00477
00478
00479 if (verboseLevel>0) G4cout << " sigmaPSSR=" << sigmaPSSR<< G4endl;
00480
00481
00482
00483 G4double pssDeltaK = (4./(systemMass*sigmaPSS*tetaK))*(sigmaPSS/velocity)*(sigmaPSS/velocity);
00484
00485
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);
00492
00493
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.)));
00498
00499
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);
00506
00507
00508 if (verboseLevel>0) G4cout << " cParameter-short=" << coulombDeflection<< G4endl;
00509
00510 G4double cParameter = 2.*coulombDeflection/(energyLoss*(energyLoss+1.));
00511
00512
00513 if (verboseLevel>0) G4cout << " cParameter-full=" << cParameter<< G4endl;
00514
00515 G4double coulombDeflectionFunction = 9.*ExpIntFunction(10,cParameter);
00516
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;
00527
00528
00529
00530
00531
00532 if (crossSection >= 0) {
00533 return crossSection * barn;
00534 }
00535 else {return 0;}
00536
00537 }
00538
00539
00540
00541 G4double G4ecpssrBaseKxsModel::FunctionFK(G4double k, G4double theta)
00542 {
00543
00544 G4double sigma = 0.;
00545 G4double valueT1 = 0;
00546 G4double valueT2 = 0;
00547 G4double valueE21 = 0;
00548 G4double valueE22 = 0;
00549 G4double valueE12 = 0;
00550 G4double valueE11 = 0;
00551 G4double xs11 = 0;
00552 G4double xs12 = 0;
00553 G4double xs21 = 0;
00554 G4double xs22 = 0;
00555
00556
00557
00558
00559 if (
00560 theta==8.66e-3 ||
00561 theta==8.66e-2 ||
00562 theta==8.66e-1 ||
00563 theta==8.66e+0 ||
00564 theta==8.66e+1
00565 ) theta=theta-1e-12;
00566
00567 if (
00568 theta==1.e-3 ||
00569 theta==1.e-2 ||
00570 theta==1.e-1 ||
00571 theta==1.e+00 ||
00572 theta==1.e+01
00573 ) theta=theta+1e-12;
00574
00575
00576
00577 std::vector<double>::iterator t2 = std::upper_bound(dummyVec.begin(),dummyVec.end(), k);
00578 std::vector<double>::iterator t1 = t2-1;
00579
00580 std::vector<double>::iterator e12 = std::upper_bound(aVecMap[(*t1)].begin(),aVecMap[(*t1)].end(), theta);
00581 std::vector<double>::iterator e11 = e12-1;
00582
00583 std::vector<double>::iterator e22 = std::upper_bound(aVecMap[(*t2)].begin(),aVecMap[(*t2)].end(), theta);
00584 std::vector<double>::iterator e21 = e22-1;
00585
00586 valueT1 =*t1;
00587 valueT2 =*t2;
00588 valueE21 =*e21;
00589 valueE22 =*e22;
00590 valueE12 =*e12;
00591 valueE11 =*e11;
00592
00593 xs11 = FKData[valueT1][valueE11];
00594 xs12 = FKData[valueT1][valueE12];
00595 xs21 = FKData[valueT2][valueE21];
00596 xs22 = FKData[valueT2][valueE22];
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
00627
00628 if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) return (0.);
00629
00630 if (xsProduct != 0.)
00631 {
00632 sigma = QuadInterpolator( valueE11, valueE12,
00633 valueE21, valueE22,
00634 xs11, xs12,
00635 xs21, xs22,
00636 valueT1, valueT2,
00637 k, theta );
00638 }
00639
00640 return sigma;
00641 }
00642
00643
00644
00645 G4double G4ecpssrBaseKxsModel::LinLogInterpolate(G4double e1,
00646 G4double e2,
00647 G4double e,
00648 G4double xs1,
00649 G4double xs2)
00650 {
00651 G4double d1 = std::log(xs1);
00652 G4double d2 = std::log(xs2);
00653 G4double value = std::exp(d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
00654 return value;
00655 }
00656
00657
00658
00659 G4double G4ecpssrBaseKxsModel::LogLogInterpolate(G4double e1,
00660 G4double e2,
00661 G4double e,
00662 G4double xs1,
00663 G4double xs2)
00664 {
00665 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
00666 G4double b = std::log10(xs2) - a*std::log10(e2);
00667 G4double sigma = a*std::log10(e) + b;
00668 G4double value = (std::pow(10.,sigma));
00669 return value;
00670 }
00671
00672
00673
00674 G4double G4ecpssrBaseKxsModel::QuadInterpolator(G4double e11, G4double e12,
00675 G4double e21, G4double e22,
00676 G4double xs11, G4double xs12,
00677 G4double xs21, G4double xs22,
00678 G4double t1, G4double t2,
00679 G4double t, G4double e)
00680 {
00681
00682 G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
00683 G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
00684 G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
00685
00686
00687
00688
00689
00690
00691
00692 return value;
00693 }
00694
00695
00696
00697
00698