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
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044 #include "G4OrlicLiXsModel.hh"
00045
00046 #include "globals.hh"
00047 #include "G4PhysicalConstants.hh"
00048 #include "G4SystemOfUnits.hh"
00049 #include "G4Proton.hh"
00050
00051
00052 G4OrlicLiXsModel::G4OrlicLiXsModel()
00053 {
00054
00055 transitionManager = G4AtomicTransitionManager::Instance();
00056
00057 }
00058
00059 G4OrlicLiXsModel::~G4OrlicLiXsModel()
00060 {
00061
00062 }
00063
00064
00065
00066
00067
00068
00069
00070 G4double G4OrlicLiXsModel::CalculateL1CrossSection(G4int zTarget, G4double energyIncident)
00071
00072 {
00073
00074 if ( zTarget < 41 )
00075
00076 {
00077
00078
00079
00080 return 0;
00081 }
00082
00083
00084
00085
00086
00087
00088
00089
00090 G4double l1BindingEnergy = transitionManager->Shell(zTarget,1)->BindingEnergy()/keV;
00091
00092
00093 G4double lamda = 1836.109;
00094
00095 G4double normalizedEnergy = (energyIncident/keV)/(lamda*l1BindingEnergy);
00096
00097 G4double x = std::log(normalizedEnergy);
00098
00099 G4double a0 = 0.;
00100 G4double a1 = 0.;
00101 G4double a2 = 0.;
00102 G4double a3 = 0.;
00103 G4double a4 = 0.;
00104 G4double a5 = 0.;
00105 G4double a6 = 0.;
00106 G4double a7 = 0.;
00107 G4double a8 = 0.;
00108 G4double a9 = 0.;
00109
00110
00111 if ( (zTarget>=41 && zTarget<=50) && (normalizedEnergy>=0.013 && normalizedEnergy<=1) )
00112 {
00113
00114 G4cout << "Energy1 (keV) = " << normalizedEnergy * lamda*l1BindingEnergy << G4endl;
00115
00116 a0=11.274881;
00117 a1=-0.187401;
00118 a2=-0.943341;
00119 a3=-1.47817;
00120 a4=-1.282343;
00121 a5=-0.386544;
00122 a6=-0.037932;
00123 a7=0.;
00124 a8=0.;
00125 a9=0.;
00126 }
00127
00128 else if ( (zTarget>=51 && zTarget<=60) && (normalizedEnergy>=0.012 && normalizedEnergy<=0.95))
00129 {
00130
00131
00132
00133 a0=11.242637;
00134 a1=-0.162515;
00135 a2=1.035774;
00136 a3=3.970908;
00137 a4=3.968233;
00138 a5=1.655714;
00139 a6=0.058885;
00140 a7=-0.155743;
00141 a8=-0.042228;
00142 a9=-0.003371;
00143 }
00144
00145 else if ( (zTarget>=61 && zTarget<=70) && (normalizedEnergy>=0.01 && normalizedEnergy<=0.6) )
00146 {
00147
00148
00149
00150 a0=6.476722;
00151 a1=-25.804787;
00152 a2=-54.061629;
00153 a3=-56.684589;
00154 a4=-33.223367;
00155 a5=-11.034979;
00156 a6=-2.042851;
00157 a7=-0.194075;
00158 a8=-0.007252;
00159 a9=0.;
00160 }
00161 else if ( (zTarget>=71 && zTarget<=80) && (normalizedEnergy>=0.01 && normalizedEnergy<=0.45) )
00162 {
00163
00164
00165
00166 a0=12.776794;
00167 a1=6.562907;
00168 a2=10.158703;
00169 a3=7.432592;
00170 a4=2.332036;
00171 a5=0.317946;
00172 a6=0.014479;
00173 a7=0.;
00174 a8=0.;
00175 a9=0.;
00176 }
00177 else if ( (zTarget>=81 && zTarget<=92) && (normalizedEnergy>=0.008 && normalizedEnergy<=0.3) )
00178 {
00179
00180
00181
00182 a0=28.243087;
00183 a1=50.199585;
00184 a2=58.281684;
00185 a3=34.130538;
00186 a4=10.268531;
00187 a5=1.525302;
00188 a6=0.08835;
00189 a7=0.;
00190 a8=0.;
00191 a9=0.;
00192 }
00193 else {return 0;}
00194
00195
00196 G4double analyticalFunction = a0 + (a1*x)+(a2*x*x)+(a3*std::pow(x,3))+(a4*std::pow(x,4))+(a5*std::pow(x,5))+(a6*std::pow(x,6))+
00197 (a7*std::pow(x,7))+(a8*std::pow(x,8))+(a9*std::pow(x,9));
00198
00199
00200
00201 G4double L1crossSection = std::exp(analyticalFunction)/(l1BindingEnergy*l1BindingEnergy);
00202
00203
00204 if (L1crossSection >= 0) {
00205 return L1crossSection * barn;
00206 }
00207 else {return 0;}
00208
00209 }
00210
00211
00212
00213
00214 G4double G4OrlicLiXsModel::CalculateL2CrossSection(G4int zTarget, G4double energyIncident)
00215
00216 {
00217
00218
00219 if ( zTarget < 41)
00220
00221 {
00222
00223
00224
00225 return 0;
00226 }
00227
00228
00229
00230 G4double massIncident;
00231
00232 G4Proton* aProtone = G4Proton::Proton();
00233
00234 massIncident = aProtone->GetPDGMass();
00235
00236 G4double L2crossSection;
00237
00238 G4double l2BindingEnergy = (transitionManager->Shell(zTarget,2)->BindingEnergy())/keV;
00239
00240 G4double lamda = massIncident/electron_mass_c2;
00241
00242 G4double normalizedEnergy = (energyIncident/keV)/(lamda*l2BindingEnergy);
00243
00244 G4double x = std::log(normalizedEnergy);
00245
00246 G4double a0 = 0.;
00247 G4double a1 = 0.;
00248 G4double a2 = 0.;
00249 G4double a3 = 0.;
00250 G4double a4 = 0.;
00251 G4double a5 = 0.;
00252
00253 if ( (zTarget>=41 && zTarget<=50) && (normalizedEnergy>=0.015 && normalizedEnergy<=1.5))
00254 {
00255 a0=11.194798;
00256 a1=0.178807;
00257 a2=-0.449865;
00258 a3=-0.063528;
00259 a4=-0.015364;
00260 a5=0.;
00261 }
00262
00263 else if ( (zTarget>=51 && zTarget<=60) && (normalizedEnergy>=0.012 && normalizedEnergy<=1.0))
00264 {
00265 a0=11.241409;
00266 a1=0.149635;
00267 a2=-0.633269;
00268 a3=-0.17834;
00269 a4=-0.034743;
00270 a5=0.006474;
00271
00272 }
00273
00274 else if ( (zTarget>=61 && zTarget<=70) && (normalizedEnergy>=0.01 && normalizedEnergy<=0.65))
00275 {
00276 a0=11.247424;
00277 a1=0.203051;
00278 a2=-0.219083;
00279 a3=0.164514;
00280 a4=0.058692;
00281 a5=0.007866;
00282 }
00283 else if ( (zTarget>=71 && zTarget<=80) && (normalizedEnergy>=0.01 && normalizedEnergy<=0.47))
00284 {
00285 a0=11.229924;
00286 a1=-0.087241;
00287 a2=-0.753908;
00288 a3=-0.181546;
00289 a4=-0.030406;
00290 a5=0.;
00291 }
00292 else if ( (zTarget>=81 && zTarget<=92) && (normalizedEnergy>=0.01 && normalizedEnergy<=0.35))
00293 {
00294 a0=11.586671;
00295 a1=0.730838;
00296 a2=-0.056713;
00297 a3=0.053262;
00298 a4=-0.003672;
00299 a5=0.;
00300 }
00301 else {return 0;}
00302
00303 G4double analyticalFunction = a0 + (a1*x)+(a2*x*x)+(a3*std::pow(x,3))+(a4*std::pow(x,4))+(a5*std::pow(x,5));
00304
00305
00306 L2crossSection = std::exp(analyticalFunction)/(l2BindingEnergy*l2BindingEnergy);
00307
00308
00309
00310 if (L2crossSection >= 0) {
00311 return L2crossSection * barn;
00312 }
00313 else {return 0;}
00314
00315 }
00316
00317
00318
00319
00320 G4double G4OrlicLiXsModel::CalculateL3CrossSection(G4int zTarget, G4double energyIncident)
00321
00322 {
00323
00324 if ( zTarget < 41)
00325
00326 {
00327
00328
00329
00330 return 0;
00331 }
00332
00333 G4double massIncident;
00334
00335 G4Proton* aProtone = G4Proton::Proton();
00336
00337 massIncident = aProtone->GetPDGMass();
00338
00339
00340 G4double L3crossSection;
00341
00342 G4double l3BindingEnergy = (transitionManager->Shell(zTarget,3)->BindingEnergy())/keV;
00343
00344
00345 G4double lamda = massIncident/electron_mass_c2;
00346
00347 G4double normalizedEnergy = (energyIncident/keV)/(lamda*l3BindingEnergy);
00348
00349 G4double x = std::log(normalizedEnergy);
00350
00351
00352 G4double a0 = 0.;
00353 G4double a1 = 0.;
00354 G4double a2 = 0.;
00355 G4double a3 = 0.;
00356 G4double a4 = 0.;
00357 G4double a5 = 0.;
00358
00359 if ( (zTarget>=41 && zTarget<=50 ) && (normalizedEnergy>=0.015 && normalizedEnergy<=1.5))
00360 {
00361 a0=11.91837;
00362 a1=0.03064;
00363 a2=-0.657644;
00364 a3=-0.14532;
00365 a4=-0.026059;
00366
00367
00368 }
00369
00370 else if ( (zTarget>=51 && zTarget<=60 ) && (normalizedEnergy>=0.013 && normalizedEnergy<=1.1))
00371 {
00372 a0=11.909485;
00373 a1=0.15918;
00374 a2=-0.588004;
00375 a3=-0.159466;
00376 a4=-0.033184;
00377 }
00378
00379 else if ( (zTarget>=61 && zTarget<=70 ) && (normalizedEnergy>=0.01 && normalizedEnergy<=0.67))
00380 {
00381 a0=11.878472;
00382 a1=-0.137007;
00383 a2=-0.959475;
00384 a3=-0.316505;
00385 a4=-0.054154;
00386 }
00387 else if ( (zTarget>=71 && zTarget<=80 ) && (normalizedEnergy>=0.013 && normalizedEnergy<=0.5))
00388 {
00389 a0=11.802538;
00390 a1=-0.371796;
00391 a2=-1.052238;
00392 a3=-0.28766;
00393 a4=-0.042608;
00394 }
00395 else if ( (zTarget>=81 && zTarget<=92 ) && (normalizedEnergy>=0.01 && normalizedEnergy<=0.35))
00396 {
00397 a0=11.423712;
00398 a1=-1.428823;
00399 a2=-1.946979;
00400 a3=-0.585198;
00401 a4=-0.076467;
00402 }
00403 else {return 0;}
00404
00405
00406 G4double analyticalFunction = a0 + (a1*x)+(a2*x*x)+(a3*std::pow(x,3))+(a4*std::pow(x,4))+(a5*std::pow(x,5));
00407
00408
00409 L3crossSection = std::exp(analyticalFunction)/(l3BindingEnergy*l3BindingEnergy);
00410
00411
00412 if (L3crossSection >= 0) {
00413 return L3crossSection * barn;
00414 }
00415 else {return 0;}
00416
00417
00418 }