G4PenelopePhotoElectricModel Class Reference

#include <G4PenelopePhotoElectricModel.hh>

Inheritance diagram for G4PenelopePhotoElectricModel:

Public Member Functions
	G4PenelopePhotoElectricModel (const G4ParticleDefinition *p=0, const G4String &processName="PenPhotoElec")
virtual	~G4PenelopePhotoElectricModel ()
virtual void	Initialise (const G4ParticleDefinition *, const G4DataVector &)
virtual G4double	ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
virtual void	SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
void	SetVerbosityLevel (G4int lev)
G4int	GetVerbosityLevel ()
size_t	GetNumberOfShellXS (G4int)
G4double	GetShellCrossSection (G4int Z, size_t shellID, G4double energy)
Protected Attributes
G4ParticleChangeForGamma *	fParticleChange

---

Detailed Description

Definition at line 58 of file G4PenelopePhotoElectricModel.hh.

---

Constructor & Destructor Documentation

G4PenelopePhotoElectricModel::G4PenelopePhotoElectricModel	(	const G4ParticleDefinition *	p = 0,
		const G4String &	processName = "PenPhotoElec"
	)

Definition at line 60 of file G4PenelopePhotoElectricModel.cc.

References G4AtomicTransitionManager::Instance(), G4VEmModel::SetDeexcitationFlag(), and G4VEmModel::SetHighEnergyLimit().

  :G4VEmModel(nam),fParticleChange(0),isInitialised(false),fAtomDeexcitation(0),
   logAtomicShellXS(0)
{
  fIntrinsicLowEnergyLimit = 100.0*eV;
  fIntrinsicHighEnergyLimit = 100.0*GeV;
  //  SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
  SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
  //
  verboseLevel= 0;
  // Verbosity scale:
  // 0 = nothing 
  // 1 = warning for energy non-conservation 
  // 2 = details of energy budget
  // 3 = calculation of cross sections, file openings, sampling of atoms
  // 4 = entering in methods

  //Mark this model as "applicable" for atomic deexcitation
  SetDeexcitationFlag(true);

  fTransitionManager = G4AtomicTransitionManager::Instance();
}

G4PenelopePhotoElectricModel::~G4PenelopePhotoElectricModel

(

)

[virtual]

Definition at line 86 of file G4PenelopePhotoElectricModel.cc.

{  
  std::map <G4int,G4PhysicsTable*>::iterator i;
  if (logAtomicShellXS)
    {
      for (i=logAtomicShellXS->begin();i != logAtomicShellXS->end();i++)
        {
          G4PhysicsTable* tab = i->second;
          tab->clearAndDestroy();
          delete tab;
        }
    }
  delete logAtomicShellXS;
}

---

Member Function Documentation

G4double G4PenelopePhotoElectricModel::ComputeCrossSectionPerAtom	(	const G4ParticleDefinition *	,
		G4double	kinEnergy,
		G4double	Z,
		G4double	A = 0,
		G4double	cut = 0,
		G4double	emax = DBL_MAX
	)			[virtual]

Reimplemented from G4VEmModel.

Definition at line 141 of file G4PenelopePhotoElectricModel.cc.

References FatalException, G4cout, G4endl, G4Exception(), and G4PhysicsVector::Value().

{
  //
  // Penelope model v2008
  //

  if (verboseLevel > 3)
    G4cout << "Calling ComputeCrossSectionPerAtom() of G4PenelopePhotoElectricModel" << G4endl;

  G4int iZ = (G4int) Z;

  //read data files
  if (!logAtomicShellXS->count(iZ))
    ReadDataFile(iZ);
  //now it should be ok
  if (!logAtomicShellXS->count(iZ))     
    G4Exception("G4PenelopePhotoElectricModel::ComputeCrossSectionPerAtom()",
                "em2038",FatalException,
                "Unable to retrieve the shell cross section table");     
  
  G4double cross = 0;

  G4PhysicsTable* theTable =  logAtomicShellXS->find(iZ)->second;
  G4PhysicsFreeVector* totalXSLog = (G4PhysicsFreeVector*) (*theTable)[0];

   if (!totalXSLog)
     {
       G4Exception("G4PenelopePhotoElectricModel::ComputeCrossSectionPerAtom()",
                   "em2039",FatalException,
                   "Unable to retrieve the total cross section table");       
       return 0;
     }
   G4double logene = std::log(energy);
   G4double logXS = totalXSLog->Value(logene);
   cross = std::exp(logXS);
 
  if (verboseLevel > 2)
    G4cout << "Photoelectric cross section at " << energy/MeV << " MeV for Z=" << Z <<
      " = " << cross/barn << " barn" << G4endl;
  return cross;
}

