292{
293
294
295
296
297
298
299
300
301
302
303
304
305 if (fVerboseLevel > 3)
306 G4cout <<
"Calling SamplingSecondaries() of G4PenelopeGammaConversionModel" <<
G4endl;
307
309
310
313
314 if (photonEnergy <= fIntrinsicLowEnergyLimit)
315 {
317 return ;
318 }
319
322
323
324
325 if (!fEffectiveCharge)
326 {
327
328
329 fLocalTable = true;
330 fEffectiveCharge = new std::map<const G4Material*,G4double>;
331 fMaterialInvScreeningRadius = new std::map<const G4Material*,G4double>;
332 fScreeningFunction = new std::map<const G4Material*,std::pair<G4double,G4double> >;
333 }
334
335 if (!fEffectiveCharge->count(mat))
336 {
337
338
339 if (fVerboseLevel > 0)
340 {
341
343 ed << "Unable to allocate the EffectiveCharge data for " <<
345 ed <<
"This can happen only in Unit Tests" <<
G4endl;
346 G4Exception(
"G4PenelopeGammaConversionModel::SampleSecondaries()",
348 }
349
350 G4AutoLock lock(&PenelopeGammaConversionModelMutex);
351 InitializeScreeningFunctions(mat);
352 lock.unlock();
353 }
354
355
357 G4double eki = electron_mass_c2/photonEnergy;
358
359
360 if (photonEnergy < fSmallEnergy)
362 else
363 {
364
365 G4double effC = fEffectiveCharge->find(mat)->second;
366 G4double alz = effC*fine_structure_const;
368 G4double F00=(-1.774-1.210e1*alz+1.118e1*alz*alz)*T
369 +(8.523+7.326e1*alz-4.441e1*alz*alz)*T*T
370 -(1.352e1+1.211e2*alz-9.641e1*alz*alz)*T*T*T
371 +(8.946+6.205e1*alz-6.341e1*alz*alz)*T*T*T*T;
372
373 G4double F0b = fScreeningFunction->find(mat)->second.second;
375 G4double invRad = fMaterialInvScreeningRadius->find(mat)->second;
377 std::pair<G4double,G4double> scree = GetScreeningFunctions(bmin);
385
387
388 do{
389 loopAgain = false;
391 {
393 if (ru2m1 < 0)
394 eps = 0.5-xr*std::pow(std::abs(ru2m1),1./3.);
395 else
396 eps = 0.5+xr*std::pow(ru2m1,1./3.);
398 scree = GetScreeningFunctions(B);
399 g1 = scree.first;
400 g1 = std::max(g1+g0,0.);
402 loopAgain = true;
403 }
404 else
405 {
408 scree = GetScreeningFunctions(B);
409 g2 = scree.second;
410 g2 = std::max(g2+g0,0.);
412 loopAgain = true;
413 }
414 }while(loopAgain);
415 }
416 if (fVerboseLevel > 4)
418
419 G4double electronTotEnergy = eps*photonEnergy;
420 G4double positronTotEnergy = (1.0-eps)*photonEnergy;
421
422
423
424
425 G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
427 G4double kk = std::sqrt(electronKineEnergy*(electronKineEnergy+2.*electron_mass_c2));
428 costheta_el = (costheta_el*electronTotEnergy+kk)/(electronTotEnergy+costheta_el*kk);
430 G4double dirX_el = std::sqrt(1.-costheta_el*costheta_el) * std::cos(phi_el);
431 G4double dirY_el = std::sqrt(1.-costheta_el*costheta_el) * std::sin(phi_el);
433
434
435 G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
437 kk = std::sqrt(positronKineEnergy*(positronKineEnergy+2.*electron_mass_c2));
438 costheta_po = (costheta_po*positronTotEnergy+kk)/(positronTotEnergy+costheta_po*kk);
440 G4double dirX_po = std::sqrt(1.-costheta_po*costheta_po) * std::cos(phi_po);
441 G4double dirY_po = std::sqrt(1.-costheta_po*costheta_po) * std::sin(phi_po);
443
444
445
446
448
449 if (electronKineEnergy > 0.0)
450 {
451 G4ThreeVector electronDirection ( dirX_el, dirY_el, dirZ_el);
452 electronDirection.rotateUz(photonDirection);
454 electronDirection,
455 electronKineEnergy);
456 fvect->push_back(electron);
457 }
458 else
459 {
460 localEnergyDeposit += electronKineEnergy;
461 electronKineEnergy = 0;
462 }
463
464
465
466 if (positronKineEnergy < 0.0)
467 {
468 localEnergyDeposit += positronKineEnergy;
469 positronKineEnergy = 0;
470 }
472 positronDirection.rotateUz(photonDirection);
474 positronDirection, positronKineEnergy);
475 fvect->push_back(positron);
476
477
479
480 if (fVerboseLevel > 1)
481 {
482 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
483 G4cout <<
"Energy balance from G4PenelopeGammaConversion" <<
G4endl;
484 G4cout <<
"Incoming photon energy: " << photonEnergy/keV <<
" keV" <<
G4endl;
485 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
486 if (electronKineEnergy)
487 G4cout <<
"Electron (explicitly produced) " << electronKineEnergy/keV <<
" keV"
489 if (positronKineEnergy)
490 G4cout <<
"Positron (not at rest) " << positronKineEnergy/keV <<
" keV" <<
G4endl;
491 G4cout <<
"Rest masses of e+/- " << 2.0*electron_mass_c2/keV <<
" keV" <<
G4endl;
492 if (localEnergyDeposit)
493 G4cout <<
"Local energy deposit " << localEnergyDeposit/keV <<
" keV" <<
G4endl;
494 G4cout <<
"Total final state: " << (electronKineEnergy+positronKineEnergy+
495 localEnergyDeposit+2.0*electron_mass_c2)/keV <<
497 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
498 }
499 if (fVerboseLevel > 0)
500 {
501 G4double energyDiff = std::fabs(electronKineEnergy+positronKineEnergy+
502 localEnergyDeposit+2.0*electron_mass_c2-photonEnergy);
503 if (energyDiff > 0.05*keV)
504 G4cout <<
"Warning from G4PenelopeGammaConversion: problem with energy conservation: "
505 << (electronKineEnergy+positronKineEnergy+
506 localEnergyDeposit+2.0*electron_mass_c2)/keV
507 <<
" keV (final) vs. " << photonEnergy/keV <<
" keV (initial)" <<
G4endl;
508 }
509}
G4double B(G4double 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)