192{
193
195
197
201
203 {
205 }
206 else
207 {
209 {
211 }
212 else
213 {
214 G4cout <<
"*** WARNING in G4ecpssrBaseKxsModel::CalculateCrossSection : we can treat only Proton or Alpha incident particles " <<
G4endl;
215 return 0;
216 }
217 }
218
220
222
224
226
228
230
232
234
235
236 G4double screenedzTarget = zTarget-zkshell;
237
238
239 constexpr G4double rydbergMeV= 13.6056923e-6;
240
241 G4double tetaK = kBindingEnergy/((screenedzTarget*screenedzTarget)*rydbergMeV);
242
243
245
246 G4double velocity =(2./(tetaK*screenedzTarget))*std::pow(((energyIncident*
electron_mass_c2)/(massIncident*rydbergMeV)),0.5);
247
248
249
250
252
254
256
257 G4double sigma0 = 8.*
pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.);
258
259
260
262
263 const G4double kAnalyticalApproximation= 1.5;
264 G4double x = kAnalyticalApproximation/velocity;
265
266
267
269
271
272
273
274
275 if ((0.< x) && (x <= 0.035))
276 {
277 electrIonizationEnergy= 0.75*
pi*(std::log(1./(x*x))-1.);
278 }
279 else
280 {
281 if ( (0.035 < x) && (x <=3.))
282 {
283 electrIonizationEnergy =
G4Exp(-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));
284 }
285
286 else
287 {
288 if ( (3.< x) && (x<=11.))
289 {
290 electrIonizationEnergy =2.*
G4Exp(-2.*x)/std::pow(x,1.6);
291 }
292
293 else electrIonizationEnergy =0.;
294 }
295 }
296
298
299 G4double hFunction =(electrIonizationEnergy*2.)/(tetaK*std::pow(velocity,3));
300
301
303
304 G4double gFunction = (1.+(9.*velocity)+(31.*velocity*velocity)+(98.*std::pow(velocity,3.))+(12.*std::pow(velocity,4.))+(25.*std::pow(velocity,5.))
305 +(4.2*std::pow(velocity,6.))+(0.515*std::pow(velocity,7.)))/std::pow(1.+velocity,9.);
306
307
309
310
311
312 G4double sigmaPSS = 1.+(((2.*zIncident)/(screenedzTarget*tetaK))*(gFunction-hFunction));
313
314
315
316
318
320
321
322
324
326
327 G4double ykFormula=0.4*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(velocity/sigmaPSS);
328
329
330
331
333
334 G4double relativityCorrection = std::pow((1.+(1.1*ykFormula*ykFormula)),0.5)+ykFormula;
335
336
337
339
340 G4double reducedVelocity = velocity*std::pow(relativityCorrection,0.5);
341
342
343
344
346
347 G4double etaOverTheta2 = (energyIncident*
electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget)
348 /(sigmaPSS*tetaK)/(sigmaPSS*tetaK);
349
350
351
353
355
356
357
358 if ( velocity < 1. )
359
360
361
362
363
364 {
366
367 universalFunction = (std::pow(2.,9.)/45.)*std::pow(reducedVelocity/sigmaPSS,8.)*std::pow((1.+(1.72*(reducedVelocity/sigmaPSS)*(reducedVelocity/sigmaPSS))),-4.);
368
369
371
372 }
373 else
374 {
375
376 if ( etaOverTheta2 > 86.6 && (sigmaPSS*tetaK) > 0.4 && (sigmaPSS*tetaK) < 2.9996 )
377 {
378
379
381
383
387
391
393
394
396
397 G4double etaT = (sigmaPSS*tetaK)*(sigmaPSS*tetaK)*(86.6);
398
399
401
403
404
406 {
408 "*** WARNING in G4ecpssrBaseKxsModel::CalculateCrossSection : unable to interpolate FK function in high velocity region ! ***" <<
G4endl;
409 return 0;
410 }
411
413
415
417
419
421
423
425
427
429
431
432 G4double universalFunction3= fKK/(etaK/tetaK);
433
434
436
437 universalFunction=universalFunction3;
438
439 }
440 else if ( etaOverTheta2 >= 1.e-3 && etaOverTheta2 <= 86.6 && (sigmaPSS*tetaK) >= 0.4 && (sigmaPSS*tetaK) <= 2.9996 )
441 {
442
443
445
447
448 if (universalFunction2<=0)
449 {
451 "*** WARNING : G4ecpssrBaseKxsModel::CalculateCrossSection is unable to interpolate FK function in medium velocity region ! ***" <<
G4endl;
452 return 0;
453 }
454
455 if (
verboseLevel>0)
G4cout <<
" universalFunction2=" << universalFunction2 <<
" for theta=" << sigmaPSS*tetaK <<
" and etaOverTheta2=" << etaOverTheta2 <<
G4endl;
456
457 universalFunction=universalFunction2;
458 }
459
460 }
461
462
463
464 G4double sigmaPSSR = (sigma0/(sigmaPSS*tetaK))*universalFunction;
465
466
468
469
470
471 G4double pssDeltaK = (4./(systemMass*sigmaPSS*tetaK))*(sigmaPSS/velocity)*(sigmaPSS/velocity);
472
473
474
476
477 if (pssDeltaK>1) return 0.;
478
479 G4double energyLoss = std::pow(1-pssDeltaK,0.5);
480
481
482
484
485 G4double energyLossFunction = (std::pow(2.,-9)/8.)*((((9.*energyLoss)-1.)*std::pow(1.+energyLoss,9.))+(((9.*energyLoss)+1.)*std::pow(1.-energyLoss,9.)));
486
487
488
490
491
492
493 G4double coulombDeflection = (4.*
pi*zIncident/systemMass)*std::pow(tetaK*sigmaPSS,-2.)*std::pow(velocity/sigmaPSS,-3.)*(zTarget/screenedzTarget);
494
495
497
498 G4double cParameter = 2.*coulombDeflection/(energyLoss*(energyLoss+1.));
499
500
502
504
505
507
509
510
511
513
514 crossSection = energyLossFunction* coulombDeflectionFunction*sigmaPSSR;
515
516
517
518
519
520 if (crossSection >= 0) {
521 return crossSection *
barn;
522 }
523 else {return 0;}
524
525}
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
static constexpr double barn
static constexpr double eplus
static constexpr double eV
static constexpr double pi
G4GLOB_DLL std::ostream G4cout
G4double BindingEnergy() const
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const
static G4AtomicTransitionManager * Instance()
virtual G4double FindValue(G4double e, G4int componentId=0) const
static G4NistManager * Instance()
G4double GetAtomicMassAmu(const G4String &symb) const
G4double GetPDGMass() const
G4double GetPDGCharge() const
static G4Proton * Proton()
G4double FunctionFK(G4double k, G4double theta)
G4double ExpIntFunction(G4int n, G4double x)