size_t G4PenelopePhotoElectricModel::GetNumberOfShellXS

(

G4int

)

Definition at line 620 of file G4PenelopePhotoElectricModel.cc.

References FatalException, and G4Exception().

Referenced by GetShellCrossSection().

{
  //read data files
  if (!logAtomicShellXS->count(Z))
    ReadDataFile(Z);
  //now it should be ok
  if (!logAtomicShellXS->count(Z))
     {
       G4ExceptionDescription ed;
       ed << "Cannot find shell cross section data for Z=" << Z << G4endl;
       G4Exception("G4PenelopePhotoElectricModel::GetNumberOfShellXS()",
                   "em2038",FatalException,ed);
     }
  //one vector is allocated for the _total_ cross section
  size_t nEntries = logAtomicShellXS->find(Z)->second->entries();
  return  (nEntries-1);
}

G4double G4PenelopePhotoElectricModel::GetShellCrossSection	(	G4int	Z,
		size_t	shellID,
		G4double	energy
	)

Definition at line 640 of file G4PenelopePhotoElectricModel.cc.

References FatalException, G4cout, G4Exception(), GetNumberOfShellXS(), and G4PhysicsVector::Value().

{
  //this forces also the loading of the data
  size_t entries = GetNumberOfShellXS(Z);

  if (shellID >= entries)
    {
      G4cout << "Element Z=" << Z << " has data for " << entries << " shells only" << G4endl;
      G4cout << "so shellID should be from 0 to " << entries-1 << G4endl;
      return 0;
    }
  
  G4PhysicsTable* theTable =  logAtomicShellXS->find(Z)->second;
  //[0] is the total XS, shellID is in the element [shellID+1]
  G4PhysicsFreeVector* totalXSLog = (G4PhysicsFreeVector*) (*theTable)[shellID+1];
 
  if (!totalXSLog)
     {
       G4Exception("G4PenelopePhotoElectricModel::GetShellCrossSection()",
                   "em2039",FatalException,
                   "Unable to retrieve the total cross section table");
       return 0;
     }
   G4double logene = std::log(energy);
   G4double logXS = totalXSLog->Value(logene);
   G4double cross = std::exp(logXS);
   if (cross < 2e-40*cm2) cross = 0;
   return cross;
}

G4int G4PenelopePhotoElectricModel::GetVerbosityLevel

(

)

[inline]

Definition at line 85 of file G4PenelopePhotoElectricModel.hh.

{return verboseLevel;};

void G4PenelopePhotoElectricModel::Initialise	(	const G4ParticleDefinition *	,
		const G4DataVector &
	)			[virtual]

Implements G4VEmModel.

Definition at line 103 of file G4PenelopePhotoElectricModel.cc.

References G4LossTableManager::AtomDeexcitation(), fParticleChange, G4cout, G4endl, G4VEmModel::GetParticleChangeForGamma(), G4VEmModel::HighEnergyLimit(), G4VEmModel::InitialiseElementSelectors(), G4LossTableManager::Instance(), and G4VEmModel::LowEnergyLimit().

