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
00045
00046
00047
00048
00049
00050
00051
00052
00053 #include "G4eIonisationSpectrum.hh"
00054 #include "G4AtomicTransitionManager.hh"
00055 #include "G4AtomicShell.hh"
00056 #include "G4DataVector.hh"
00057 #include "Randomize.hh"
00058 #include "G4PhysicalConstants.hh"
00059 #include "G4SystemOfUnits.hh"
00060
00061
00062 G4eIonisationSpectrum::G4eIonisationSpectrum():G4VEnergySpectrum(),
00063 lowestE(0.1*eV),
00064 factor(1.3),
00065 iMax(24),
00066 verbose(0)
00067 {
00068 theParam = new G4eIonisationParameters();
00069 }
00070
00071
00072 G4eIonisationSpectrum::~G4eIonisationSpectrum()
00073 {
00074 delete theParam;
00075 }
00076
00077
00078 G4double G4eIonisationSpectrum::Probability(G4int Z,
00079 G4double tMin,
00080 G4double tMax,
00081 G4double e,
00082 G4int shell,
00083 const G4ParticleDefinition* ) const
00084 {
00085
00086
00087
00088
00089 G4double eMax = MaxEnergyOfSecondaries(e);
00090 G4double t0 = std::max(tMin, lowestE);
00091 G4double tm = std::min(tMax, eMax);
00092 if(t0 >= tm) return 0.0;
00093
00094 G4double bindingEnergy = (G4AtomicTransitionManager::Instance())->
00095 Shell(Z, shell)->BindingEnergy();
00096
00097 if(e <= bindingEnergy) return 0.0;
00098
00099 G4double energy = e + bindingEnergy;
00100
00101 G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
00102 G4double x2 = std::min(0.5,(tm + bindingEnergy)/energy);
00103
00104 if(verbose > 1 || (Z==4 && e>= 1.0 && e<= 0.0)) {
00105 G4cout << "G4eIonisationSpectrum::Probability: Z= " << Z
00106 << "; shell= " << shell
00107 << "; E(keV)= " << e/keV
00108 << "; Eb(keV)= " << bindingEnergy/keV
00109 << "; x1= " << x1
00110 << "; x2= " << x2
00111 << G4endl;
00112
00113 }
00114
00115 G4DataVector p;
00116
00117
00118 for (G4int i=0; i<iMax; i++)
00119 {
00120 G4double x = theParam->Parameter(Z, shell, i, e);
00121 if(i<4) x /= energy;
00122 p.push_back(x);
00123 }
00124
00125 if(p[3] > 0.5) p[3] = 0.5;
00126
00127 G4double gLocal = energy/electron_mass_c2 + 1.;
00128 p.push_back((2.0*gLocal - 1.0)/(gLocal*gLocal));
00129
00130
00131
00132 if (p[3] > 0)
00133 p[iMax-1] = Function(p[3], p);
00134 else
00135 {
00136 G4cout << "WARNING: G4eIonisationSpectrum::Probability "
00137 << "parameter p[3] <= 0. G4LEDATA dabatase might be corrupted for Z = "
00138 << Z << ". Please check and/or update it " << G4endl;
00139 }
00140
00141 if(e >= 1. && e <= 0. && Z == 4) p.push_back(0.0);
00142
00143
00144 G4double val = IntSpectrum(x1, x2, p);
00145 G4double x0 = (lowestE + bindingEnergy)/energy;
00146 G4double nor = IntSpectrum(x0, 0.5, p);
00147
00148 if(verbose > 1 || (Z==4 && e>= 1.0 && e<= 0.0)) {
00149 G4cout << "tcut= " << tMin
00150 << "; tMax= " << tMax
00151 << "; x0= " << x0
00152 << "; x1= " << x1
00153 << "; x2= " << x2
00154 << "; val= " << val
00155 << "; nor= " << nor
00156 << "; sum= " << p[0]
00157 << "; a= " << p[1]
00158 << "; b= " << p[2]
00159 << "; c= " << p[3]
00160 << G4endl;
00161 if(shell == 1) G4cout << "============" << G4endl;
00162 }
00163
00164 p.clear();
00165
00166 if(nor > 0.0) val /= nor;
00167 else val = 0.0;
00168
00169 return val;
00170 }
00171
00172
00173 G4double G4eIonisationSpectrum::AverageEnergy(G4int Z,
00174 G4double tMin,
00175 G4double tMax,
00176 G4double e,
00177 G4int shell,
00178 const G4ParticleDefinition* ) const
00179 {
00180
00181
00182
00183
00184 G4double eMax = MaxEnergyOfSecondaries(e);
00185 G4double t0 = std::max(tMin, lowestE);
00186 G4double tm = std::min(tMax, eMax);
00187 if(t0 >= tm) return 0.0;
00188
00189 G4double bindingEnergy = (G4AtomicTransitionManager::Instance())->
00190 Shell(Z, shell)->BindingEnergy();
00191
00192 if(e <= bindingEnergy) return 0.0;
00193
00194 G4double energy = e + bindingEnergy;
00195
00196 G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
00197 G4double x2 = std::min(0.5,(tm + bindingEnergy)/energy);
00198
00199 if(verbose > 1) {
00200 G4cout << "G4eIonisationSpectrum::AverageEnergy: Z= " << Z
00201 << "; shell= " << shell
00202 << "; E(keV)= " << e/keV
00203 << "; bindingE(keV)= " << bindingEnergy/keV
00204 << "; x1= " << x1
00205 << "; x2= " << x2
00206 << G4endl;
00207 }
00208
00209 G4DataVector p;
00210
00211
00212 for (G4int i=0; i<iMax; i++)
00213 {
00214 G4double x = theParam->Parameter(Z, shell, i, e);
00215 if(i<4) x /= energy;
00216 p.push_back(x);
00217 }
00218
00219 if(p[3] > 0.5) p[3] = 0.5;
00220
00221 G4double gLocal2 = energy/electron_mass_c2 + 1.;
00222 p.push_back((2.0*gLocal2 - 1.0)/(gLocal2*gLocal2));
00223
00224
00225
00226
00227 if (p[3] > 0)
00228 p[iMax-1] = Function(p[3], p);
00229 else
00230 {
00231 G4cout << "WARNING: G4eIonisationSpectrum::AverageEnergy "
00232 << "parameter p[3] <= 0. G4LEDATA dabatase might be corrupted for Z = "
00233 << Z << ". Please check and/or update it " << G4endl;
00234 }
00235
00236 G4double val = AverageValue(x1, x2, p);
00237 G4double x0 = (lowestE + bindingEnergy)/energy;
00238 G4double nor = IntSpectrum(x0, 0.5, p);
00239 val *= energy;
00240
00241 if(verbose > 1) {
00242 G4cout << "tcut(MeV)= " << tMin/MeV
00243 << "; tMax(MeV)= " << tMax/MeV
00244 << "; x0= " << x0
00245 << "; x1= " << x1
00246 << "; x2= " << x2
00247 << "; val= " << val
00248 << "; nor= " << nor
00249 << "; sum= " << p[0]
00250 << "; a= " << p[1]
00251 << "; b= " << p[2]
00252 << "; c= " << p[3]
00253 << G4endl;
00254 }
00255
00256 p.clear();
00257
00258 if(nor > 0.0) val /= nor;
00259 else val = 0.0;
00260
00261 return val;
00262 }
00263
00264
00265 G4double G4eIonisationSpectrum::SampleEnergy(G4int Z,
00266 G4double tMin,
00267 G4double tMax,
00268 G4double e,
00269 G4int shell,
00270 const G4ParticleDefinition* ) const
00271 {
00272
00273 G4double tDelta = 0.0;
00274 G4double t0 = std::max(tMin, lowestE);
00275 G4double tm = std::min(tMax, MaxEnergyOfSecondaries(e));
00276 if(t0 > tm) return tDelta;
00277
00278 G4double bindingEnergy = (G4AtomicTransitionManager::Instance())->
00279 Shell(Z, shell)->BindingEnergy();
00280
00281 if(e <= bindingEnergy) return 0.0;
00282
00283 G4double energy = e + bindingEnergy;
00284
00285 G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
00286 G4double x2 = std::min(0.5,(tm + bindingEnergy)/energy);
00287 if(x1 >= x2) return tDelta;
00288
00289 if(verbose > 1) {
00290 G4cout << "G4eIonisationSpectrum::SampleEnergy: Z= " << Z
00291 << "; shell= " << shell
00292 << "; E(keV)= " << e/keV
00293 << G4endl;
00294 }
00295
00296
00297 G4DataVector p;
00298
00299
00300 for (G4int i=0; i<iMax; i++)
00301 {
00302 G4double x = theParam->Parameter(Z, shell, i, e);
00303 if(i<4) x /= energy;
00304 p.push_back(x);
00305 }
00306
00307 if(p[3] > 0.5) p[3] = 0.5;
00308
00309 G4double gLocal3 = energy/electron_mass_c2 + 1.;
00310 p.push_back((2.0*gLocal3 - 1.0)/(gLocal3*gLocal3));
00311
00312
00313
00314
00315 if (p[3] > 0)
00316 p[iMax-1] = Function(p[3], p);
00317 else
00318 {
00319 G4cout << "WARNING: G4eIonisationSpectrum::SampleSpectrum "
00320 << "parameter p[3] <= 0. G4LEDATA dabatase might be corrupted for Z = "
00321 << Z << ". Please check and/or update it " << G4endl;
00322 }
00323
00324 G4double aria1 = 0.0;
00325 G4double a1 = std::max(x1,p[1]);
00326 G4double a2 = std::min(x2,p[3]);
00327 if(a1 < a2) aria1 = IntSpectrum(a1, a2, p);
00328 G4double aria2 = 0.0;
00329 G4double a3 = std::max(x1,p[3]);
00330 G4double a4 = x2;
00331 if(a3 < a4) aria2 = IntSpectrum(a3, a4, p);
00332
00333 G4double aria = (aria1 + aria2)*G4UniformRand();
00334 G4double amaj, fun, q, x, z1, z2, dx, dx1;
00335
00336
00337
00338 if(aria <= aria1) {
00339
00340 amaj = p[4];
00341 for (G4int j=5; j<iMax; j++) {
00342 if(p[j] > amaj) amaj = p[j];
00343 }
00344
00345 a1 = 1./a1;
00346 a2 = 1./a2;
00347
00348 G4int i;
00349 do {
00350
00351 x = 1./(a2 + G4UniformRand()*(a1 - a2));
00352 z1 = p[1];
00353 z2 = p[3];
00354 dx = (p[2] - p[1]) / 3.0;
00355 dx1= std::exp(std::log(p[3]/p[2]) / 16.0);
00356 for (i=4; i<iMax-1; i++) {
00357
00358 if (i < 7) {
00359 z2 = z1 + dx;
00360 } else if(iMax-2 == i) {
00361 z2 = p[3];
00362 break;
00363 } else {
00364 z2 = z1*dx1;
00365 }
00366 if(x >= z1 && x <= z2) break;
00367 z1 = z2;
00368 }
00369 fun = p[i] + (x - z1) * (p[i+1] - p[i])/(z2 - z1);
00370
00371 if(fun > amaj) {
00372 G4cout << "WARNING in G4eIonisationSpectrum::SampleEnergy:"
00373 << " Majoranta " << amaj
00374 << " < " << fun
00375 << " in the first aria at x= " << x
00376 << G4endl;
00377 }
00378
00379 q = amaj*G4UniformRand();
00380
00381 } while (q >= fun);
00382
00383
00384
00385 } else {
00386
00387 amaj = std::max(p[iMax-1], Function(0.5, p)) * factor;
00388 a1 = 1./a3;
00389 a2 = 1./a4;
00390
00391 do {
00392
00393 x = 1./(a2 + G4UniformRand()*(a1 - a2));
00394 fun = Function(x, p);
00395
00396 if(fun > amaj) {
00397 G4cout << "WARNING in G4eIonisationSpectrum::SampleEnergy:"
00398 << " Majoranta " << amaj
00399 << " < " << fun
00400 << " in the second aria at x= " << x
00401 << G4endl;
00402 }
00403
00404 q = amaj*G4UniformRand();
00405
00406 } while (q >= fun);
00407
00408 }
00409
00410 p.clear();
00411
00412 tDelta = x*energy - bindingEnergy;
00413
00414 if(verbose > 1) {
00415 G4cout << "tcut(MeV)= " << tMin/MeV
00416 << "; tMax(MeV)= " << tMax/MeV
00417 << "; x1= " << x1
00418 << "; x2= " << x2
00419 << "; a1= " << a1
00420 << "; a2= " << a2
00421 << "; x= " << x
00422 << "; be= " << bindingEnergy
00423 << "; e= " << e
00424 << "; tDelta= " << tDelta
00425 << G4endl;
00426 }
00427
00428
00429 return tDelta;
00430 }
00431
00432
00433 G4double G4eIonisationSpectrum::IntSpectrum(G4double xMin,
00434 G4double xMax,
00435 const G4DataVector& p) const
00436 {
00437
00438 G4double sum = 0.0;
00439 if(xMin >= xMax) return sum;
00440
00441 G4double x1, x2, xs1, xs2, y1, y2, ys1, ys2, q;
00442
00443
00444 if(xMin < p[3]) {
00445
00446 x1 = p[1];
00447 y1 = p[4];
00448
00449 G4double dx = (p[2] - p[1]) / 3.0;
00450 G4double dx1= std::exp(std::log(p[3]/p[2]) / 16.0);
00451
00452 for (size_t i=0; i<19; i++) {
00453
00454 q = 0.0;
00455 if (i < 3) {
00456 x2 = x1 + dx;
00457 } else if(18 == i) {
00458 x2 = p[3];
00459 } else {
00460 x2 = x1*dx1;
00461 }
00462
00463 y2 = p[5 + i];
00464
00465 if (xMax <= x1) {
00466 break;
00467 } else if (xMin < x2) {
00468
00469 xs1 = x1;
00470 xs2 = x2;
00471 ys1 = y1;
00472 ys2 = y2;
00473
00474 if (x2 > x1) {
00475 if (xMin > x1) {
00476 xs1 = xMin;
00477 ys1 += (xs1 - x1)*(y2 - y1)/(x2 - x1);
00478 }
00479 if (xMax < x2) {
00480 xs2 = xMax;
00481 ys2 += (xs2 - x2)*(y1 - y2)/(x1 - x2);
00482 }
00483 if (xs2 > xs1) {
00484 q = (ys1*xs2 - ys2*xs1)/(xs1*xs2)
00485 + std::log(xs2/xs1)*(ys2 - ys1)/(xs2 - xs1);
00486 sum += q;
00487 if(p.size() == 26) G4cout << "i= " << i << " q= " << q << " sum= " << sum << G4endl;
00488 }
00489 }
00490 }
00491 x1 = x2;
00492 y1 = y2;
00493 }
00494 }
00495
00496
00497
00498 x1 = std::max(xMin, p[3]);
00499 if(x1 >= xMax) return sum;
00500 x2 = xMax;
00501
00502 xs1 = 1./x1;
00503 xs2 = 1./x2;
00504 q = (xs1 - xs2)*(1.0 - p[0])
00505 - p[iMax]*std::log(x2/x1)
00506 + (1. - p[iMax])*(x2 - x1)
00507 + 1./(1. - x2) - 1./(1. - x1)
00508 + p[iMax]*std::log((1. - x2)/(1. - x1))
00509 + 0.25*p[0]*(xs1*xs1 - xs2*xs2);
00510 sum += q;
00511 if(p.size() == 26) G4cout << "param... q= " << q << " sum= " << sum << G4endl;
00512
00513 return sum;
00514 }
00515
00516
00517 G4double G4eIonisationSpectrum::AverageValue(G4double xMin,
00518 G4double xMax,
00519 const G4DataVector& p) const
00520 {
00521 G4double sum = 0.0;
00522 if(xMin >= xMax) return sum;
00523
00524 G4double x1, x2, xs1, xs2, y1, y2, ys1, ys2;
00525
00526
00527 if(xMin < p[3]) {
00528
00529 x1 = p[1];
00530 y1 = p[4];
00531
00532 G4double dx = (p[2] - p[1]) / 3.0;
00533 G4double dx1= std::exp(std::log(p[3]/p[2]) / 16.0);
00534
00535 for (size_t i=0; i<19; i++) {
00536
00537 if (i < 3) {
00538 x2 = x1 + dx;
00539 } else if(18 == i) {
00540 x2 = p[3];
00541 } else {
00542 x2 = x1*dx1;
00543 }
00544
00545 y2 = p[5 + i];
00546
00547 if (xMax <= x1) {
00548 break;
00549 } else if (xMin < x2) {
00550
00551 xs1 = x1;
00552 xs2 = x2;
00553 ys1 = y1;
00554 ys2 = y2;
00555
00556 if (x2 > x1) {
00557 if (xMin > x1) {
00558 xs1 = xMin;
00559 ys1 += (xs1 - x1)*(y2 - y1)/(x2 - x1);
00560 }
00561 if (xMax < x2) {
00562 xs2 = xMax;
00563 ys2 += (xs2 - x2)*(y1 - y2)/(x1 - x2);
00564 }
00565 if (xs2 > xs1) {
00566 sum += std::log(xs2/xs1)*(ys1*xs2 - ys2*xs1)/(xs2 - xs1)
00567 + ys2 - ys1;
00568 }
00569 }
00570 }
00571 x1 = x2;
00572 y1 = y2;
00573
00574 }
00575 }
00576
00577
00578
00579 x1 = std::max(xMin, p[3]);
00580 if(x1 >= xMax) return sum;
00581 x2 = xMax;
00582
00583 xs1 = 1./x1;
00584 xs2 = 1./x2;
00585
00586 sum += std::log(x2/x1)*(1.0 - p[0])
00587 + 0.5*(1. - p[iMax])*(x2*x2 - x1*x1)
00588 + 1./(1. - x2) - 1./(1. - x1)
00589 + (1. + p[iMax])*std::log((1. - x2)/(1. - x1))
00590 + 0.5*p[0]*(xs1 - xs2);
00591
00592 return sum;
00593 }
00594
00595
00596 void G4eIonisationSpectrum::PrintData() const
00597 {
00598 theParam->PrintData();
00599 }
00600
00601 G4double G4eIonisationSpectrum::MaxEnergyOfSecondaries(G4double kineticEnergy,
00602 G4int,
00603 const G4ParticleDefinition* ) const
00604 {
00605 return 0.5 * kineticEnergy;
00606 }