273{
274
275
276
277
278
279
280 if (verboseLevel > 3)
281 G4cout <<
"Calling SampleSecondaries() of G4LivermorePolarizedGammaConversionModel" <<
G4endl;
282
284
285
286 if(photonEnergy <= lowEnergyLimit)
287 {
290 return;
291 }
292
293
296
297
298
299
300 if(!(gammaPolarization0.
isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.
mag()==0))
301 {
302 gammaPolarization0 = GetRandomPolarization(gammaDirection0);
303 }
304 else
305 {
307 {
308 gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0);
309 }
310 }
311
312
313
314
316 G4double epsilon0Local = electron_mass_c2 / photonEnergy ;
317
318
319
320 if (photonEnergy < smallEnergy )
321 {
323 }
324 else
325 {
326
327
328
331
332
333 if (element == 0)
334 {
335 G4cout <<
"G4LivermorePolarizedGammaConversionModel::SampleSecondaries - element = 0" <<
G4endl;
336 return;
337 }
338
339
341
342 if (ionisation == 0)
343 {
344 G4cout <<
"G4LivermorePolarizedGammaConversionModel::SampleSecondaries - ionisation = 0" <<
G4endl;
345 return;
346 }
347
348
349
351 if (photonEnergy > 50. * MeV) fZ += 8. * (element->
GetfCoulomb());
352
353
356 G4double screenMin = std::min(4.*screenFactor,screenMax) ;
357
358
359 G4double epsilon1 = 0.5 - 0.5 * sqrt(1. - screenMin / screenMax) ;
360 G4double epsilonMin = std::max(epsilon0Local,epsilon1);
361 G4double epsilonRange = 0.5 - epsilonMin ;
362
363
366
367 G4double f10 = ScreenFunction1(screenMin) - fZ;
368 G4double f20 = ScreenFunction2(screenMin) - fZ;
369 G4double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.);
370 G4double normF2 = std::max(1.5 * f20,0.);
371
372 do {
374 {
377 gReject = (ScreenFunction1(screen) - fZ) / f10 ;
378 }
379 else
380 {
383 gReject = (ScreenFunction2(screen) - fZ) / f20 ;
384
385
386 }
388
389 }
390
391
392
395
396
397
399 {
400 electronTotEnergy = (1. -
epsilon) * photonEnergy;
401 positronTotEnergy =
epsilon * photonEnergy;
402 }
403 else
404 {
405 positronTotEnergy = (1. -
epsilon) * photonEnergy;
406 electronTotEnergy =
epsilon * photonEnergy;
407 }
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428 G4double Ene = electronTotEnergy/electron_mass_c2;
429
432
433 SetTheta(&cosTheta,&sinTheta,Ene);
434
435
436
437
439
440
441
442
443
444 phi = SetPhi(photonEnergy);
445 psi = SetPsi(photonEnergy,phi);
446
447
448
449
450
451
452
453 Psi = psi;
454 Phi = phi;
455
456
457
462
463 if (choice2 <= 0.5)
464 {
465
466
467 }
468 else
469 {
470 phi = -phi;
471 }
472
473 if (choice <= 0.5)
474 {
475 phie = psi;
476 phip = phie+phi;
477 }
478 else
479 {
480
481
482 phip = psi;
483 phie = phip + phi;
484 }
485
486
487
488
493
494
495
496
497
498
499
500 G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
501
502 SystemOfRefChange(gammaDirection0,electronDirection,
503 gammaPolarization0);
504
506 electronDirection,
507 electronKineEnergy);
508
509
510
511 Ene = positronTotEnergy/electron_mass_c2;
512
513 cosTheta = 0.;
514 sinTheta = 0.;
515
516 SetTheta(&cosTheta,&sinTheta,Ene);
517
518
519
520 dirX = sinTheta*cos(phip);
521 dirY = sinTheta*sin(phip);
522 dirZ = cosTheta;
524
525 G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
526 SystemOfRefChange(gammaDirection0,positronDirection,
527 gammaPolarization0);
528
529
531 positronDirection, positronKineEnergy);
532
533
534 fvect->push_back(particle1);
535 fvect->push_back(particle2);
536
537
540
541}
double epsilon(double density, double temperature)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
bool isOrthogonal(const Hep3Vector &v, double epsilon=tolerance) const
double howOrthogonal(const Hep3Vector &v) const
const G4ThreeVector & GetMomentumDirection() const
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
const G4ThreeVector & GetPolarization() const
static G4Electron * Electron()
G4double GetfCoulomb() const
G4IonisParamElm * GetIonisation() const
G4double GetlogZ3() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
static G4Positron * Positron()
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
void ProposeTrackStatus(G4TrackStatus status)