{
  if (verboseLevel > 3)
    G4cout << "Calling  G4PenelopePhotoElectricModel::Initialise()" << G4endl;

  // logAtomicShellXS is created only once, since it is  never cleared
  if (!logAtomicShellXS)
    logAtomicShellXS = new std::map<G4int,G4PhysicsTable*>;

  InitialiseElementSelectors(particle,cuts);

  fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
  //Issue warning if the AtomicDeexcitation has not been declared
  if (!fAtomDeexcitation)
    {
      G4cout << G4endl;
      G4cout << "WARNING from G4PenelopePhotoElectricModel " << G4endl;
      G4cout << "Atomic de-excitation module is not instantiated, so there will not be ";
      G4cout << "any fluorescence/Auger emission." << G4endl;
      G4cout << "Please make sure this is intended" << G4endl;
    }

  if (verboseLevel > 0) { 
    G4cout << "Penelope Photo-Electric model v2008 is initialized " << G4endl
           << "Energy range: "
           << LowEnergyLimit() / MeV << " MeV - "
           << HighEnergyLimit() / GeV << " GeV";
  }

  if(isInitialised) return;
  fParticleChange = GetParticleChangeForGamma();
  isInitialised = true;

}

void G4PenelopePhotoElectricModel::SampleSecondaries	(	std::vector< G4DynamicParticle * > *	,
		const G4MaterialCutsCouple *	,
		const G4DynamicParticle *	,
		G4double	tmin,
		G4double	maxEnergy
	)			[virtual]

Implements G4VEmModel.

Definition at line 189 of file G4PenelopePhotoElectricModel.cc.

References G4AtomicShell::BindingEnergy(), G4InuclSpecialFunctions::bindingEnergy(), G4VAtomDeexcitation::CheckDeexcitationActiveRegion(), G4Electron::Definition(), G4Gamma::Definition(), G4Electron::Electron(), fParticleChange, fStopAndKill, G4cout, G4endl, G4UniformRand, G4Gamma::GammaDefinition(), G4VAtomDeexcitation::GenerateParticles(), G4MaterialCutsCouple::GetIndex(), G4DynamicParticle::GetKineticEnergy(), G4MaterialCutsCouple::GetMaterial(), G4DynamicParticle::GetMomentumDirection(), G4Element::GetName(), G4Material::GetName(), G4Element::GetZ(), G4AtomicTransitionManager::Instance(), G4AtomicTransitionManager::NumberOfShells(), G4VParticleChange::ProposeLocalEnergyDeposit(), G4VParticleChange::ProposeTrackStatus(), G4VEmModel::SelectRandomAtom(), G4ParticleChangeForGamma::SetProposedKineticEnergy(), and G4AtomicTransitionManager::Shell().

