File indexing completed on 2025-01-31 09:22:07
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050 #include "G4RDeIonisationSpectrum.hh"
0051 #include "G4RDAtomicTransitionManager.hh"
0052 #include "G4RDAtomicShell.hh"
0053 #include "G4PhysicalConstants.hh"
0054 #include "G4SystemOfUnits.hh"
0055 #include "G4DataVector.hh"
0056 #include "Randomize.hh"
0057
0058
0059 G4RDeIonisationSpectrum::G4RDeIonisationSpectrum():G4RDVEnergySpectrum(),
0060 lowestE(0.1*eV),
0061 factor(1.3),
0062 iMax(24),
0063 verbose(0)
0064 {
0065 theParam = new G4RDeIonisationParameters();
0066 }
0067
0068
0069 G4RDeIonisationSpectrum::~G4RDeIonisationSpectrum()
0070 {
0071 delete theParam;
0072 }
0073
0074
0075 G4double G4RDeIonisationSpectrum::Probability(G4int Z,
0076 G4double tMin,
0077 G4double tMax,
0078 G4double e,
0079 G4int shell,
0080 const G4ParticleDefinition* ) const
0081 {
0082
0083
0084
0085
0086 G4double eMax = MaxEnergyOfSecondaries(e);
0087 G4double t0 = std::max(tMin, lowestE);
0088 G4double tm = std::min(tMax, eMax);
0089 if(t0 >= tm) return 0.0;
0090
0091 G4double bindingEnergy = (G4RDAtomicTransitionManager::Instance())->
0092 Shell(Z, shell)->BindingEnergy();
0093
0094 if(e <= bindingEnergy) return 0.0;
0095
0096 G4double energy = e + bindingEnergy;
0097
0098 G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
0099 G4double x2 = std::min(0.5,(tm + bindingEnergy)/energy);
0100
0101 if(verbose > 1 || (Z==4 && e>= 1.0 && e<= 0.0)) {
0102 G4cout << "G4RDeIonisationSpectrum::Probability: Z= " << Z
0103 << "; shell= " << shell
0104 << "; E(keV)= " << e/keV
0105 << "; Eb(keV)= " << bindingEnergy/keV
0106 << "; x1= " << x1
0107 << "; x2= " << x2
0108 << G4endl;
0109
0110 }
0111
0112 G4DataVector p;
0113
0114
0115 for (G4int i=0; i<iMax; i++)
0116 {
0117 G4double x = theParam->Parameter(Z, shell, i, e);
0118 if(i<4) x /= energy;
0119 p.push_back(x);
0120 }
0121
0122 if(p[3] > 0.5) p[3] = 0.5;
0123
0124 G4double g = energy/electron_mass_c2 + 1.;
0125 p.push_back((2.0*g - 1.0)/(g*g));
0126
0127 p[iMax-1] = Function(p[3], p);
0128
0129 if(e >= 1. && e <= 0. && Z == 4) p.push_back(0.0);
0130
0131
0132 G4double val = IntSpectrum(x1, x2, p);
0133 G4double x0 = (lowestE + bindingEnergy)/energy;
0134 G4double nor = IntSpectrum(x0, 0.5, p);
0135
0136 if(verbose > 1 || (Z==4 && e>= 1.0 && e<= 0.0)) {
0137 G4cout << "tcut= " << tMin
0138 << "; tMax= " << tMax
0139 << "; x0= " << x0
0140 << "; x1= " << x1
0141 << "; x2= " << x2
0142 << "; val= " << val
0143 << "; nor= " << nor
0144 << "; sum= " << p[0]
0145 << "; a= " << p[1]
0146 << "; b= " << p[2]
0147 << "; c= " << p[3]
0148 << G4endl;
0149 if(shell == 1) G4cout << "============" << G4endl;
0150 }
0151
0152 p.clear();
0153
0154 if(nor > 0.0) val /= nor;
0155 else val = 0.0;
0156
0157 return val;
0158 }
0159
0160
0161 G4double G4RDeIonisationSpectrum::AverageEnergy(G4int Z,
0162 G4double tMin,
0163 G4double tMax,
0164 G4double e,
0165 G4int shell,
0166 const G4ParticleDefinition* ) const
0167 {
0168
0169
0170
0171
0172 G4double eMax = MaxEnergyOfSecondaries(e);
0173 G4double t0 = std::max(tMin, lowestE);
0174 G4double tm = std::min(tMax, eMax);
0175 if(t0 >= tm) return 0.0;
0176
0177 G4double bindingEnergy = (G4RDAtomicTransitionManager::Instance())->
0178 Shell(Z, shell)->BindingEnergy();
0179
0180 if(e <= bindingEnergy) return 0.0;
0181
0182 G4double energy = e + bindingEnergy;
0183
0184 G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
0185 G4double x2 = std::min(0.5,(tm + bindingEnergy)/energy);
0186
0187 if(verbose > 1) {
0188 G4cout << "G4RDeIonisationSpectrum::AverageEnergy: Z= " << Z
0189 << "; shell= " << shell
0190 << "; E(keV)= " << e/keV
0191 << "; bindingE(keV)= " << bindingEnergy/keV
0192 << "; x1= " << x1
0193 << "; x2= " << x2
0194 << G4endl;
0195 }
0196
0197 G4DataVector p;
0198
0199
0200 for (G4int i=0; i<iMax; i++)
0201 {
0202 G4double x = theParam->Parameter(Z, shell, i, e);
0203 if(i<4) x /= energy;
0204 p.push_back(x);
0205 }
0206
0207 if(p[3] > 0.5) p[3] = 0.5;
0208
0209 G4double g = energy/electron_mass_c2 + 1.;
0210 p.push_back((2.0*g - 1.0)/(g*g));
0211
0212 p[iMax-1] = Function(p[3], p);
0213
0214 G4double val = AverageValue(x1, x2, p);
0215 G4double x0 = (lowestE + bindingEnergy)/energy;
0216 G4double nor = IntSpectrum(x0, 0.5, p);
0217 val *= energy;
0218
0219 if(verbose > 1) {
0220 G4cout << "tcut(MeV)= " << tMin/MeV
0221 << "; tMax(MeV)= " << tMax/MeV
0222 << "; x0= " << x0
0223 << "; x1= " << x1
0224 << "; x2= " << x2
0225 << "; val= " << val
0226 << "; nor= " << nor
0227 << "; sum= " << p[0]
0228 << "; a= " << p[1]
0229 << "; b= " << p[2]
0230 << "; c= " << p[3]
0231 << G4endl;
0232 }
0233
0234 p.clear();
0235
0236 if(nor > 0.0) val /= nor;
0237 else val = 0.0;
0238
0239 return val;
0240 }
0241
0242
0243 G4double G4RDeIonisationSpectrum::SampleEnergy(G4int Z,
0244 G4double tMin,
0245 G4double tMax,
0246 G4double e,
0247 G4int shell,
0248 const G4ParticleDefinition* ) const
0249 {
0250
0251 G4double tDelta = 0.0;
0252 G4double t0 = std::max(tMin, lowestE);
0253 G4double tm = std::min(tMax, MaxEnergyOfSecondaries(e));
0254 if(t0 > tm) return tDelta;
0255
0256 G4double bindingEnergy = (G4RDAtomicTransitionManager::Instance())->
0257 Shell(Z, shell)->BindingEnergy();
0258
0259 if(e <= bindingEnergy) return 0.0;
0260
0261 G4double energy = e + bindingEnergy;
0262
0263 G4double x1 = std::min(0.5,(t0 + bindingEnergy)/energy);
0264 G4double x2 = std::min(0.5,(tm + bindingEnergy)/energy);
0265 if(x1 >= x2) return tDelta;
0266
0267 if(verbose > 1) {
0268 G4cout << "G4RDeIonisationSpectrum::SampleEnergy: Z= " << Z
0269 << "; shell= " << shell
0270 << "; E(keV)= " << e/keV
0271 << G4endl;
0272 }
0273
0274
0275 G4DataVector p;
0276
0277
0278 for (G4int i=0; i<iMax; i++)
0279 {
0280 G4double x = theParam->Parameter(Z, shell, i, e);
0281 if(i<4) x /= energy;
0282 p.push_back(x);
0283 }
0284
0285 if(p[3] > 0.5) p[3] = 0.5;
0286
0287 G4double g = energy/electron_mass_c2 + 1.;
0288 p.push_back((2.0*g - 1.0)/(g*g));
0289
0290 p[iMax-1] = Function(p[3], p);
0291
0292 G4double aria1 = 0.0;
0293 G4double a1 = std::max(x1,p[1]);
0294 G4double a2 = std::min(x2,p[3]);
0295 if(a1 < a2) aria1 = IntSpectrum(a1, a2, p);
0296 G4double aria2 = 0.0;
0297 G4double a3 = std::max(x1,p[3]);
0298 G4double a4 = x2;
0299 if(a3 < a4) aria2 = IntSpectrum(a3, a4, p);
0300
0301 G4double aria = (aria1 + aria2)*G4UniformRand();
0302 G4double amaj, fun, q, x, z1, z2, dx, dx1;
0303
0304
0305
0306 if(aria <= aria1) {
0307
0308 amaj = p[4];
0309 for (G4int j=5; j<iMax; j++) {
0310 if(p[j] > amaj) amaj = p[j];
0311 }
0312
0313 a1 = 1./a1;
0314 a2 = 1./a2;
0315
0316 G4int i;
0317 do {
0318
0319 x = 1./(a2 + G4UniformRand()*(a1 - a2));
0320 z1 = p[1];
0321 z2 = p[3];
0322 dx = (p[2] - p[1]) / 3.0;
0323 dx1= std::exp(std::log(p[3]/p[2]) / 16.0);
0324 for (i=4; i<iMax-1; i++) {
0325
0326 if (i < 7) {
0327 z2 = z1 + dx;
0328 } else if(iMax-2 == i) {
0329 z2 = p[3];
0330 break;
0331 } else {
0332 z2 = z1*dx1;
0333 }
0334 if(x >= z1 && x <= z2) break;
0335 z1 = z2;
0336 }
0337 fun = p[i] + (x - z1) * (p[i+1] - p[i])/(z2 - z1);
0338
0339 if(fun > amaj) {
0340 G4cout << "WARNING in G4RDeIonisationSpectrum::SampleEnergy:"
0341 << " Majoranta " << amaj
0342 << " < " << fun
0343 << " in the first aria at x= " << x
0344 << G4endl;
0345 }
0346
0347 q = amaj*G4UniformRand();
0348
0349 } while (q >= fun);
0350
0351
0352
0353 } else {
0354
0355 amaj = std::max(p[iMax-1], Function(0.5, p)) * factor;
0356 a1 = 1./a3;
0357 a2 = 1./a4;
0358
0359 do {
0360
0361 x = 1./(a2 + G4UniformRand()*(a1 - a2));
0362 fun = Function(x, p);
0363
0364 if(fun > amaj) {
0365 G4cout << "WARNING in G4RDeIonisationSpectrum::SampleEnergy:"
0366 << " Majoranta " << amaj
0367 << " < " << fun
0368 << " in the second aria at x= " << x
0369 << G4endl;
0370 }
0371
0372 q = amaj*G4UniformRand();
0373
0374 } while (q >= fun);
0375
0376 }
0377
0378 p.clear();
0379
0380 tDelta = x*energy - bindingEnergy;
0381
0382 if(verbose > 1) {
0383 G4cout << "tcut(MeV)= " << tMin/MeV
0384 << "; tMax(MeV)= " << tMax/MeV
0385 << "; x1= " << x1
0386 << "; x2= " << x2
0387 << "; a1= " << a1
0388 << "; a2= " << a2
0389 << "; x= " << x
0390 << "; be= " << bindingEnergy
0391 << "; e= " << e
0392 << "; tDelta= " << tDelta
0393 << G4endl;
0394 }
0395
0396
0397 return tDelta;
0398 }
0399
0400
0401 G4double G4RDeIonisationSpectrum::IntSpectrum(G4double xMin,
0402 G4double xMax,
0403 const G4DataVector& p) const
0404 {
0405
0406 G4double sum = 0.0;
0407 if(xMin >= xMax) return sum;
0408
0409 G4double x1, x2, xs1, xs2, y1, y2, ys1, ys2, q;
0410
0411
0412 if(xMin < p[3]) {
0413
0414 x1 = p[1];
0415 y1 = p[4];
0416
0417 G4double dx = (p[2] - p[1]) / 3.0;
0418 G4double dx1= std::exp(std::log(p[3]/p[2]) / 16.0);
0419
0420 for (size_t i=0; i<19; i++) {
0421
0422 q = 0.0;
0423 if (i < 3) {
0424 x2 = x1 + dx;
0425 } else if(18 == i) {
0426 x2 = p[3];
0427 } else {
0428 x2 = x1*dx1;
0429 }
0430
0431 y2 = p[5 + i];
0432
0433 if (xMax <= x1) {
0434 break;
0435 } else if (xMin < x2) {
0436
0437 xs1 = x1;
0438 xs2 = x2;
0439 ys1 = y1;
0440 ys2 = y2;
0441
0442 if (x2 > x1) {
0443 if (xMin > x1) {
0444 xs1 = xMin;
0445 ys1 += (xs1 - x1)*(y2 - y1)/(x2 - x1);
0446 }
0447 if (xMax < x2) {
0448 xs2 = xMax;
0449 ys2 += (xs2 - x2)*(y1 - y2)/(x1 - x2);
0450 }
0451 if (xs2 > xs1) {
0452 q = (ys1*xs2 - ys2*xs1)/(xs1*xs2)
0453 + std::log(xs2/xs1)*(ys2 - ys1)/(xs2 - xs1);
0454 sum += q;
0455 if(p.size() == 26) G4cout << "i= " << i << " q= " << q << " sum= " << sum << G4endl;
0456 }
0457 }
0458 }
0459 x1 = x2;
0460 y1 = y2;
0461 }
0462 }
0463
0464
0465
0466 x1 = std::max(xMin, p[3]);
0467 if(x1 >= xMax) return sum;
0468 x2 = xMax;
0469
0470 xs1 = 1./x1;
0471 xs2 = 1./x2;
0472 q = (xs1 - xs2)*(1.0 - p[0])
0473 - p[iMax]*std::log(x2/x1)
0474 + (1. - p[iMax])*(x2 - x1)
0475 + 1./(1. - x2) - 1./(1. - x1)
0476 + p[iMax]*std::log((1. - x2)/(1. - x1))
0477 + 0.25*p[0]*(xs1*xs1 - xs2*xs2);
0478 sum += q;
0479 if(p.size() == 26) G4cout << "param... q= " << q << " sum= " << sum << G4endl;
0480
0481 return sum;
0482 }
0483
0484
0485 G4double G4RDeIonisationSpectrum::AverageValue(G4double xMin,
0486 G4double xMax,
0487 const G4DataVector& p) const
0488 {
0489 G4double sum = 0.0;
0490 if(xMin >= xMax) return sum;
0491
0492 G4double x1, x2, xs1, xs2, y1, y2, ys1, ys2;
0493
0494
0495 if(xMin < p[3]) {
0496
0497 x1 = p[1];
0498 y1 = p[4];
0499
0500 G4double dx = (p[2] - p[1]) / 3.0;
0501 G4double dx1= std::exp(std::log(p[3]/p[2]) / 16.0);
0502
0503 for (size_t i=0; i<19; i++) {
0504
0505 if (i < 3) {
0506 x2 = x1 + dx;
0507 } else if(18 == i) {
0508 x2 = p[3];
0509 } else {
0510 x2 = x1*dx1;
0511 }
0512
0513 y2 = p[5 + i];
0514
0515 if (xMax <= x1) {
0516 break;
0517 } else if (xMin < x2) {
0518
0519 xs1 = x1;
0520 xs2 = x2;
0521 ys1 = y1;
0522 ys2 = y2;
0523
0524 if (x2 > x1) {
0525 if (xMin > x1) {
0526 xs1 = xMin;
0527 ys1 += (xs1 - x1)*(y2 - y1)/(x2 - x1);
0528 }
0529 if (xMax < x2) {
0530 xs2 = xMax;
0531 ys2 += (xs2 - x2)*(y1 - y2)/(x1 - x2);
0532 }
0533 if (xs2 > xs1) {
0534 sum += std::log(xs2/xs1)*(ys1*xs2 - ys2*xs1)/(xs2 - xs1)
0535 + ys2 - ys1;
0536 }
0537 }
0538 }
0539 x1 = x2;
0540 y1 = y2;
0541
0542 }
0543 }
0544
0545
0546
0547 x1 = std::max(xMin, p[3]);
0548 if(x1 >= xMax) return sum;
0549 x2 = xMax;
0550
0551 xs1 = 1./x1;
0552 xs2 = 1./x2;
0553
0554 sum += std::log(x2/x1)*(1.0 - p[0])
0555 + 0.5*(1. - p[iMax])*(x2*x2 - x1*x1)
0556 + 1./(1. - x2) - 1./(1. - x1)
0557 + (1. + p[iMax])*std::log((1. - x2)/(1. - x1))
0558 + 0.5*p[0]*(xs1 - xs2);
0559
0560 return sum;
0561 }
0562
0563
0564 void G4RDeIonisationSpectrum::PrintData() const
0565 {
0566 theParam->PrintData();
0567 }
0568
0569 G4double G4RDeIonisationSpectrum::MaxEnergyOfSecondaries(G4double kineticEnergy,
0570 G4int,
0571 const G4ParticleDefinition* ) const
0572 {
0573 return 0.5 * kineticEnergy;
0574 }
