G4PenelopeIonisationModel Class Reference

#include <G4PenelopeIonisationModel.hh>

Inheritance diagram for G4PenelopeIonisationModel:

Public Member Functions
	G4PenelopeIonisationModel (const G4ParticleDefinition *p=0, const G4String &processName="PenIoni")
virtual	~G4PenelopeIonisationModel ()
virtual void	Initialise (const G4ParticleDefinition *, const G4DataVector &)
virtual G4double	ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double, G4double, G4double, G4double, G4double)
virtual G4double	CrossSectionPerVolume (const G4Material *material, const G4ParticleDefinition *theParticle, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy=DBL_MAX)
virtual void	SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
virtual G4double	ComputeDEDXPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy)
virtual G4double	MinEnergyCut (const G4ParticleDefinition *, const G4MaterialCutsCouple *)
void	SetVerbosityLevel (G4int lev)
G4int	GetVerbosityLevel ()
Protected Attributes
G4ParticleChangeForLoss *	fParticleChange

---

Detailed Description

Definition at line 64 of file G4PenelopeIonisationModel.hh.

---

Constructor & Destructor Documentation

G4PenelopeIonisationModel::G4PenelopeIonisationModel	(	const G4ParticleDefinition *	p = 0,
		const G4String &	processName = "PenIoni"
	)

Definition at line 70 of file G4PenelopeIonisationModel.cc.

References G4PenelopeOscillatorManager::GetOscillatorManager(), G4VEmModel::SetDeexcitationFlag(), and G4VEmModel::SetHighEnergyLimit().

  :G4VEmModel(nam),fParticleChange(0),isInitialised(false),
   fAtomDeexcitation(0),kineticEnergy1(0.*eV),
   cosThetaPrimary(1.0),energySecondary(0.*eV),
   cosThetaSecondary(0.0),targetOscillator(-1),
   theCrossSectionHandler(0)
{
  fIntrinsicLowEnergyLimit = 100.0*eV;
  fIntrinsicHighEnergyLimit = 100.0*GeV;
  //  SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
  SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
  nBins = 200;
  //
  oscManager = G4PenelopeOscillatorManager::GetOscillatorManager();
  //
  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

  // Atomic deexcitation model activated by default
  SetDeexcitationFlag(true);
}

G4PenelopeIonisationModel::~G4PenelopeIonisationModel

(

)

[virtual]

Definition at line 101 of file G4PenelopeIonisationModel.cc.

{
  if (theCrossSectionHandler)
    delete theCrossSectionHandler;
}

---

Member Function Documentation

G4double G4PenelopeIonisationModel::ComputeCrossSectionPerAtom	(	const G4ParticleDefinition *	,
		G4double	,
		G4double	,
		G4double	,
		G4double	,
		G4double
	)			[virtual]

Reimplemented from G4VEmModel.

Definition at line 225 of file G4PenelopeIonisationModel.cc.

References G4cout, and G4endl.

{
  G4cout << "*** G4PenelopeIonisationModel -- WARNING ***" << G4endl;
  G4cout << "Penelope Ionisation model v2008 does not calculate cross section _per atom_ " << G4endl;
  G4cout << "so the result is always zero. For physics values, please invoke " << G4endl;
  G4cout << "GetCrossSectionPerVolume() or GetMeanFreePath() via the G4EmCalculator" << G4endl;
  return 0;
}

G4double G4PenelopeIonisationModel::ComputeDEDXPerVolume	(	const G4Material *	,
		const G4ParticleDefinition *	,
		G4double	kineticEnergy,
		G4double	cutEnergy
	)			[virtual]

Reimplemented from G4VEmModel.

Definition at line 241 of file G4PenelopeIonisationModel.cc.

References G4cout, G4endl, G4PenelopeOscillatorManager::GetAtomsPerMolecule(), G4PenelopeIonisationXSHandler::GetCrossSectionTableForCouple(), G4PenelopeCrossSection::GetSoftStoppingPower(), and G4Material::GetTotNbOfAtomsPerVolume().