{
  //
  // Photoelectric effect, Penelope model v2008
  //
  // The target atom and the target shell are sampled according to the Livermore 
  // database 
  //  D.E. Cullen et al., Report UCRL-50400 (1989)
  // The angular distribution of the electron in the final state is sampled 
  // according to the Sauter distribution from 
  //  F. Sauter, Ann. Phys. 11 (1931) 454
  // The energy of the final electron is given by the initial photon energy minus 
  // the binding energy. Fluorescence de-excitation is subsequently produced 
  // (to fill the vacancy) according to the general Geant4 G4DeexcitationManager:
  //  J. Stepanek, Comp. Phys. Comm. 1206 pp 1-1-9 (1997)

  if (verboseLevel > 3)
    G4cout << "Calling SamplingSecondaries() of G4PenelopePhotoElectricModel" << G4endl;

  G4double photonEnergy = aDynamicGamma->GetKineticEnergy();

  // always kill primary
  fParticleChange->ProposeTrackStatus(fStopAndKill);
  fParticleChange->SetProposedKineticEnergy(0.);

  if (photonEnergy <= fIntrinsicLowEnergyLimit)
    {
      fParticleChange->ProposeLocalEnergyDeposit(photonEnergy);
      return ;
    }

  G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();

  // Select randomly one element in the current material
  if (verboseLevel > 2)
    G4cout << "Going to select element in " << couple->GetMaterial()->GetName() << G4endl;

  // atom can be selected efficiently if element selectors are initialised
  const G4Element* anElement =
    SelectRandomAtom(couple,G4Gamma::GammaDefinition(),photonEnergy);
  G4int Z = (G4int) anElement->GetZ();
  if (verboseLevel > 2)
    G4cout << "Selected " << anElement->GetName() << G4endl;
  
  // Select the ionised shell in the current atom according to shell cross sections
  //shellIndex = 0 --> K shell
  //             1-3 --> L shells
  //             4-8 --> M shells
  //             9 --> outer shells cumulatively
  //
  size_t shellIndex = SelectRandomShell(Z,photonEnergy);

  if (verboseLevel > 2)
    G4cout << "Selected shell " << shellIndex << " of element " << anElement->GetName() << G4endl;

  // Retrieve the corresponding identifier and binding energy of the selected shell
  const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();

  //The number of shell cross section possibly reported in the Penelope database 
  //might be different from the number of shells in the G4AtomicTransitionManager
  //(namely, Penelope may contain more shell, especially for very light elements).
  //In order to avoid a warning message from the G4AtomicTransitionManager, I 
  //add this protection. Results are anyway changed, because when G4AtomicTransitionManager
  //has a shellID>maxID, it sets the shellID to the last valid shell. 
  size_t numberOfShells = (size_t) transitionManager->NumberOfShells(Z);
  if (shellIndex >= numberOfShells)
    shellIndex = numberOfShells-1;

  const G4AtomicShell* shell = fTransitionManager->Shell(Z,shellIndex);
  G4double bindingEnergy = shell->BindingEnergy();
  //G4int shellId = shell->ShellId();

  //Penelope considers only K, L and M shells. Cross sections of outer shells are 
  //not included in the Penelope database. If SelectRandomShell() returns 
  //shellIndex = 9, it means that an outer shell was ionized. In this case the 
  //Penelope recipe is to set bindingEnergy = 0 (the energy is entirely assigned 
  //to the electron) and to disregard fluorescence.
  if (shellIndex == 9)
    bindingEnergy = 0.*eV;


  G4double localEnergyDeposit = 0.0;
  G4double cosTheta = 1.0;

  // Primary outcoming electron
  G4double eKineticEnergy = photonEnergy - bindingEnergy;
  
  // There may be cases where the binding energy of the selected shell is > photon energy
  // In such cases do not generate secondaries
  if (eKineticEnergy > 0.)
    {
      // The electron is created
      // Direction sampled from the Sauter distribution
      cosTheta = SampleElectronDirection(eKineticEnergy);
      G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
      G4double phi = twopi * G4UniformRand() ;
      G4double dirx = sinTheta * std::cos(phi);
      G4double diry = sinTheta * std::sin(phi);
      G4double dirz = cosTheta ;
      G4ThreeVector electronDirection(dirx,diry,dirz); //electron direction
      electronDirection.rotateUz(photonDirection);
      G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(), 
                                                           electronDirection, 
                                                           eKineticEnergy);
      fvect->push_back(electron);
    } 
  else    
    bindingEnergy = photonEnergy;


  G4double energyInFluorescence = 0; //testing purposes
  G4double energyInAuger = 0; //testing purposes
 
  //Now, take care of fluorescence, if required. According to the Penelope
  //recipe, I have to skip fluoresence completely if shellIndex == 9 
  //(= sampling of a shell outer than K,L,M)
  if (fAtomDeexcitation && shellIndex<9)
    {      
      G4int index = couple->GetIndex();
      if (fAtomDeexcitation->CheckDeexcitationActiveRegion(index))
        {       
          size_t nBefore = fvect->size();
          fAtomDeexcitation->GenerateParticles(fvect,shell,Z,index);
          size_t nAfter = fvect->size();

          if (nAfter > nBefore) //actual production of fluorescence
            {
              for (size_t j=nBefore;j<nAfter;j++) //loop on products
                {
                  G4double itsEnergy = ((*fvect)[j])->GetKineticEnergy();
                  bindingEnergy -= itsEnergy;
                  if (((*fvect)[j])->GetParticleDefinition() == G4Gamma::Definition())
                    energyInFluorescence += itsEnergy;
                  else if (((*fvect)[j])->GetParticleDefinition() == G4Electron::Definition())
                    energyInAuger += itsEnergy;
                  
                }
            }

        }
    }

  //Residual energy is deposited locally
  localEnergyDeposit += bindingEnergy;
      
  if (localEnergyDeposit < 0)
    {
      G4cout << "WARNING - "
             << "G4PenelopePhotoElectricModel::SampleSecondaries() - Negative energy deposit"
             << G4endl;
      localEnergyDeposit = 0;
    }

  fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);

  if (verboseLevel > 1)
    {
      G4cout << "-----------------------------------------------------------" << G4endl;
      G4cout << "Energy balance from G4PenelopePhotoElectric" << G4endl;
      G4cout << "Selected shell: " << WriteTargetShell(shellIndex) << " of element " << 
        anElement->GetName() << G4endl;
      G4cout << "Incoming photon energy: " << photonEnergy/keV << " keV" << G4endl;
      G4cout << "-----------------------------------------------------------" << G4endl;
      if (eKineticEnergy)
        G4cout << "Outgoing electron " << eKineticEnergy/keV << " keV" << G4endl;
      if (energyInFluorescence)
        G4cout << "Fluorescence x-rays: " << energyInFluorescence/keV << " keV" << G4endl;
      if (energyInAuger)
        G4cout << "Auger electrons: " << energyInAuger/keV << " keV" << G4endl;
      G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
      G4cout << "Total final state: " << 
        (eKineticEnergy+energyInFluorescence+localEnergyDeposit+energyInAuger)/keV << 
        " keV" << G4endl;
      G4cout << "-----------------------------------------------------------" << G4endl;
    }
  if (verboseLevel > 0)
    {
      G4double energyDiff = 
        std::fabs(eKineticEnergy+energyInFluorescence+localEnergyDeposit+energyInAuger-photonEnergy);
      if (energyDiff > 0.05*keV)
        {
          G4cout << "Warning from G4PenelopePhotoElectric: problem with energy conservation: " << 
            (eKineticEnergy+energyInFluorescence+localEnergyDeposit+energyInAuger)/keV 
                 << " keV (final) vs. " << 
            photonEnergy/keV << " keV (initial)" << G4endl;
          G4cout << "-----------------------------------------------------------" << G4endl;
          G4cout << "Energy balance from G4PenelopePhotoElectric" << G4endl;
          G4cout << "Selected shell: " << WriteTargetShell(shellIndex) << " of element " << 
            anElement->GetName() << G4endl;
          G4cout << "Incoming photon energy: " << photonEnergy/keV << " keV" << G4endl;
          G4cout << "-----------------------------------------------------------" << G4endl;
          if (eKineticEnergy)
            G4cout << "Outgoing electron " << eKineticEnergy/keV << " keV" << G4endl;
          if (energyInFluorescence)
            G4cout << "Fluorescence x-rays: " << energyInFluorescence/keV << " keV" << G4endl;
          if (energyInAuger)
            G4cout << "Auger electrons: " << energyInAuger/keV << " keV" << G4endl;
          G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
          G4cout << "Total final state: " << 
            (eKineticEnergy+energyInFluorescence+localEnergyDeposit+energyInAuger)/keV << 
            " keV" << G4endl;
          G4cout << "-----------------------------------------------------------" << G4endl;
        }
    }
}

void G4PenelopePhotoElectricModel::SetVerbosityLevel

(

G4int

lev

)

[inline]

Definition at line 84 of file G4PenelopePhotoElectricModel.hh.

{verboseLevel = lev;};

---

Field Documentation

G4ParticleChangeForGamma* G4PenelopePhotoElectricModel::fParticleChange [protected]

Definition at line 93 of file G4PenelopePhotoElectricModel.hh.

Referenced by Initialise(), and SampleSecondaries().

---

The documentation for this class was generated from the following files:

• G4PenelopePhotoElectricModel.hh
• G4PenelopePhotoElectricModel.cc

---