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 #include <cmath>
00028 #include <iostream>
00029
00030 #include "G4ecpssrBaseLixsModel.hh"
00031
00032 #include "globals.hh"
00033 #include "G4PhysicalConstants.hh"
00034 #include "G4SystemOfUnits.hh"
00035 #include "G4AtomicTransitionManager.hh"
00036 #include "G4NistManager.hh"
00037 #include "G4Proton.hh"
00038 #include "G4Alpha.hh"
00039 #include "G4LinLogInterpolation.hh"
00040
00041
00042
00043 G4ecpssrBaseLixsModel::G4ecpssrBaseLixsModel()
00044 {
00045 verboseLevel=0;
00046
00047
00048
00049 char *path = getenv("G4LEDATA");
00050
00051 if (!path) {
00052 G4Exception("G4ecpssrLCrossSection::G4ecpssrBaseLixsModel()","em0006", FatalException ,"G4LEDATA environment variable not set");
00053 return;
00054 }
00055 std::ostringstream fileName1;
00056 std::ostringstream fileName2;
00057
00058 fileName1 << path << "/pixe/uf/FL1.dat";
00059 fileName2 << path << "/pixe/uf/FL2.dat";
00060
00061
00062
00063 std::ifstream FL1(fileName1.str().c_str());
00064 if (!FL1) G4Exception("G4ecpssrLCrossSection::G4ecpssrBaseLixsModel()","em0003",FatalException, "error opening FL1 data file");
00065
00066 dummyVec1.push_back(0.);
00067
00068 while(!FL1.eof())
00069 {
00070 double x1;
00071 double y1;
00072
00073 FL1>>x1>>y1;
00074
00075
00076 if (x1 != dummyVec1.back())
00077 {
00078 dummyVec1.push_back(x1);
00079 aVecMap1[x1].push_back(-1.);
00080 }
00081
00082 FL1>>FL1Data[x1][y1];
00083
00084 if (y1 != aVecMap1[x1].back()) aVecMap1[x1].push_back(y1);
00085 }
00086
00087
00088
00089 std::ifstream FL2(fileName2.str().c_str());
00090 if (!FL2) G4Exception("G4ecpssrLCrossSection::G4ecpssrBaseLixsModel()","em0003", FatalException," error opening FL2 data file");
00091
00092 dummyVec2.push_back(0.);
00093
00094 while(!FL2.eof())
00095 {
00096 double x2;
00097 double y2;
00098
00099 FL2>>x2>>y2;
00100
00101
00102 if (x2 != dummyVec2.back())
00103 {
00104 dummyVec2.push_back(x2);
00105 aVecMap2[x2].push_back(-1.);
00106 }
00107
00108 FL2>>FL2Data[x2][y2];
00109
00110 if (y2 != aVecMap2[x2].back()) aVecMap2[x2].push_back(y2);
00111 }
00112
00113 }
00114
00115
00116
00117 G4ecpssrBaseLixsModel::~G4ecpssrBaseLixsModel()
00118 {}
00119
00120
00121
00122 G4double G4ecpssrBaseLixsModel::ExpIntFunction(G4int n,G4double x)
00123
00124 {
00125
00126
00127 G4int i;
00128 G4int ii;
00129 G4int nm1;
00130 G4double a;
00131 G4double b;
00132 G4double c;
00133 G4double d;
00134 G4double del;
00135 G4double fact;
00136 G4double h;
00137 G4double psi;
00138 G4double ans = 0;
00139 const G4double euler= 0.5772156649;
00140 const G4int maxit= 100;
00141 const G4double fpmin = 1.0e-30;
00142 const G4double eps = 1.0e-7;
00143 nm1=n-1;
00144 if (n<0 || x<0.0 || (x==0.0 && (n==0 || n==1)))
00145 G4cout << "*** WARNING in G4ecpssrBaseLixsModel::ExpIntFunction: bad arguments in ExpIntFunction"
00146 << G4endl;
00147 else {
00148 if (n==0) ans=std::exp(-x)/x;
00149 else {
00150 if (x==0.0) ans=1.0/nm1;
00151 else {
00152 if (x > 1.0) {
00153 b=x+n;
00154 c=1.0/fpmin;
00155 d=1.0/b;
00156 h=d;
00157 for (i=1;i<=maxit;i++) {
00158 a=-i*(nm1+i);
00159 b +=2.0;
00160 d=1.0/(a*d+b);
00161 c=b+a/c;
00162 del=c*d;
00163 h *=del;
00164 if (std::fabs(del-1.0) < eps) {
00165 ans=h*std::exp(-x);
00166 return ans;
00167 }
00168 }
00169 } else {
00170 ans = (nm1!=0 ? 1.0/nm1 : -std::log(x)-euler);
00171 fact=1.0;
00172 for (i=1;i<=maxit;i++) {
00173 fact *=-x/i;
00174 if (i !=nm1) del = -fact/(i-nm1);
00175 else {
00176 psi = -euler;
00177 for (ii=1;ii<=nm1;ii++) psi +=1.0/ii;
00178 del=fact*(-std::log(x)+psi);
00179 }
00180 ans += del;
00181 if (std::fabs(del) < std::fabs(ans)*eps) return ans;
00182 }
00183 }
00184 }
00185 }
00186 }
00187 return ans;
00188 }
00189
00190
00191
00192 G4double G4ecpssrBaseLixsModel::CalculateL1CrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
00193 {
00194
00195 if (zTarget <=4) return 0.;
00196
00197
00198
00199
00200 G4NistManager* massManager = G4NistManager::Instance();
00201
00202 G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
00203
00204 G4double zIncident = 0;
00205 G4Proton* aProtone = G4Proton::Proton();
00206 G4Alpha* aAlpha = G4Alpha::Alpha();
00207
00208 if (massIncident == aProtone->GetPDGMass() )
00209
00210 zIncident = (aProtone->GetPDGCharge())/eplus;
00211
00212 else
00213 {
00214 if (massIncident == aAlpha->GetPDGMass())
00215
00216 zIncident = (aAlpha->GetPDGCharge())/eplus;
00217
00218 else
00219 {
00220 G4cout << "*** WARNING in G4ecpssrBaseLixsModel::CalculateL1CrossSection : Proton or Alpha incident particles only. " << G4endl;
00221 G4cout << massIncident << ", " << aAlpha->GetPDGMass() << " (alpha)" << aProtone->GetPDGMass() << " (proton)" << G4endl;
00222 return 0;
00223 }
00224 }
00225
00226 G4double l1BindingEnergy = transitionManager->Shell(zTarget,1)->BindingEnergy();
00227
00228 G4double massTarget = (massManager->GetAtomicMassAmu(zTarget))*amu_c2;
00229
00230 G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2;
00231
00232 const G4double zlshell= 4.15;
00233
00234
00235 G4double screenedzTarget = zTarget-zlshell;
00236
00237 const G4double rydbergMeV= 13.6056923e-6;
00238
00239 const G4double nl= 2.;
00240
00241
00242 G4double tetal1 = (l1BindingEnergy*nl*nl)/((screenedzTarget*screenedzTarget)*rydbergMeV);
00243
00244
00245 if (verboseLevel>0) G4cout << " tetal1=" << tetal1<< G4endl;
00246
00247 G4double reducedEnergy = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
00248
00249
00250
00251 const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/barn ;
00252
00253 G4double sigma0 = 8.*pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.);
00254
00255
00256
00257 G4double velocityl1 = CalculateVelocity(1, zTarget, massIncident, energyIncident);
00258
00259 if (verboseLevel>0) G4cout << " velocityl1=" << velocityl1<< G4endl;
00260
00261 const G4double l1AnalyticalApproximation= 1.5;
00262 G4double x1 =(nl*l1AnalyticalApproximation)/velocityl1;
00263
00264
00265
00266 if (verboseLevel>0) G4cout << " x1=" << x1<< G4endl;
00267
00268 G4double electrIonizationEnergyl1=0.;
00269
00270
00271
00272
00273 if ( x1<=0.035) electrIonizationEnergyl1= 0.75*pi*(std::log(1./(x1*x1))-1.);
00274 else
00275 {
00276 if ( x1<=3.)
00277 electrIonizationEnergyl1 =std::exp(-2.*x1)/(0.031+(0.213*std::pow(x1,0.5))+(0.005*x1)-(0.069*std::pow(x1,3./2.))+(0.324*x1*x1));
00278 else
00279 {if ( x1<=11.) electrIonizationEnergyl1 =2.*std::exp(-2.*x1)/std::pow(x1,1.6);}
00280 }
00281
00282 G4double hFunctionl1 =(electrIonizationEnergyl1*2.*nl)/(tetal1*std::pow(velocityl1,3));
00283
00284
00285 if (verboseLevel>0) G4cout << " hFunctionl1=" << hFunctionl1<< G4endl;
00286
00287 G4double gFunctionl1 = (1.+(9.*velocityl1)+(31.*velocityl1*velocityl1)+(49.*std::pow(velocityl1,3.))+(162.*std::pow(velocityl1,4.))+(63.*std::pow(velocityl1,5.))+(18.*std::pow(velocityl1,6.))+(1.97*std::pow(velocityl1,7.)))/std::pow(1.+velocityl1,9.);
00288
00289
00290 if (verboseLevel>0) G4cout << " gFunctionl1=" << gFunctionl1<< G4endl;
00291
00292 G4double sigmaPSS_l1 = 1.+(((2.*zIncident)/(screenedzTarget*tetal1))*(gFunctionl1-hFunctionl1));
00293
00294
00295
00296
00297 if (verboseLevel>0) G4cout << "sigmaPSS_l1 =" << sigmaPSS_l1<< G4endl;
00298
00299 const G4double cNaturalUnit= 137.;
00300
00301 G4double yl1Formula=0.4*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(nl*velocityl1/sigmaPSS_l1);
00302
00303
00304
00305
00306 G4double l1relativityCorrection = std::pow((1.+(1.1*yl1Formula*yl1Formula)),0.5)+yl1Formula;
00307
00308
00309
00310
00311
00312 G4double L1etaOverTheta2;
00313
00314 G4double universalFunction_l1 = 0.;
00315
00316 G4double sigmaPSSR_l1;
00317
00318
00319
00320 if ( velocityl1 <20. )
00321 {
00322
00323 L1etaOverTheta2 =(reducedEnergy* l1relativityCorrection)/((tetal1*sigmaPSS_l1)*(tetal1*sigmaPSS_l1));
00324
00325
00326
00327
00328
00329 if ( ((tetal1*sigmaPSS_l1) >=0.2) && ((tetal1*sigmaPSS_l1) <=2.6670) && (L1etaOverTheta2>=0.1e-3) && (L1etaOverTheta2<=0.866e2) )
00330
00331 universalFunction_l1 = FunctionFL1((tetal1*sigmaPSS_l1),L1etaOverTheta2);
00332
00333 if (verboseLevel>0) G4cout << "at low velocity range, universalFunction_l1 =" << universalFunction_l1 << G4endl;
00334
00335 sigmaPSSR_l1 = (sigma0/(tetal1*sigmaPSS_l1))*universalFunction_l1;
00336
00337
00338 if (verboseLevel>0) G4cout << " at low velocity range, sigma PWBA L1 CS = " << sigmaPSSR_l1<< G4endl;
00339
00340 }
00341
00342 else
00343
00344 {
00345
00346 L1etaOverTheta2 = reducedEnergy/(tetal1*tetal1);
00347
00348
00349
00350
00351
00352 if ( (tetal1 >=0.2) && (tetal1 <=2.6670) && (L1etaOverTheta2>=0.1e-3) && (L1etaOverTheta2<=0.866e2) )
00353
00354 universalFunction_l1 = FunctionFL1(tetal1,L1etaOverTheta2);
00355
00356 if (verboseLevel>0) G4cout << "at medium and high velocity range, universalFunction_l1 =" << universalFunction_l1 << G4endl;
00357
00358 sigmaPSSR_l1 = (sigma0/tetal1)*universalFunction_l1;
00359
00360
00361 if (verboseLevel>0) G4cout << " sigma PWBA L1 CS at medium and high velocity range = " << sigmaPSSR_l1<< G4endl;
00362 }
00363
00364 G4double pssDeltal1 = (4./(systemMass *sigmaPSS_l1*tetal1))*(sigmaPSS_l1/velocityl1)*(sigmaPSS_l1/velocityl1);
00365
00366
00367
00368 if (verboseLevel>0) G4cout << " pssDeltal1=" << pssDeltal1<< G4endl;
00369
00370 if (pssDeltal1>1) return 0.;
00371
00372 G4double energyLossl1 = std::pow(1-pssDeltal1,0.5);
00373
00374
00375
00376 if (verboseLevel>0) G4cout << " energyLossl1=" << energyLossl1<< G4endl;
00377
00378 G4double coulombDeflectionl1 =
00379 (8.*pi*zIncident/systemMass)*std::pow(tetal1*sigmaPSS_l1,-2.)*std::pow(velocityl1/sigmaPSS_l1,-3.)*(zTarget/screenedzTarget);
00380
00381
00382
00383 G4double cParameterl1 =2.* coulombDeflectionl1/(energyLossl1*(energyLossl1+1.));
00384
00385
00386 G4double coulombDeflectionFunction_l1 = 9.*ExpIntFunction(10,cParameterl1);
00387
00388 if (verboseLevel>0) G4cout << " coulombDeflectionFunction_l1 =" << coulombDeflectionFunction_l1 << G4endl;
00389
00390 G4double crossSection_L1 = coulombDeflectionFunction_l1 * sigmaPSSR_l1;
00391
00392
00393
00394
00395 if (verboseLevel>0) G4cout << " crossSection_L1 =" << crossSection_L1 << G4endl;
00396
00397 if (crossSection_L1 >= 0)
00398 {
00399 return crossSection_L1 * barn;
00400 }
00401
00402 else {return 0;}
00403 }
00404
00405
00406
00407 G4double G4ecpssrBaseLixsModel::CalculateL2CrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
00408
00409 {
00410 if (zTarget <=13 ) return 0.;
00411
00412
00413
00414
00415 G4NistManager* massManager = G4NistManager::Instance();
00416
00417 G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
00418
00419 G4double zIncident = 0;
00420
00421 G4Proton* aProtone = G4Proton::Proton();
00422 G4Alpha* aAlpha = G4Alpha::Alpha();
00423
00424 if (massIncident == aProtone->GetPDGMass() )
00425
00426 zIncident = (aProtone->GetPDGCharge())/eplus;
00427
00428 else
00429 {
00430 if (massIncident == aAlpha->GetPDGMass())
00431
00432 zIncident = (aAlpha->GetPDGCharge())/eplus;
00433
00434 else
00435 {
00436 G4cout << "*** WARNING in G4ecpssrBaseLixsModel::CalculateL2CrossSection : Proton or Alpha incident particles only. " << G4endl;
00437 G4cout << massIncident << ", " << aAlpha->GetPDGMass() << " (alpha)" << aProtone->GetPDGMass() << " (proton)" << G4endl;
00438 return 0;
00439 }
00440 }
00441
00442 G4double l2BindingEnergy = transitionManager->Shell(zTarget,2)->BindingEnergy();
00443
00444 G4double massTarget = (massManager->GetAtomicMassAmu(zTarget))*amu_c2;
00445
00446 G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2;
00447
00448 const G4double zlshell= 4.15;
00449
00450 G4double screenedzTarget = zTarget-zlshell;
00451
00452 const G4double rydbergMeV= 13.6056923e-6;
00453
00454 const G4double nl= 2.;
00455
00456 G4double tetal2 = (l2BindingEnergy*nl*nl)/((screenedzTarget*screenedzTarget)*rydbergMeV);
00457
00458 if (verboseLevel>0) G4cout << " tetal2=" << tetal2<< G4endl;
00459
00460 G4double reducedEnergy = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
00461
00462 const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/barn ;
00463
00464 G4double sigma0 = 8.*pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.);
00465
00466 G4double velocityl2 = CalculateVelocity(2, zTarget, massIncident, energyIncident);
00467
00468 if (verboseLevel>0) G4cout << " velocityl2=" << velocityl2<< G4endl;
00469
00470 const G4double l23AnalyticalApproximation= 1.25;
00471
00472 G4double x2 = (nl*l23AnalyticalApproximation)/velocityl2;
00473
00474 if (verboseLevel>0) G4cout << " x2=" << x2<< G4endl;
00475
00476 G4double electrIonizationEnergyl2=0.;
00477
00478 if ( x2<=0.035) electrIonizationEnergyl2= 0.75*pi*(std::log(1./(x2*x2))-1.);
00479 else
00480 {
00481 if ( x2<=3.)
00482 electrIonizationEnergyl2 =std::exp(-2.*x2)/(0.031+(0.213*std::pow(x2,0.5))+(0.005*x2)-(0.069*std::pow(x2,3./2.))+(0.324*x2*x2));
00483 else
00484 {if ( x2<=11.) electrIonizationEnergyl2 =2.*std::exp(-2.*x2)/std::pow(x2,1.6); }
00485 }
00486
00487 G4double hFunctionl2 =(electrIonizationEnergyl2*2.*nl)/(tetal2*std::pow(velocityl2,3));
00488
00489 if (verboseLevel>0) G4cout << " hFunctionl2=" << hFunctionl2<< G4endl;
00490
00491 G4double gFunctionl2 = (1.+(10.*velocityl2)+(45.*velocityl2*velocityl2)+(102.*std::pow(velocityl2,3.))+(331.*std::pow(velocityl2,4.))+(6.7*std::pow(velocityl2,5.))+(58.*std::pow(velocityl2,6.))+(7.8*std::pow(velocityl2,7.))+ (0.888*std::pow(velocityl2,8.)) )/std::pow(1.+velocityl2,10.);
00492
00493
00494 if (verboseLevel>0) G4cout << " gFunctionl2=" << gFunctionl2<< G4endl;
00495
00496 G4double sigmaPSS_l2 = 1.+(((2.*zIncident)/(screenedzTarget*tetal2))*(gFunctionl2-hFunctionl2));
00497
00498 if (verboseLevel>0) G4cout << " sigmaPSS_l2=" << sigmaPSS_l2<< G4endl;
00499
00500 const G4double cNaturalUnit= 137.;
00501
00502 G4double yl2Formula=0.15*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(velocityl2/sigmaPSS_l2);
00503
00504 G4double l2relativityCorrection = std::pow((1.+(1.1*yl2Formula*yl2Formula)),0.5)+yl2Formula;
00505
00506 G4double L2etaOverTheta2;
00507
00508 G4double universalFunction_l2 = 0.;
00509
00510 G4double sigmaPSSR_l2 ;
00511
00512 if ( velocityl2 < 20. )
00513 {
00514
00515 L2etaOverTheta2 = (reducedEnergy*l2relativityCorrection)/((sigmaPSS_l2*tetal2)*(sigmaPSS_l2*tetal2));
00516
00517 if ( (tetal2*sigmaPSS_l2>=0.2) && (tetal2*sigmaPSS_l2<=2.6670) && (L2etaOverTheta2>=0.1e-3) && (L2etaOverTheta2<=0.866e2) )
00518
00519 universalFunction_l2 = FunctionFL2((tetal2*sigmaPSS_l2),L2etaOverTheta2);
00520
00521 sigmaPSSR_l2 = (sigma0/(tetal2*sigmaPSS_l2))*universalFunction_l2;
00522
00523 if (verboseLevel>0) G4cout << " sigma PWBA L2 CS at low velocity range = " << sigmaPSSR_l2<< G4endl;
00524 }
00525 else
00526 {
00527
00528 L2etaOverTheta2 = reducedEnergy /(tetal2*tetal2);
00529
00530 if ( (tetal2>=0.2) && (tetal2<=2.6670) && (L2etaOverTheta2>=0.1e-3) && (L2etaOverTheta2<=0.866e2) )
00531
00532 universalFunction_l2 = FunctionFL2((tetal2),L2etaOverTheta2);
00533
00534 sigmaPSSR_l2 = (sigma0/tetal2)*universalFunction_l2;
00535
00536 if (verboseLevel>0) G4cout << " sigma PWBA L2 CS at medium and high velocity range = " << sigmaPSSR_l2<< G4endl;
00537
00538 }
00539
00540 G4double pssDeltal2 = (4./(systemMass*sigmaPSS_l2*tetal2))*(sigmaPSS_l2/velocityl2)*(sigmaPSS_l2/velocityl2);
00541
00542 if (pssDeltal2>1) return 0.;
00543
00544 G4double energyLossl2 = std::pow(1-pssDeltal2,0.5);
00545
00546 if (verboseLevel>0) G4cout << " energyLossl2=" << energyLossl2<< G4endl;
00547
00548 G4double coulombDeflectionl2
00549 =(8.*pi*zIncident/systemMass)*std::pow(tetal2*sigmaPSS_l2,-2.)*std::pow(velocityl2/sigmaPSS_l2,-3.)*(zTarget/screenedzTarget);
00550
00551 G4double cParameterl2 = 2.*coulombDeflectionl2/(energyLossl2*(energyLossl2+1.));
00552
00553 G4double coulombDeflectionFunction_l2 = 11.*ExpIntFunction(12,cParameterl2);
00554
00555
00556 if (verboseLevel>0) G4cout << " coulombDeflectionFunction_l2 =" << coulombDeflectionFunction_l2 << G4endl;
00557
00558 G4double crossSection_L2 = coulombDeflectionFunction_l2 * sigmaPSSR_l2;
00559
00560
00561
00562 if (verboseLevel>0) G4cout << " crossSection_L2 =" << crossSection_L2 << G4endl;
00563
00564 if (crossSection_L2 >= 0)
00565 {
00566 return crossSection_L2 * barn;
00567 }
00568 else {return 0;}
00569 }
00570
00571
00572
00573
00574 G4double G4ecpssrBaseLixsModel::CalculateL3CrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
00575
00576 {
00577 if (zTarget <=13) return 0.;
00578
00579
00580
00581
00582 G4NistManager* massManager = G4NistManager::Instance();
00583
00584 G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
00585
00586 G4double zIncident = 0;
00587
00588 G4Proton* aProtone = G4Proton::Proton();
00589 G4Alpha* aAlpha = G4Alpha::Alpha();
00590
00591 if (massIncident == aProtone->GetPDGMass() )
00592
00593 zIncident = (aProtone->GetPDGCharge())/eplus;
00594
00595 else
00596 {
00597 if (massIncident == aAlpha->GetPDGMass())
00598
00599 zIncident = (aAlpha->GetPDGCharge())/eplus;
00600
00601 else
00602 {
00603 G4cout << "*** WARNING in G4ecpssrBaseLixsModel::CalculateL3CrossSection : Proton or Alpha incident particles only. " << G4endl;
00604 G4cout << massIncident << ", " << aAlpha->GetPDGMass() << " (alpha)" << aProtone->GetPDGMass() << " (proton)" << G4endl;
00605 return 0;
00606 }
00607 }
00608
00609 G4double l3BindingEnergy = transitionManager->Shell(zTarget,3)->BindingEnergy();
00610
00611 G4double massTarget = (massManager->GetAtomicMassAmu(zTarget))*amu_c2;
00612
00613 G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2;
00614
00615 const G4double zlshell= 4.15;
00616
00617 G4double screenedzTarget = zTarget-zlshell;
00618
00619 const G4double rydbergMeV= 13.6056923e-6;
00620
00621 const G4double nl= 2.;
00622
00623 G4double tetal3 = (l3BindingEnergy*nl*nl)/((screenedzTarget*screenedzTarget)*rydbergMeV);
00624
00625 if (verboseLevel>0) G4cout << " tetal3=" << tetal3<< G4endl;
00626
00627 G4double reducedEnergy = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
00628
00629 const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/barn ;
00630
00631 G4double sigma0 = 8.*pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.);
00632
00633 G4double velocityl3 = CalculateVelocity(3, zTarget, massIncident, energyIncident);
00634
00635 if (verboseLevel>0) G4cout << " velocityl3=" << velocityl3<< G4endl;
00636
00637 const G4double l23AnalyticalApproximation= 1.25;
00638
00639 G4double x3 = (nl*l23AnalyticalApproximation)/velocityl3;
00640
00641 if (verboseLevel>0) G4cout << " x3=" << x3<< G4endl;
00642
00643 G4double electrIonizationEnergyl3=0.;
00644
00645 if ( x3<=0.035) electrIonizationEnergyl3= 0.75*pi*(std::log(1./(x3*x3))-1.);
00646 else
00647 {
00648 if ( x3<=3.) electrIonizationEnergyl3 =std::exp(-2.*x3)/(0.031+(0.213*std::pow(x3,0.5))+(0.005*x3)-(0.069*std::pow(x3,3./2.))+(0.324*x3*x3));
00649 else
00650 {
00651 if ( x3<=11.) electrIonizationEnergyl3 =2.*std::exp(-2.*x3)/std::pow(x3,1.6);}
00652 }
00653
00654 G4double hFunctionl3 =(electrIonizationEnergyl3*2.*nl)/(tetal3*std::pow(velocityl3,3));
00655
00656 if (verboseLevel>0) G4cout << " hFunctionl3=" << hFunctionl3<< G4endl;
00657
00658 G4double gFunctionl3 = (1.+(10.*velocityl3)+(45.*velocityl3*velocityl3)+(102.*std::pow(velocityl3,3.))+(331.*std::pow(velocityl3,4.))+(6.7*std::pow(velocityl3,5.))+(58.*std::pow(velocityl3,6.))+(7.8*std::pow(velocityl3,7.))+ (0.888*std::pow(velocityl3,8.)) )/std::pow(1.+velocityl3,10.);
00659
00660
00661 if (verboseLevel>0) G4cout << " gFunctionl3=" << gFunctionl3<< G4endl;
00662
00663 G4double sigmaPSS_l3 = 1.+(((2.*zIncident)/(screenedzTarget*tetal3))*(gFunctionl3-hFunctionl3));
00664
00665 if (verboseLevel>0) G4cout << "sigmaPSS_l3 =" << sigmaPSS_l3<< G4endl;
00666
00667 const G4double cNaturalUnit= 137.;
00668
00669 G4double yl3Formula=0.15*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(velocityl3/sigmaPSS_l3);
00670
00671 G4double l3relativityCorrection = std::pow((1.+(1.1*yl3Formula*yl3Formula)),0.5)+yl3Formula;
00672
00673 G4double L3etaOverTheta2;
00674
00675 G4double universalFunction_l3 = 0.;
00676
00677 G4double sigmaPSSR_l3;
00678
00679 if ( velocityl3 < 20. )
00680 {
00681
00682 L3etaOverTheta2 = (reducedEnergy* l3relativityCorrection)/((sigmaPSS_l3*tetal3)*(sigmaPSS_l3*tetal3));
00683
00684 if ( (tetal3*sigmaPSS_l3>=0.2) && (tetal3*sigmaPSS_l3<=2.6670) && (L3etaOverTheta2>=0.1e-3) && (L3etaOverTheta2<=0.866e2) )
00685
00686 universalFunction_l3 = 2.*FunctionFL2((tetal3*sigmaPSS_l3), L3etaOverTheta2 );
00687
00688 sigmaPSSR_l3 = (sigma0/(tetal3*sigmaPSS_l3))*universalFunction_l3;
00689
00690 if (verboseLevel>0) G4cout << " sigma PWBA L3 CS at low velocity range = " << sigmaPSSR_l3<< G4endl;
00691
00692 }
00693
00694 else
00695
00696 {
00697
00698 L3etaOverTheta2 = reducedEnergy/(tetal3*tetal3);
00699
00700 if ( (tetal3>=0.2) && (tetal3<=2.6670) && (L3etaOverTheta2>=0.1e-3) && (L3etaOverTheta2<=0.866e2) )
00701
00702 universalFunction_l3 = 2.*FunctionFL2(tetal3, L3etaOverTheta2 );
00703
00704 sigmaPSSR_l3 = (sigma0/tetal3)*universalFunction_l3;
00705
00706 if (verboseLevel>0) G4cout << " sigma PWBA L3 CS at medium and high velocity range = " << sigmaPSSR_l3<< G4endl;
00707 }
00708
00709 G4double pssDeltal3 = (4./(systemMass*sigmaPSS_l3*tetal3))*(sigmaPSS_l3/velocityl3)*(sigmaPSS_l3/velocityl3);
00710
00711 if (verboseLevel>0) G4cout << " pssDeltal3=" << pssDeltal3<< G4endl;
00712
00713 if (pssDeltal3>1) return 0.;
00714
00715 G4double energyLossl3 = std::pow(1-pssDeltal3,0.5);
00716
00717 if (verboseLevel>0) G4cout << " energyLossl3=" << energyLossl3<< G4endl;
00718
00719 G4double coulombDeflectionl3 =
00720 (8.*pi*zIncident/systemMass)*std::pow(tetal3*sigmaPSS_l3,-2.)*std::pow(velocityl3/sigmaPSS_l3,-3.)*(zTarget/screenedzTarget);
00721
00722 G4double cParameterl3 = 2.*coulombDeflectionl3/(energyLossl3*(energyLossl3+1.));
00723
00724 G4double coulombDeflectionFunction_l3 = 11.*ExpIntFunction(12,cParameterl3);
00725
00726
00727 if (verboseLevel>0) G4cout << " coulombDeflectionFunction_l3 =" << coulombDeflectionFunction_l3 << G4endl;
00728
00729 G4double crossSection_L3 = coulombDeflectionFunction_l3 * sigmaPSSR_l3;
00730
00731
00732
00733 if (verboseLevel>0) G4cout << " crossSection_L3 =" << crossSection_L3 << G4endl;
00734
00735 if (crossSection_L3 >= 0)
00736 {
00737 return crossSection_L3 * barn;
00738 }
00739 else {return 0;}
00740 }
00741
00742
00743
00744 G4double G4ecpssrBaseLixsModel::CalculateVelocity(G4int subShell, G4int zTarget, G4double massIncident, G4double energyIncident)
00745
00746 {
00747
00748 G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
00749
00750 G4double liBindingEnergy = transitionManager->Shell(zTarget,subShell)->BindingEnergy();
00751
00752 G4Proton* aProtone = G4Proton::Proton();
00753 G4Alpha* aAlpha = G4Alpha::Alpha();
00754
00755 if (!((massIncident == aProtone->GetPDGMass()) || (massIncident == aAlpha->GetPDGMass())))
00756 {
00757 G4cout << "*** WARNING in G4ecpssrBaseLixsModel::CalculateVelocity : Proton or Alpha incident particles only. " << G4endl;
00758 G4cout << massIncident << ", " << aAlpha->GetPDGMass() << " (alpha)" << aProtone->GetPDGMass() << " (proton)" << G4endl;
00759 return 0;
00760 }
00761
00762 const G4double zlshell= 4.15;
00763
00764 G4double screenedzTarget = zTarget- zlshell;
00765
00766 const G4double rydbergMeV= 13.6056923e-6;
00767
00768 const G4double nl= 2.;
00769
00770 G4double tetali = (liBindingEnergy*nl*nl)/(screenedzTarget*screenedzTarget*rydbergMeV);
00771
00772 G4double reducedEnergy = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
00773
00774 G4double velocity = 2.*nl*std::pow(reducedEnergy,0.5)/tetali;
00775
00776
00777 return velocity;
00778 }
00779
00780
00781
00782 G4double G4ecpssrBaseLixsModel::FunctionFL1(G4double k, G4double theta)
00783 {
00784
00785 G4double sigma = 0.;
00786 G4double valueT1 = 0;
00787 G4double valueT2 = 0;
00788 G4double valueE21 = 0;
00789 G4double valueE22 = 0;
00790 G4double valueE12 = 0;
00791 G4double valueE11 = 0;
00792 G4double xs11 = 0;
00793 G4double xs12 = 0;
00794 G4double xs21 = 0;
00795 G4double xs22 = 0;
00796
00797
00798
00799 if (
00800 theta==8.66e-4 ||
00801 theta==8.66e-3 ||
00802 theta==8.66e-2 ||
00803 theta==8.66e-1 ||
00804 theta==8.66e+00 ||
00805 theta==8.66e+01
00806 ) theta=theta-1e-12;
00807
00808 if (
00809 theta==1.e-4 ||
00810 theta==1.e-3 ||
00811 theta==1.e-2 ||
00812 theta==1.e-1 ||
00813 theta==1.e+00 ||
00814 theta==1.e+01
00815 ) theta=theta+1e-12;
00816
00817
00818
00819 std::vector<double>::iterator t2 = std::upper_bound(dummyVec1.begin(),dummyVec1.end(), k);
00820 std::vector<double>::iterator t1 = t2-1;
00821
00822 std::vector<double>::iterator e12 = std::upper_bound(aVecMap1[(*t1)].begin(),aVecMap1[(*t1)].end(), theta);
00823 std::vector<double>::iterator e11 = e12-1;
00824
00825 std::vector<double>::iterator e22 = std::upper_bound(aVecMap1[(*t2)].begin(),aVecMap1[(*t2)].end(), theta);
00826 std::vector<double>::iterator e21 = e22-1;
00827
00828 valueT1 =*t1;
00829 valueT2 =*t2;
00830 valueE21 =*e21;
00831 valueE22 =*e22;
00832 valueE12 =*e12;
00833 valueE11 =*e11;
00834
00835 xs11 = FL1Data[valueT1][valueE11];
00836 xs12 = FL1Data[valueT1][valueE12];
00837 xs21 = FL1Data[valueT2][valueE21];
00838 xs22 = FL1Data[valueT2][valueE22];
00839
00840 if (verboseLevel>0)
00841 G4cout
00842 << valueT1 << " "
00843 << valueT2 << " "
00844 << valueE11 << " "
00845 << valueE12 << " "
00846 << valueE21 << " "
00847 << valueE22 << " "
00848 << xs11 << " "
00849 << xs12 << " "
00850 << xs21 << " "
00851 << xs22 << " "
00852 << G4endl;
00853
00854 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
00855
00856 if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) return (0.);
00857
00858 if (xsProduct != 0.)
00859 {
00860 sigma = QuadInterpolator( valueE11, valueE12,
00861 valueE21, valueE22,
00862 xs11, xs12,
00863 xs21, xs22,
00864 valueT1, valueT2,
00865 k, theta );
00866 }
00867
00868 return sigma;
00869 }
00870
00871
00872
00873 G4double G4ecpssrBaseLixsModel::FunctionFL2(G4double k, G4double theta)
00874 {
00875
00876 G4double sigma = 0.;
00877 G4double valueT1 = 0;
00878 G4double valueT2 = 0;
00879 G4double valueE21 = 0;
00880 G4double valueE22 = 0;
00881 G4double valueE12 = 0;
00882 G4double valueE11 = 0;
00883 G4double xs11 = 0;
00884 G4double xs12 = 0;
00885 G4double xs21 = 0;
00886 G4double xs22 = 0;
00887
00888
00889
00890 if (
00891 theta==8.66e-4 ||
00892 theta==8.66e-3 ||
00893 theta==8.66e-2 ||
00894 theta==8.66e-1 ||
00895 theta==8.66e+00 ||
00896 theta==8.66e+01
00897 ) theta=theta-1e-12;
00898
00899 if (
00900 theta==1.e-4 ||
00901 theta==1.e-3 ||
00902 theta==1.e-2 ||
00903 theta==1.e-1 ||
00904 theta==1.e+00 ||
00905 theta==1.e+01
00906 ) theta=theta+1e-12;
00907
00908
00909
00910 std::vector<double>::iterator t2 = std::upper_bound(dummyVec2.begin(),dummyVec2.end(), k);
00911 std::vector<double>::iterator t1 = t2-1;
00912
00913 std::vector<double>::iterator e12 = std::upper_bound(aVecMap2[(*t1)].begin(),aVecMap2[(*t1)].end(), theta);
00914 std::vector<double>::iterator e11 = e12-1;
00915
00916 std::vector<double>::iterator e22 = std::upper_bound(aVecMap2[(*t2)].begin(),aVecMap2[(*t2)].end(), theta);
00917 std::vector<double>::iterator e21 = e22-1;
00918
00919 valueT1 =*t1;
00920 valueT2 =*t2;
00921 valueE21 =*e21;
00922 valueE22 =*e22;
00923 valueE12 =*e12;
00924 valueE11 =*e11;
00925
00926 xs11 = FL2Data[valueT1][valueE11];
00927 xs12 = FL2Data[valueT1][valueE12];
00928 xs21 = FL2Data[valueT2][valueE21];
00929 xs22 = FL2Data[valueT2][valueE22];
00930
00931 if (verboseLevel>0)
00932 G4cout
00933 << valueT1 << " "
00934 << valueT2 << " "
00935 << valueE11 << " "
00936 << valueE12 << " "
00937 << valueE21 << " "
00938 << valueE22 << " "
00939 << xs11 << " "
00940 << xs12 << " "
00941 << xs21 << " "
00942 << xs22 << " "
00943 << G4endl;
00944
00945 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
00946
00947 if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) return (0.);
00948
00949 if (xsProduct != 0.)
00950 {
00951 sigma = QuadInterpolator( valueE11, valueE12,
00952 valueE21, valueE22,
00953 xs11, xs12,
00954 xs21, xs22,
00955 valueT1, valueT2,
00956 k, theta );
00957 }
00958
00959 return sigma;
00960 }
00961
00962
00963
00964 G4double G4ecpssrBaseLixsModel::LinLinInterpolate(G4double e1,
00965 G4double e2,
00966 G4double e,
00967 G4double xs1,
00968 G4double xs2)
00969 {
00970 G4double value = xs1 + (xs2 - xs1)*(e - e1)/ (e2 - e1);
00971 return value;
00972 }
00973
00974
00975
00976 G4double G4ecpssrBaseLixsModel::LinLogInterpolate(G4double e1,
00977 G4double e2,
00978 G4double e,
00979 G4double xs1,
00980 G4double xs2)
00981 {
00982 G4double d1 = std::log(xs1);
00983 G4double d2 = std::log(xs2);
00984 G4double value = std::exp(d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
00985 return value;
00986 }
00987
00988
00989
00990 G4double G4ecpssrBaseLixsModel::LogLogInterpolate(G4double e1,
00991 G4double e2,
00992 G4double e,
00993 G4double xs1,
00994 G4double xs2)
00995 {
00996 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
00997 G4double b = std::log10(xs2) - a*std::log10(e2);
00998 G4double sigma = a*std::log10(e) + b;
00999 G4double value = (std::pow(10.,sigma));
01000 return value;
01001 }
01002
01003
01004
01005 G4double G4ecpssrBaseLixsModel::QuadInterpolator(G4double e11, G4double e12,
01006 G4double e21, G4double e22,
01007 G4double xs11, G4double xs12,
01008 G4double xs21, G4double xs22,
01009 G4double t1, G4double t2,
01010 G4double t, G4double e)
01011 {
01012
01013 G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
01014 G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
01015 G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
01016
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030 return value;
01031
01032 }
01033