67 :
G4VEmModel(nam),fParticleChange(0),fParticle(0),
68 isInitialised(
false),fAtomDeexcitation(0),logAtomicShellXS(0),fLocalTable(
false)
118 G4cout <<
"Calling G4PenelopePhotoElectricModel::Initialise()" <<
G4endl;
125 G4cout <<
"WARNING from G4PenelopePhotoElectricModel " <<
G4endl;
126 G4cout <<
"Atomic de-excitation module is not instantiated, so there will not be ";
128 G4cout <<
"Please make sure this is intended" <<
G4endl;
152 G4int iZ = (
G4int) theElementVector->at(j)->GetZ();
163 G4cout <<
"Penelope Photo-Electric model v2008 is initialized " << G4endl
180 G4cout <<
"Calling G4PenelopePhotoElectricModel::InitialiseLocal()" <<
G4endl;
216 G4cout <<
"Calling ComputeCrossSectionPerAtom() of G4PenelopePhotoElectricModel" <<
G4endl;
239 ed <<
"Unable to retrieve the shell cross section table for Z=" << iZ <<
G4endl;
240 ed <<
"This can happen only in Unit Tests or via G4EmCalculator" <<
G4endl;
241 G4Exception(
"G4PenelopePhotoElectricModel::ComputeCrossSectionPerAtom()",
245 G4AutoLock lock(&PenelopePhotoElectricModelMutex);
258 G4Exception(
"G4PenelopePhotoElectricModel::ComputeCrossSectionPerAtom()",
260 "Unable to retrieve the total cross section table");
265 cross =
G4Exp(logXS);
268 G4cout <<
"Photoelectric cross section at " << energy/
MeV <<
" MeV for Z=" << Z <<
296 G4cout <<
"Calling SamplingSecondaries() of G4PenelopePhotoElectricModel" <<
G4endl;
332 G4cout <<
"Selected shell " << shellIndex <<
" of element " << anElement->
GetName() <<
G4endl;
343 size_t numberOfShells = (size_t) transitionManager->
NumberOfShells(Z);
344 if (shellIndex >= numberOfShells)
345 shellIndex = numberOfShells-1;
357 bindingEnergy = 0.*
eV;
368 if (eKineticEnergy > 0.)
373 G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
375 G4double dirx = sinTheta * std::cos(phi);
376 G4double diry = sinTheta * std::sin(phi);
379 electronDirection.
rotateUz(photonDirection);
383 fvect->push_back(electron);
386 bindingEnergy = photonEnergy;
400 size_t nBefore = fvect->size();
402 size_t nAfter = fvect->size();
404 if (nAfter > nBefore)
406 for (
size_t j=nBefore;j<nAfter;j++)
408 G4double itsEnergy = ((*fvect)[j])->GetKineticEnergy();
409 if (itsEnergy < bindingEnergy)
411 bindingEnergy -= itsEnergy;
413 energyInFluorescence += itsEnergy;
415 energyInAuger += itsEnergy;
420 (*fvect)[j] =
nullptr;
430 if (localEnergyDeposit < 0)
432 G4Exception(
"G4PenelopePhotoElectricModel::SampleSecondaries()",
433 "em2099",
JustWarning,
"WARNING: Negative local energy deposit");
434 localEnergyDeposit = 0;
441 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
442 G4cout <<
"Energy balance from G4PenelopePhotoElectric" <<
G4endl;
445 G4cout <<
"Incoming photon energy: " << photonEnergy/
keV <<
" keV" <<
G4endl;
446 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
448 G4cout <<
"Outgoing electron " << eKineticEnergy/
keV <<
" keV" <<
G4endl;
449 if (energyInFluorescence)
450 G4cout <<
"Fluorescence x-rays: " << energyInFluorescence/
keV <<
" keV" <<
G4endl;
452 G4cout <<
"Auger electrons: " << energyInAuger/
keV <<
" keV" <<
G4endl;
453 G4cout <<
"Local energy deposit " << localEnergyDeposit/
keV <<
" keV" <<
G4endl;
454 G4cout <<
"Total final state: " <<
455 (eKineticEnergy+energyInFluorescence+localEnergyDeposit+energyInAuger)/
keV <<
457 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
462 std::fabs(eKineticEnergy+energyInFluorescence+localEnergyDeposit+energyInAuger-photonEnergy);
463 if (energyDiff > 0.05*
keV)
465 G4cout <<
"Warning from G4PenelopePhotoElectric: problem with energy conservation: " <<
466 (eKineticEnergy+energyInFluorescence+localEnergyDeposit+energyInAuger)/
keV
467 <<
" keV (final) vs. " <<
468 photonEnergy/
keV <<
" keV (initial)" << G4endl;
469 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
470 G4cout <<
"Energy balance from G4PenelopePhotoElectric" <<
G4endl;
473 G4cout <<
"Incoming photon energy: " << photonEnergy/
keV <<
" keV" <<
G4endl;
474 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
476 G4cout <<
"Outgoing electron " << eKineticEnergy/
keV <<
" keV" <<
G4endl;
477 if (energyInFluorescence)
478 G4cout <<
"Fluorescence x-rays: " << energyInFluorescence/
keV <<
" keV" <<
G4endl;
480 G4cout <<
"Auger electrons: " << energyInAuger/
keV <<
" keV" <<
G4endl;
481 G4cout <<
"Local energy deposit " << localEnergyDeposit/
keV <<
" keV" <<
G4endl;
482 G4cout <<
"Total final state: " <<
483 (eKineticEnergy+energyInFluorescence+localEnergyDeposit+energyInAuger)/
keV <<
485 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
495 if (energy>1*
GeV)
return costheta;
502 G4double beta = std::sqrt((gamma2-1.0)/gamma2);
507 G4double a1 = 0.5*beta*gamma*(gamma-1.0)*(gamma-2.0);
519 tsam = 2.0*ac * (2.0*rand + a2*std::sqrt(rand)) / (a2*a2 - 4.0*rand);
520 gtr = (2.0 - tsam) * (a1 + 1.0/(ac+tsam));
534 G4Exception(
"G4PenelopePhotoElectricModel::ReadDataFile()",
539 G4cout <<
"G4PenelopePhotoElectricModel::ReadDataFile()" <<
G4endl;
540 G4cout <<
"Going to read PhotoElectric data files for Z=" << Z <<
G4endl;
543 char* path = std::getenv(
"G4LEDATA");
546 G4String excep =
"G4PenelopePhotoElectricModel - G4LEDATA environment variable not set!";
547 G4Exception(
"G4PenelopePhotoElectricModel::ReadDataFile()",
555 std::ostringstream ost;
557 ost << path <<
"/penelope/photoelectric/pdgph" << Z <<
".p08";
559 ost << path <<
"/penelope/photoelectric/pdgph0" << Z <<
".p08";
560 std::ifstream
file(ost.str().c_str());
563 G4String excep =
"G4PenelopePhotoElectricModel - data file " +
G4String(ost.str()) +
" not found!";
564 G4Exception(
"G4PenelopePhotoElectricModel::ReadDataFile()",
571 while( getline(
file, line) )
578 file.open(ost.str().c_str());
582 file >> readZ >> nShells;
585 G4cout <<
"Element Z=" << Z <<
" , nShells = " << nShells <<
G4endl;
588 if (readZ != Z || nShells <= 0 || nShells > 50)
591 ed <<
"Corrupted data file for Z=" << Z <<
G4endl;
592 G4Exception(
"G4PenelopePhotoElectricModel::ReadDataFile()",
604 for (
size_t i=0;i<nShells+1;i++)
608 for (k=0;k<ndata && !
file.eof();k++)
616 for (
size_t i=0;i<nShells+1;i++)
621 if (aValue < 1
e-40*
cm2)
623 theVec->
PutValue(k,logene,std::log(aValue));
629 G4cout <<
"G4PenelopePhotoElectricModel: read " << k <<
" points for element Z = "
645 G4Exception(
"G4PenelopePhotoElectricModel::GetNumberOfShellXS()",
655 ed <<
"Cannot find shell cross section data for Z=" << Z <<
G4endl;
656 G4Exception(
"G4PenelopePhotoElectricModel::GetNumberOfShellXS()",
671 if (shellID >= entries)
673 G4cout <<
"Element Z=" << Z <<
" has data for " << entries <<
" shells only" <<
G4endl;
674 G4cout <<
"so shellID should be from 0 to " << entries-1 <<
G4endl;
684 G4Exception(
"G4PenelopePhotoElectricModel::GetShellCrossSection()",
686 "Unable to retrieve the total cross section table");
692 if (cross < 2
e-40*
cm2) cross = 0;
703 else if (shellID == 1)
705 else if (shellID == 2)
707 else if (shellID == 3)
709 else if (shellID == 4)
711 else if (shellID == 5)
713 else if (shellID == 6)
715 else if (shellID == 7)
717 else if (shellID == 8)
736 G4double logEnergy = std::log(energy);
742 ed <<
"Cannot find shell cross section data for Z=" << Z <<
G4endl;
743 G4Exception(
"G4PenelopePhotoElectricModel::SelectRandomShell()",