287{
288
289
290
291
292
293
294
295
296
297
298
299
300 if (verboseLevel > 3)
301 G4cout <<
"Calling SamplingSecondaries() of G4PenelopeGammaConversionModel" <<
G4endl;
302
304
305
308
309 if (photonEnergy <= fIntrinsicLowEnergyLimit)
310 {
312 return ;
313 }
314
317
318
319
320 if (!fEffectiveCharge)
321 {
322
323
324 fLocalTable = true;
325 fEffectiveCharge = new std::map<const G4Material*,G4double>;
326 fMaterialInvScreeningRadius = new std::map<const G4Material*,G4double>;
327 fScreeningFunction = new std::map<const G4Material*,std::pair<G4double,G4double> >;
328 }
329
330 if (!fEffectiveCharge->count(mat))
331 {
332
333
334 if (verboseLevel > 0)
335 {
336
338 ed << "Unable to allocate the EffectiveCharge data for " <<
340 ed <<
"This can happen only in Unit Tests" <<
G4endl;
341 G4Exception(
"G4PenelopeGammaConversionModel::SampleSecondaries()",
343 }
344
345 G4AutoLock lock(&PenelopeGammaConversionModelMutex);
346 InitializeScreeningFunctions(mat);
347 lock.unlock();
348 }
349
350
352 G4double eki = electron_mass_c2/photonEnergy;
353
354
355 if (photonEnergy < fSmallEnergy)
357 else
358 {
359
360 G4double effC = fEffectiveCharge->find(mat)->second;
361 G4double alz = effC*fine_structure_const;
363 G4double F00=(-1.774-1.210e1*alz+1.118e1*alz*alz)*T
364 +(8.523+7.326e1*alz-4.441e1*alz*alz)*T*T
365 -(1.352e1+1.211e2*alz-9.641e1*alz*alz)*T*T*T
366 +(8.946+6.205e1*alz-6.341e1*alz*alz)*T*T*T*T;
367
368 G4double F0b = fScreeningFunction->find(mat)->second.second;
370 G4double invRad = fMaterialInvScreeningRadius->find(mat)->second;
372 std::pair<G4double,G4double> scree = GetScreeningFunctions(bmin);
380
382
383 do{
384 loopAgain = false;
386 {
388 if (ru2m1 < 0)
389 eps = 0.5-xr*std::pow(std::abs(ru2m1),1./3.);
390 else
391 eps = 0.5+xr*std::pow(ru2m1,1./3.);
393 scree = GetScreeningFunctions(B);
394 g1 = scree.first;
395 g1 = std::max(g1+g0,0.);
397 loopAgain = true;
398 }
399 else
400 {
403 scree = GetScreeningFunctions(B);
404 g2 = scree.second;
405 g2 = std::max(g2+g0,0.);
407 loopAgain = true;
408 }
409 }while(loopAgain);
410
411 }
412 if (verboseLevel > 4)
414
415 G4double electronTotEnergy = eps*photonEnergy;
416 G4double positronTotEnergy = (1.0-eps)*photonEnergy;
417
418
419
420
421 G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
423 G4double kk = std::sqrt(electronKineEnergy*(electronKineEnergy+2.*electron_mass_c2));
424 costheta_el = (costheta_el*electronTotEnergy+kk)/(electronTotEnergy+costheta_el*kk);
426 G4double dirX_el = std::sqrt(1.-costheta_el*costheta_el) * std::cos(phi_el);
427 G4double dirY_el = std::sqrt(1.-costheta_el*costheta_el) * std::sin(phi_el);
429
430
431 G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
433 kk = std::sqrt(positronKineEnergy*(positronKineEnergy+2.*electron_mass_c2));
434 costheta_po = (costheta_po*positronTotEnergy+kk)/(positronTotEnergy+costheta_po*kk);
436 G4double dirX_po = std::sqrt(1.-costheta_po*costheta_po) * std::cos(phi_po);
437 G4double dirY_po = std::sqrt(1.-costheta_po*costheta_po) * std::sin(phi_po);
439
440
441
442
444
445 if (electronKineEnergy > 0.0)
446 {
447 G4ThreeVector electronDirection ( dirX_el, dirY_el, dirZ_el);
448 electronDirection.rotateUz(photonDirection);
450 electronDirection,
451 electronKineEnergy);
452 fvect->push_back(electron);
453 }
454 else
455 {
456 localEnergyDeposit += electronKineEnergy;
457 electronKineEnergy = 0;
458 }
459
460
461
462 if (positronKineEnergy < 0.0)
463 {
464 localEnergyDeposit += positronKineEnergy;
465 positronKineEnergy = 0;
466 }
468 positronDirection.rotateUz(photonDirection);
470 positronDirection, positronKineEnergy);
471 fvect->push_back(positron);
472
473
475
476 if (verboseLevel > 1)
477 {
478 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
479 G4cout <<
"Energy balance from G4PenelopeGammaConversion" <<
G4endl;
480 G4cout <<
"Incoming photon energy: " << photonEnergy/keV <<
" keV" <<
G4endl;
481 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
482 if (electronKineEnergy)
483 G4cout <<
"Electron (explicitly produced) " << electronKineEnergy/keV <<
" keV"
485 if (positronKineEnergy)
486 G4cout <<
"Positron (not at rest) " << positronKineEnergy/keV <<
" keV" <<
G4endl;
487 G4cout <<
"Rest masses of e+/- " << 2.0*electron_mass_c2/keV <<
" keV" <<
G4endl;
488 if (localEnergyDeposit)
489 G4cout <<
"Local energy deposit " << localEnergyDeposit/keV <<
" keV" <<
G4endl;
490 G4cout <<
"Total final state: " << (electronKineEnergy+positronKineEnergy+
491 localEnergyDeposit+2.0*electron_mass_c2)/keV <<
493 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
494 }
495 if (verboseLevel > 0)
496 {
497 G4double energyDiff = std::fabs(electronKineEnergy+positronKineEnergy+
498 localEnergyDeposit+2.0*electron_mass_c2-photonEnergy);
499 if (energyDiff > 0.05*keV)
500 G4cout <<
"Warning from G4PenelopeGammaConversion: problem with energy conservation: "
501 << (electronKineEnergy+positronKineEnergy+
502 localEnergyDeposit+2.0*electron_mass_c2)/keV
503 <<
" keV (final) vs. " << photonEnergy/keV <<
" keV (initial)" <<
G4endl;
504 }
505}
double B(double temperature)
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
static G4Electron * Electron()
const G4String & GetName() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
static G4Positron * Positron()
void ProposeTrackStatus(G4TrackStatus status)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)