67 :
G4VEmModel(nam),fParticleChange(0),fParticle(0),
68 isInitialised(false),fAtomDeexcitation(0),logAtomicShellXS(0),fLocalTable(false)
99 std::map <G4int,G4PhysicsTable*>::iterator i;
119 G4cout <<
"Calling G4PenelopePhotoElectricModel::Initialise()" <<
G4endl;
126 G4cout <<
"WARNING from G4PenelopePhotoElectricModel " <<
G4endl;
127 G4cout <<
"Atomic de-excitation module is not instantiated, so there will not be ";
129 G4cout <<
"Please make sure this is intended" <<
G4endl;
153 G4int iZ = (
G4int) theElementVector->at(j)->GetZ();
164 G4cout <<
"Penelope Photo-Electric model v2008 is initialized " << G4endl
181 G4cout <<
"Calling G4PenelopePhotoElectricModel::InitialiseLocal()" <<
G4endl;
217 G4cout <<
"Calling ComputeCrossSectionPerAtom() of G4PenelopePhotoElectricModel" <<
G4endl;
240 ed <<
"Unable to retrieve the shell cross section table for Z=" << iZ <<
G4endl;
241 ed <<
"This can happen only in Unit Tests or via G4EmCalculator" <<
G4endl;
242 G4Exception(
"G4PenelopePhotoElectricModel::ComputeCrossSectionPerAtom()",
246 G4AutoLock lock(&PenelopePhotoElectricModelMutex);
259 G4Exception(
"G4PenelopePhotoElectricModel::ComputeCrossSectionPerAtom()",
261 "Unable to retrieve the total cross section table");
266 cross = std::exp(logXS);
269 G4cout <<
"Photoelectric cross section at " << energy/
MeV <<
" MeV for Z=" << Z <<
297 G4cout <<
"Calling SamplingSecondaries() of G4PenelopePhotoElectricModel" <<
G4endl;
333 G4cout <<
"Selected shell " << shellIndex <<
" of element " << anElement->
GetName() <<
G4endl;
344 size_t numberOfShells = (size_t) transitionManager->
NumberOfShells(Z);
345 if (shellIndex >= numberOfShells)
346 shellIndex = numberOfShells-1;
358 bindingEnergy = 0.*
eV;
369 if (eKineticEnergy > 0.)
374 G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
376 G4double dirx = sinTheta * std::cos(phi);
377 G4double diry = sinTheta * std::sin(phi);
380 electronDirection.rotateUz(photonDirection);
384 fvect->push_back(electron);
387 bindingEnergy = photonEnergy;
401 size_t nBefore = fvect->size();
403 size_t nAfter = fvect->size();
405 if (nAfter > nBefore)
407 for (
size_t j=nBefore;j<nAfter;j++)
409 G4double itsEnergy = ((*fvect)[j])->GetKineticEnergy();
410 bindingEnergy -= itsEnergy;
412 energyInFluorescence += itsEnergy;
414 energyInAuger += itsEnergy;
423 if (localEnergyDeposit < 0)
426 <<
"G4PenelopePhotoElectricModel::SampleSecondaries() - Negative energy deposit"
428 localEnergyDeposit = 0;
435 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
436 G4cout <<
"Energy balance from G4PenelopePhotoElectric" <<
G4endl;
439 G4cout <<
"Incoming photon energy: " << photonEnergy/
keV <<
" keV" <<
G4endl;
440 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
442 G4cout <<
"Outgoing electron " << eKineticEnergy/
keV <<
" keV" <<
G4endl;
443 if (energyInFluorescence)
444 G4cout <<
"Fluorescence x-rays: " << energyInFluorescence/
keV <<
" keV" <<
G4endl;
446 G4cout <<
"Auger electrons: " << energyInAuger/
keV <<
" keV" <<
G4endl;
447 G4cout <<
"Local energy deposit " << localEnergyDeposit/
keV <<
" keV" <<
G4endl;
448 G4cout <<
"Total final state: " <<
449 (eKineticEnergy+energyInFluorescence+localEnergyDeposit+energyInAuger)/
keV <<
451 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
456 std::fabs(eKineticEnergy+energyInFluorescence+localEnergyDeposit+energyInAuger-photonEnergy);
457 if (energyDiff > 0.05*
keV)
459 G4cout <<
"Warning from G4PenelopePhotoElectric: problem with energy conservation: " <<
460 (eKineticEnergy+energyInFluorescence+localEnergyDeposit+energyInAuger)/
keV
461 <<
" keV (final) vs. " <<
462 photonEnergy/
keV <<
" keV (initial)" << G4endl;
463 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
464 G4cout <<
"Energy balance from G4PenelopePhotoElectric" <<
G4endl;
467 G4cout <<
"Incoming photon energy: " << photonEnergy/
keV <<
" keV" <<
G4endl;
468 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
470 G4cout <<
"Outgoing electron " << eKineticEnergy/
keV <<
" keV" <<
G4endl;
471 if (energyInFluorescence)
472 G4cout <<
"Fluorescence x-rays: " << energyInFluorescence/
keV <<
" keV" <<
G4endl;
474 G4cout <<
"Auger electrons: " << energyInAuger/
keV <<
" keV" <<
G4endl;
475 G4cout <<
"Local energy deposit " << localEnergyDeposit/
keV <<
" keV" <<
G4endl;
476 G4cout <<
"Total final state: " <<
477 (eKineticEnergy+energyInFluorescence+localEnergyDeposit+energyInAuger)/
keV <<
479 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
489 if (energy>1*
GeV)
return costheta;
494 G4double gamma = 1.0 + energy/electron_mass_c2;
496 G4double beta = std::sqrt((gamma2-1.0)/gamma2);
501 G4double a1 = 0.5*beta*gamma*(gamma-1.0)*(gamma-2.0);
513 tsam = 2.0*ac * (2.0*rand + a2*std::sqrt(rand)) / (a2*a2 - 4.0*rand);
514 gtr = (2.0 - tsam) * (a1 + 1.0/(ac+tsam));
528 G4Exception(
"G4PenelopePhotoElectricModel::ReadDataFile()",
533 G4cout <<
"G4PenelopePhotoElectricModel::ReadDataFile()" <<
G4endl;
534 G4cout <<
"Going to read PhotoElectric data files for Z=" << Z <<
G4endl;
537 char* path = getenv(
"G4LEDATA");
540 G4String excep =
"G4PenelopePhotoElectricModel - G4LEDATA environment variable not set!";
541 G4Exception(
"G4PenelopePhotoElectricModel::ReadDataFile()",
549 std::ostringstream ost;
551 ost << path <<
"/penelope/photoelectric/pdgph" << Z <<
".p08";
553 ost << path <<
"/penelope/photoelectric/pdgph0" << Z <<
".p08";
554 std::ifstream file(ost.str().c_str());
557 G4String excep =
"G4PenelopePhotoElectricModel - data file " +
G4String(ost.str()) +
" not found!";
558 G4Exception(
"G4PenelopePhotoElectricModel::ReadDataFile()",
565 while( getline(file, line) )
572 file.open(ost.str().c_str());
576 file >> readZ >> nShells;
579 G4cout <<
"Element Z=" << Z <<
" , nShells = " << nShells <<
G4endl;
582 if (readZ != Z || nShells <= 0 || nShells > 50)
585 ed <<
"Corrupted data file for Z=" << Z <<
G4endl;
586 G4Exception(
"G4PenelopePhotoElectricModel::ReadDataFile()",
598 for (
size_t i=0;i<nShells+1;i++)
602 for (k=0;k<ndata && !file.eof();k++)
610 for (
size_t i=0;i<nShells+1;i++)
615 if (aValue < 1e-40*
cm2)
617 theVec->
PutValue(k,logene,std::log(aValue));
623 G4cout <<
"G4PenelopePhotoElectricModel: read " << k <<
" points for element Z = "
639 G4Exception(
"G4PenelopePhotoElectricModel::GetNumberOfShellXS()",
649 ed <<
"Cannot find shell cross section data for Z=" << Z <<
G4endl;
650 G4Exception(
"G4PenelopePhotoElectricModel::GetNumberOfShellXS()",
665 if (shellID >= entries)
667 G4cout <<
"Element Z=" << Z <<
" has data for " << entries <<
" shells only" <<
G4endl;
668 G4cout <<
"so shellID should be from 0 to " << entries-1 <<
G4endl;
678 G4Exception(
"G4PenelopePhotoElectricModel::GetShellCrossSection()",
680 "Unable to retrieve the total cross section table");
686 if (cross < 2e-40*
cm2) cross = 0;
697 else if (shellID == 1)
699 else if (shellID == 2)
701 else if (shellID == 3)
703 else if (shellID == 4)
705 else if (shellID == 5)
707 else if (shellID == 6)
709 else if (shellID == 7)
711 else if (shellID == 8)
730 G4double logEnergy = std::log(energy);
736 ed <<
"Cannot find shell cross section data for Z=" << Z <<
G4endl;
737 G4Exception(
"G4PenelopePhotoElectricModel::SelectRandomShell()",
759 for (
size_t k=1;k<theTable->
entries();k++)
764 G4double partialXS = std::exp(logXSLocal);
G4double LowEnergyLimit() const
G4bool CheckDeexcitationActiveRegion(G4int coupleIndex)
static G4LossTableManager * Instance()
G4int NumberOfShells(G4int Z) const
std::vector< G4Element * > G4ElementVector
void PutValue(size_t binNumber, G4double binValue, G4double dataValue)
std::ostringstream G4ExceptionDescription
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
G4double HighEnergyLimit() const
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
G4ParticleChangeForGamma * fParticleChange
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel)
const G4String & GetName() const
G4double GetShellCrossSection(G4int Z, size_t shellID, G4double energy)
void push_back(G4PhysicsVector *)
static G4Electron * Definition()
G4double BindingEnergy() const
G4PenelopePhotoElectricModel(const G4ParticleDefinition *p=0, const G4String &processName="PenPhotoElec")
const G4ElementVector * GetElementVector() const
std::map< G4int, G4PhysicsTable * > * logAtomicShellXS
const G4AtomicTransitionManager * fTransitionManager
#define G4MUTEX_INITIALIZER
G4double fIntrinsicLowEnergyLimit
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void SetHighEnergyLimit(G4double)
G4GLOB_DLL std::ostream G4cout
size_t GetTableSize() const
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
std::ostream & tab(std::ostream &)
const G4ThreeVector & GetMomentumDirection() const
static const double twopi
void ReadDataFile(G4int Z)
std::vector< G4EmElementSelector * > * GetElementSelectors()
G4double Value(G4double theEnergy, size_t &lastidx) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
static G4ProductionCutsTable * GetProductionCutsTable()
void SetElementSelectors(std::vector< G4EmElementSelector * > *)
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
virtual ~G4PenelopePhotoElectricModel()
G4double energy(const ThreeVector &p, const G4double m)
G4VAtomDeexcitation * fAtomDeexcitation
void SetParticle(const G4ParticleDefinition *)
const G4ParticleDefinition * fParticle
size_t SelectRandomShell(G4int Z, G4double energy)
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
static G4Electron * Electron()
void SetProposedKineticEnergy(G4double proposedKinEnergy)
size_t GetNumberOfElements() const
G4VAtomDeexcitation * AtomDeexcitation()
void ProposeTrackStatus(G4TrackStatus status)
const G4String & GetName() const
size_t GetNumberOfShellXS(G4int)
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
void SetDeexcitationFlag(G4bool val)
G4double fIntrinsicHighEnergyLimit
G4ThreeVector G4ParticleMomentum
G4double bindingEnergy(G4int A, G4int Z)
static G4AtomicTransitionManager * Instance()
G4double SampleElectronDirection(G4double energy)
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
G4ParticleChangeForGamma * GetParticleChangeForGamma()
static G4Gamma * Definition()
const G4Material * GetMaterial() const
static G4Gamma * GammaDefinition()
G4String WriteTargetShell(size_t shellID)