{
  // Penelope model v2008 to calculate the stopping power for soft inelastic collisions
  // below the threshold. It makes use of the Generalised Oscillator Strength (GOS)
  // model from
  //  D. Liljequist, J. Phys. D: Appl. Phys. 16 (1983) 1567
  //
  // The stopping power is calculated analytically using the dsigma/dW cross
  // section from the GOS models, which includes separate contributions from
  // distant longitudinal collisions, distant transverse collisions and
  // close collisions. Only the atomic oscillators that come in the play
  // (according to the threshold) are considered for the calculation.
  // Differential cross sections have a different form for e+ and e-.
  //
  // Fermi density correction is calculated analytically according to
  //  U. Fano, Ann. Rev. Nucl. Sci. 13 (1963),1
 
  if (verboseLevel > 3)
    G4cout << "Calling ComputeDEDX() of G4PenelopeIonisationModel" << G4endl;
 

  G4PenelopeCrossSection* theXS = 
    theCrossSectionHandler->GetCrossSectionTableForCouple(theParticle,material,
                                                          cutEnergy);
  G4double sPowerPerMolecule = 0.0;
  if (theXS)
    sPowerPerMolecule = theXS->GetSoftStoppingPower(kineticEnergy);

  
  G4double atomDensity = material->GetTotNbOfAtomsPerVolume();
  G4double atPerMol =  oscManager->GetAtomsPerMolecule(material);
                                                                                        
  G4double moleculeDensity = 0.; 
  if (atPerMol)
    moleculeDensity = atomDensity/atPerMol;
 
  G4double sPowerPerVolume = sPowerPerMolecule*moleculeDensity;

  if (verboseLevel > 2)
    {
      G4cout << "G4PenelopeIonisationModel " << G4endl;
      G4cout << "Stopping power < " << cutEnergy/keV << " keV at " <<
        kineticEnergy/keV << " keV = " << 
        sPowerPerVolume/(keV/mm) << " keV/mm" << G4endl;
    }
  return sPowerPerVolume;
}

G4double G4PenelopeIonisationModel::CrossSectionPerVolume	(	const G4Material *	material,
		const G4ParticleDefinition *	theParticle,
		G4double	kineticEnergy,
		G4double	cutEnergy,
		G4double	maxEnergy = DBL_MAX
	)			[virtual]

Reimplemented from G4VEmModel.

Definition at line 155 of file G4PenelopeIonisationModel.cc.

References G4cout, G4endl, G4PenelopeOscillatorManager::GetAtomsPerMolecule(), G4PenelopeIonisationXSHandler::GetCrossSectionTableForCouple(), G4PenelopeCrossSection::GetHardCrossSection(), G4Material::GetName(), G4Material::GetTotNbOfAtomsPerVolume(), and G4VEmModel::SetupForMaterial().

{
  // Penelope model v2008 to calculate the cross section for inelastic collisions above the
  // threshold. It makes use of the Generalised Oscillator Strength (GOS) model from
  //  D. Liljequist, J. Phys. D: Appl. Phys. 16 (1983) 1567
  //
  // The total cross section is calculated analytically by taking
  // into account the atomic oscillators coming into the play for a given threshold.
  //
  // For incident e- the maximum energy allowed for the delta-rays is energy/2.
  // because particles are undistinghishable.
  //
  // The contribution is splitted in three parts: distant longitudinal collisions,
  // distant transverse collisions and close collisions. Each term is described by
  // its own analytical function.
  // Fermi density correction is calculated analytically according to
  //  U. Fano, Ann. Rev. Nucl. Sci. 13 (1963),1
  //
  if (verboseLevel > 3)
    G4cout << "Calling CrossSectionPerVolume() of G4PenelopeIonisationModel" << G4endl;
 
  SetupForMaterial(theParticle, material, energy);
   
  G4double totalCross = 0.0;
  G4double crossPerMolecule = 0.;

  G4PenelopeCrossSection* theXS = 
    theCrossSectionHandler->GetCrossSectionTableForCouple(theParticle,
                                                          material,
                                                          cutEnergy);

  if (theXS)
    crossPerMolecule = theXS->GetHardCrossSection(energy);

  G4double atomDensity = material->GetTotNbOfAtomsPerVolume();
  G4double atPerMol =  oscManager->GetAtomsPerMolecule(material);
 
  if (verboseLevel > 3)
    G4cout << "Material " << material->GetName() << " has " << atPerMol <<
      "atoms per molecule" << G4endl;

  G4double moleculeDensity = 0.; 
  if (atPerMol)
    moleculeDensity = atomDensity/atPerMol;
 
  G4double crossPerVolume = crossPerMolecule*moleculeDensity;

  if (verboseLevel > 2)
  {
    G4cout << "G4PenelopeIonisationModel " << G4endl;
    G4cout << "Mean free path for delta emission > " << cutEnergy/keV << " keV at " <<
      energy/keV << " keV = " << (1./crossPerVolume)/mm << " mm" << G4endl;
    if (theXS)
      totalCross = (theXS->GetTotalCrossSection(energy))*moleculeDensity;
    G4cout << "Total free path for ionisation (no threshold) at " <<
      energy/keV << " keV = " << (1./totalCross)/mm << " mm" << G4endl;
  }
  return crossPerVolume;
}

G4int G4PenelopeIonisationModel::GetVerbosityLevel

(

)

[inline]

Definition at line 108 of file G4PenelopeIonisationModel.hh.

{return verboseLevel;};

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

Implements G4VEmModel.

Definition at line 109 of file G4PenelopeIonisationModel.cc.

References G4LossTableManager::AtomDeexcitation(), fParticleChange, G4cout, G4endl, G4VEmModel::GetParticleChangeForLoss(), G4VEmModel::HighEnergyLimit(), G4LossTableManager::Instance(), G4VEmModel::LowEnergyLimit(), and G4PenelopeIonisationXSHandler::SetVerboseLevel().

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

  fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
  //Issue warning if the AtomicDeexcitation has not been declared
  if (!fAtomDeexcitation)
    {
      G4cout << G4endl;
      G4cout << "WARNING from G4PenelopeIonisationModel " << 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;
    }

  //Set the number of bins for the tables. 20 points per decade
  nBins = (size_t) (20*std::log10(HighEnergyLimit()/LowEnergyLimit()));
  nBins = std::max(nBins,(size_t)100);

   //Clear and re-build the tables
  if (theCrossSectionHandler)
    {
      delete theCrossSectionHandler;
      theCrossSectionHandler = 0;
    }
  theCrossSectionHandler = new G4PenelopeIonisationXSHandler(nBins);
  theCrossSectionHandler->SetVerboseLevel(verboseLevel);

  if (verboseLevel > 2) {
    G4cout << "Penelope Ionisation model v2008 is initialized " << G4endl
           << "Energy range: "
           << LowEnergyLimit() / keV << " keV - "
           << HighEnergyLimit() / GeV << " GeV. Using " 
           << nBins << " bins." 
           << G4endl;
  }
 
  if(isInitialised) return;
  fParticleChange = GetParticleChangeForLoss();
  isInitialised = true;
}

G4double G4PenelopeIonisationModel::MinEnergyCut	(	const G4ParticleDefinition *	,
		const G4MaterialCutsCouple *
	)			[virtual]

Definition at line 294 of file G4PenelopeIonisationModel.cc.

{
  return fIntrinsicLowEnergyLimit;
}

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

Implements G4VEmModel.

Definition at line 302 of file G4PenelopeIonisationModel.cc.

References G4AtomicShell::BindingEnergy(), G4InuclSpecialFunctions::bindingEnergy(), G4VAtomDeexcitation::CheckDeexcitationActiveRegion(), G4Electron::Definition(), G4Gamma::Definition(), G4Electron::Electron(), FatalException, fParticleChange, G4cout, G4endl, G4Exception(), G4UniformRand, G4VAtomDeexcitation::GenerateParticles(), G4DynamicParticle::GetDefinition(), G4MaterialCutsCouple::GetIndex(), G4DynamicParticle::GetKineticEnergy(), G4MaterialCutsCouple::GetMaterial(), G4DynamicParticle::GetMomentumDirection(), G4PenelopeOscillatorManager::GetOscillatorTableIonisation(), G4ParticleDefinition::GetParticleName(), G4AtomicTransitionManager::Instance(), G4INCL::Math::pi, G4Positron::Positron(), G4VParticleChange::ProposeLocalEnergyDeposit(), G4ParticleChangeForLoss::ProposeMomentumDirection(), G4ParticleChangeForLoss::SetProposedKineticEnergy(), and G4AtomicTransitionManager::Shell().

{
  // Penelope model v2008 to sample the final state following an hard inelastic interaction.
  // It makes use of the Generalised Oscillator Strength (GOS) model from
  //  D. Liljequist, J. Phys. D: Appl. Phys. 16 (1983) 1567
  //
  // The GOS model is used to calculate the individual cross sections for all
  // the atomic oscillators coming in the play, taking into account the three
  // contributions (distant longitudinal collisions, distant transverse collisions and
  // close collisions). Then the target shell and the interaction channel are
  // sampled. Final state of the delta-ray (energy, angle) are generated according
  // to the analytical distributions (dSigma/dW) for the selected interaction
  // channels.
  // For e-, the maximum energy for the delta-ray is initialEnergy/2. (because
  // particles are indistinghusbable), while it is the full initialEnergy for
  // e+.
  // The efficiency of the random sampling algorithm (e.g. for close collisions)
  // decreases when initial and cutoff energy increase (e.g. 87% for 10-keV primary
  // and 1 keV threshold, 99% for 10-MeV primary and 10-keV threshold).
  // Differential cross sections have a different form for e+ and e-.
  //
  // WARNING: The model provides an _average_ description of inelastic collisions.
  // Anyway, the energy spectrum associated to distant excitations of a given
  // atomic shell is approximated as a single resonance. The simulated energy spectra
  // show _unphysical_ narrow peaks at energies that are multiple of the shell
  // resonance energies. The spurious speaks are automatically smoothed out after
  // multiple inelastic collisions.
  //
  // The model determines also the original shell from which the delta-ray is expelled,
  // in order to produce fluorescence de-excitation (from G4DeexcitationManager)
  //
  // Fermi density correction is calculated analytically according to
  //  U. Fano, Ann. Rev. Nucl. Sci. 13 (1963),1

  if (verboseLevel > 3)
    G4cout << "Calling SamplingSecondaries() of G4PenelopeIonisationModel" << G4endl;
 
  G4double kineticEnergy0 = aDynamicParticle->GetKineticEnergy();
  const G4ParticleDefinition* theParticle = aDynamicParticle->GetDefinition();
 
  if (kineticEnergy0 <= fIntrinsicLowEnergyLimit)
  {
    fParticleChange->SetProposedKineticEnergy(0.);
    fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy0);
    return ;
  }

  const G4Material* material = couple->GetMaterial();
  G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableIonisation(material); 

  G4ParticleMomentum particleDirection0 = aDynamicParticle->GetMomentumDirection();
 
  //Initialise final-state variables. The proper values will be set by the methods
  // SampleFinalStateElectron() and SampleFinalStatePositron()
  kineticEnergy1=kineticEnergy0;
  cosThetaPrimary=1.0;
  energySecondary=0.0;
  cosThetaSecondary=1.0;
  targetOscillator = -1;
     
  if (theParticle == G4Electron::Electron())
    SampleFinalStateElectron(material,cutE,kineticEnergy0);
  else if (theParticle == G4Positron::Positron())
    SampleFinalStatePositron(material,cutE,kineticEnergy0);
  else
    {
      G4ExceptionDescription ed;
      ed << "Invalid particle " << theParticle->GetParticleName() << G4endl;
      G4Exception("G4PenelopeIonisationModel::SamplingSecondaries()",
                  "em0001",FatalException,ed);
      
    }
  if (energySecondary == 0) return;

  if (verboseLevel > 3)
  {
     G4cout << "G4PenelopeIonisationModel::SamplingSecondaries() for " << 
        theParticle->GetParticleName() << G4endl;
      G4cout << "Final eKin = " << kineticEnergy1 << " keV" << G4endl;
      G4cout << "Final cosTheta = " << cosThetaPrimary << G4endl;
      G4cout << "Delta-ray eKin = " << energySecondary << " keV" << G4endl;
      G4cout << "Delta-ray cosTheta = " << cosThetaSecondary << G4endl;
      G4cout << "Oscillator: " << targetOscillator << G4endl;
   }
   
  //Update the primary particle
  G4double sint = std::sqrt(1. - cosThetaPrimary*cosThetaPrimary);
  G4double phiPrimary  = twopi * G4UniformRand();
  G4double dirx = sint * std::cos(phiPrimary);
  G4double diry = sint * std::sin(phiPrimary);
  G4double dirz = cosThetaPrimary;
 
  G4ThreeVector electronDirection1(dirx,diry,dirz);
  electronDirection1.rotateUz(particleDirection0);
   
  if (kineticEnergy1 > 0)
    {
      fParticleChange->ProposeMomentumDirection(electronDirection1);
      fParticleChange->SetProposedKineticEnergy(kineticEnergy1);
    }
  else    
    fParticleChange->SetProposedKineticEnergy(0.);
    
  
  //Generate the delta ray
  G4double ionEnergyInPenelopeDatabase = 
    (*theTable)[targetOscillator]->GetIonisationEnergy();
  
  //Now, try to handle fluorescence
  //Notice: merged levels are indicated with Z=0 and flag=30
  G4int shFlag = (*theTable)[targetOscillator]->GetShellFlag(); 
  G4int Z = (G4int) (*theTable)[targetOscillator]->GetParentZ();

  //initialize here, then check photons created by Atomic-Deexcitation, and the final state e-
  const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
  G4double bindingEnergy = 0.*eV;
  //G4int shellId = 0;

  const G4AtomicShell* shell = 0;
  //Real level
  if (Z > 0 && shFlag<30)
    {
      shell = transitionManager->Shell(Z,shFlag-1);
      bindingEnergy = shell->BindingEnergy();
      //shellId = shell->ShellId();
    }

  //correct the energySecondary to account for the fact that the Penelope 
  //database of ionisation energies is in general (slightly) different 
  //from the fluorescence database used in Geant4.
  energySecondary += ionEnergyInPenelopeDatabase-bindingEnergy;
  
  G4double localEnergyDeposit = bindingEnergy;
  //testing purposes only
  G4double energyInFluorescence = 0;
  G4double energyInAuger = 0; 

  if (energySecondary < 0)
    {
      //It means that there was some problem/mismatch between the two databases. 
      //In this case, the available energy is ok to excite the level according 
      //to the Penelope database, but not according to the Geant4 database
      //Full residual energy is deposited locally
      localEnergyDeposit += energySecondary;
      energySecondary = 0.0;
    }

  //Notice: shell might be NULL (invalid!) if shFlag=30. Must be protected
  //Now, take care of fluorescence, if required
  if (fAtomDeexcitation && shell)
    {
      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();
                  localEnergyDeposit -= itsEnergy;
                  if (((*fvect)[j])->GetParticleDefinition() == G4Gamma::Definition())
                    energyInFluorescence += itsEnergy;
                  else if (((*fvect)[j])->GetParticleDefinition() == G4Electron::Definition())
                    energyInAuger += itsEnergy;
                }
            }
        }
    }

  // Generate the delta ray --> to be done only if above cut
  if (energySecondary > cutE)
    {
      G4DynamicParticle* electron = 0;
      G4double sinThetaE = std::sqrt(1.-cosThetaSecondary*cosThetaSecondary);
      G4double phiEl = phiPrimary+pi; //pi with respect to the primary electron/positron
      G4double xEl = sinThetaE * std::cos(phiEl);
      G4double yEl = sinThetaE * std::sin(phiEl);
      G4double zEl = cosThetaSecondary;
      G4ThreeVector eDirection(xEl,yEl,zEl); //electron direction
      eDirection.rotateUz(particleDirection0);
      electron = new G4DynamicParticle (G4Electron::Electron(),
                                        eDirection,energySecondary) ;
      fvect->push_back(electron);
    }
  else
    {
      localEnergyDeposit += energySecondary;
      energySecondary = 0;
    }

  if (localEnergyDeposit < 0)
    {
      G4cout << "WARNING-" 
             << "G4PenelopeIonisationModel::SampleSecondaries - Negative energy deposit"
             << G4endl;
      localEnergyDeposit=0.;
    }
  fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);

  if (verboseLevel > 1)
    {
      G4cout << "-----------------------------------------------------------" << G4endl;
      G4cout << "Energy balance from G4PenelopeIonisation" << G4endl;
      G4cout << "Incoming primary energy: " << kineticEnergy0/keV << " keV" << G4endl;
      G4cout << "-----------------------------------------------------------" << G4endl;
      G4cout << "Outgoing primary energy: " << kineticEnergy1/keV << " keV" << G4endl;
      G4cout << "Delta ray " << energySecondary/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: " << (energySecondary+energyInFluorescence+kineticEnergy1+
                                          localEnergyDeposit+energyInAuger)/keV <<
        " keV" << G4endl;
      G4cout << "-----------------------------------------------------------" << G4endl;
    }
      
  if (verboseLevel > 0)
    {
      G4double energyDiff = std::fabs(energySecondary+energyInFluorescence+kineticEnergy1+
                                      localEnergyDeposit+energyInAuger-kineticEnergy0);
      if (energyDiff > 0.05*keV)
        G4cout << "Warning from G4PenelopeIonisation: problem with energy conservation: " <<  
          (energySecondary+energyInFluorescence+kineticEnergy1+localEnergyDeposit+energyInAuger)/keV <<
          " keV (final) vs. " <<
          kineticEnergy0/keV << " keV (initial)" << G4endl;      
    }
    
}

void G4PenelopeIonisationModel::SetVerbosityLevel

(

G4int

lev

)

[inline]

Definition at line 107 of file G4PenelopeIonisationModel.hh.

{verboseLevel = lev;};

---

Field Documentation

G4ParticleChangeForLoss* G4PenelopeIonisationModel::fParticleChange [protected]

Definition at line 108 of file G4PenelopeIonisationModel.hh.

Referenced by Initialise(), and SampleSecondaries().

---

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

• G4PenelopeIonisationModel.hh
• G4PenelopeIonisationModel.cc

---