G4PenelopeIonisationModel Class Reference

#include <G4PenelopeIonisationModel.hh>

Inheritance diagram for G4PenelopeIonisationModel:

Collaboration 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 void	InitialiseLocal (const G4ParticleDefinition *, G4VEmModel *)
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 ()
Public Member Functions inherited from G4VEmModel
	G4VEmModel (const G4String &nam)
virtual	~G4VEmModel ()
virtual void	InitialiseForMaterial (const G4ParticleDefinition *, const G4Material *)
virtual void	InitialiseForElement (const G4ParticleDefinition *, G4int Z)
virtual G4double	GetPartialCrossSection (const G4Material *, G4int, const G4ParticleDefinition *, G4double)
virtual G4double	ComputeCrossSectionPerShell (const G4ParticleDefinition *, G4int Z, G4int shellIdx, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
virtual G4double	ChargeSquareRatio (const G4Track &)
virtual G4double	GetChargeSquareRatio (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
virtual G4double	GetParticleCharge (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
virtual void	StartTracking (G4Track *)
virtual void	CorrectionsAlongStep (const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double &eloss, G4double &niel, G4double length)
virtual G4double	Value (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
virtual G4double	MinPrimaryEnergy (const G4Material *, const G4ParticleDefinition *, G4double cut=0.0)
virtual void	SetupForMaterial (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
virtual void	DefineForRegion (const G4Region *)
virtual void	ModelDescription (std::ostream &outFile) const
void	InitialiseElementSelectors (const G4ParticleDefinition *, const G4DataVector &)
std::vector< G4EmElementSelector * > *	GetElementSelectors ()
void	SetElementSelectors (std::vector< G4EmElementSelector *> *)
G4double	ComputeDEDX (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
G4double	CrossSection (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
G4double	ComputeMeanFreePath (const G4ParticleDefinition *, G4double kineticEnergy, const G4Material *, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
G4double	ComputeCrossSectionPerAtom (const G4ParticleDefinition *, const G4Element *, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
G4int	SelectIsotopeNumber (const G4Element *)
const G4Element *	SelectRandomAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
const G4Element *	SelectRandomAtom (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
G4int	SelectRandomAtomNumber (const G4Material *)
void	SetParticleChange (G4VParticleChange *, G4VEmFluctuationModel *f=0)
void	SetCrossSectionTable (G4PhysicsTable *, G4bool isLocal)
G4ElementData *	GetElementData ()
G4PhysicsTable *	GetCrossSectionTable ()
G4VEmFluctuationModel *	GetModelOfFluctuations ()
G4VEmAngularDistribution *	GetAngularDistribution ()
void	SetAngularDistribution (G4VEmAngularDistribution *)
G4double	HighEnergyLimit () const
G4double	LowEnergyLimit () const
G4double	HighEnergyActivationLimit () const
G4double	LowEnergyActivationLimit () const
G4double	PolarAngleLimit () const
G4double	SecondaryThreshold () const
G4bool	LPMFlag () const
G4bool	DeexcitationFlag () const
G4bool	ForceBuildTableFlag () const
G4bool	UseAngularGeneratorFlag () const
void	SetAngularGeneratorFlag (G4bool)
void	SetHighEnergyLimit (G4double)
void	SetLowEnergyLimit (G4double)
void	SetActivationHighEnergyLimit (G4double)
void	SetActivationLowEnergyLimit (G4double)
G4bool	IsActive (G4double kinEnergy)
void	SetPolarAngleLimit (G4double)
void	SetSecondaryThreshold (G4double)
void	SetLPMFlag (G4bool val)
void	SetDeexcitationFlag (G4bool val)
void	SetForceBuildTable (G4bool val)
void	SetMasterThread (G4bool val)
G4bool	IsMaster () const
G4double	MaxSecondaryKinEnergy (const G4DynamicParticle *dynParticle)
const G4String &	GetName () const
void	SetCurrentCouple (const G4MaterialCutsCouple *)
const G4Element *	GetCurrentElement () const
const G4Isotope *	GetCurrentIsotope () const
G4bool	IsLocked () const
void	SetLocked (G4bool)

Protected Attributes
G4ParticleChangeForLoss *	fParticleChange
const G4ParticleDefinition *	fParticle
Protected Attributes inherited from G4VEmModel
G4ElementData *	fElementData
G4VParticleChange *	pParticleChange
G4PhysicsTable *	xSectionTable
const std::vector< G4double > *	theDensityFactor
const std::vector< G4int > *	theDensityIdx
size_t	idxTable

Private Member Functions
G4PenelopeIonisationModel &	operator= (const G4PenelopeIonisationModel &right)
	G4PenelopeIonisationModel (const G4PenelopeIonisationModel &)
void	SetParticle (const G4ParticleDefinition *)
void	SampleFinalStateElectron (const G4Material *, G4double cutEnergy, G4double kineticEnergy)
void	SampleFinalStatePositron (const G4Material *, G4double cutEnergy, G4double kineticEnergy)

Private Attributes
G4double	fIntrinsicLowEnergyLimit
G4double	fIntrinsicHighEnergyLimit
G4int	verboseLevel
G4bool	isInitialised
G4VAtomDeexcitation *	fAtomDeexcitation
G4bool	fPIXEflag
G4double	kineticEnergy1
G4double	cosThetaPrimary
G4double	energySecondary
G4double	cosThetaSecondary
G4int	targetOscillator
G4PenelopeOscillatorManager *	oscManager
G4PenelopeIonisationXSHandler *	theCrossSectionHandler
size_t	nBins
G4bool	fLocalTable

Additional Inherited Members
Protected Member Functions inherited from G4VEmModel
G4ParticleChangeForLoss *	GetParticleChangeForLoss ()
G4ParticleChangeForGamma *	GetParticleChangeForGamma ()
virtual G4double	MaxSecondaryEnergy (const G4ParticleDefinition *, G4double kineticEnergy)
const G4MaterialCutsCouple *	CurrentCouple () const
void	SetCurrentElement (const G4Element *)
Static Protected Attributes inherited from G4VEmModel
static const G4double	inveplus = 1.0/CLHEP::eplus

Detailed Description

Definition at line 66 of file G4PenelopeIonisationModel.hh.

Constructor & Destructor Documentation

G4PenelopeIonisationModel() [1/2]

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

Definition at line 73 of file G4PenelopeIonisationModel.cc.

  :G4VEmModel(nam),fParticleChange(0),fParticle(0),isInitialised(false),
   fAtomDeexcitation(0),fPIXEflag(false),kineticEnergy1(0.*eV),
   cosThetaPrimary(1.0),energySecondary(0.*eV),
   cosThetaSecondary(0.0),targetOscillator(-1),
   theCrossSectionHandler(0),fLocalTable(false)
{
  fIntrinsicLowEnergyLimit = 100.0*eV;
  fIntrinsicHighEnergyLimit = 100.0*GeV;
  //  SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
  SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
  nBins = 200;

  if (part)
    SetParticle(part);

  //
  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);
}

Here is the call graph for this function:

~G4PenelopeIonisationModel()

G4PenelopeIonisationModel::~G4PenelopeIonisationModel ( )

virtual

Definition at line 108 of file G4PenelopeIonisationModel.cc.

{
  if (IsMaster() || fLocalTable)
    {
      if (theCrossSectionHandler)
    delete theCrossSectionHandler;
    }
}

Here is the call graph for this function:

G4PenelopeIonisationModel() [2/2]

G4PenelopeIonisationModel::G4PenelopeIonisationModel ( const G4PenelopeIonisationModel &  )

private

Member Function Documentation

ComputeCrossSectionPerAtom()

G4double G4PenelopeIonisationModel::ComputeCrossSectionPerAtom ( const G4ParticleDefinition *  , G4double  , G4double  , G4double  , G4double  , G4double    )

virtual

Reimplemented from G4VEmModel.

Definition at line 356 of file G4PenelopeIonisationModel.cc.

{
  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;
}

ComputeDEDXPerVolume()

G4double G4PenelopeIonisationModel::ComputeDEDXPerVolume ( const G4Material *  material, const G4ParticleDefinition *  theParticle, G4double  kineticEnergy, G4double  cutEnergy  )

virtual

Reimplemented from G4VEmModel.

Definition at line 372 of file G4PenelopeIonisationModel.cc.

{
  // 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;
 
  //Either Initialize() was not called, or we are in a slave and InitializeLocal() was 
  //not invoked
  if (!theCrossSectionHandler)
    {
      //create a **thread-local** version of the table. Used only for G4EmCalculator and 
      //Unit Tests
      fLocalTable = true;
      theCrossSectionHandler = new G4PenelopeIonisationXSHandler(nBins);    
    }

  const G4PenelopeCrossSection* theXS = 
    theCrossSectionHandler->GetCrossSectionTableForCouple(theParticle,material,
                              cutEnergy);
  
  if (!theXS)
    {
      //If we are here, it means that Initialize() was inkoved, but the MaterialTable was 
      //not filled up. This can happen in a UnitTest or via G4EmCalculator
      if (verboseLevel > 0)
    {   
      //Issue a G4Exception (warning) only in verbose mode
      G4ExceptionDescription ed;
      ed << "Unable to retrieve the cross section table for " << theParticle->GetParticleName() <<
        " in " << material->GetName() << ", cut = " << cutEnergy/keV << " keV " << G4endl;
      ed << "This can happen only in Unit Tests or via G4EmCalculator" << G4endl;
      G4Exception("G4PenelopeIonisationModel::ComputeDEDXPerVolume()",
              "em2038",JustWarning,ed);
    }
      //protect file reading via autolock
      G4AutoLock lock(&PenelopeIonisationModelMutex);
      theCrossSectionHandler->BuildXSTable(material,cutEnergy,theParticle);
      lock.unlock();
      //now it should be ok
      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;
}

Here is the call graph for this function:

CrossSectionPerVolume()

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

virtual

Reimplemented from G4VEmModel.

Definition at line 251 of file G4PenelopeIonisationModel.cc.

{
  // 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.;

  //Either Initialize() was not called, or we are in a slave and InitializeLocal() was 
  //not invoked
  if (!theCrossSectionHandler)
    {
      //create a **thread-local** version of the table. Used only for G4EmCalculator and 
      //Unit Tests
      fLocalTable = true;
      theCrossSectionHandler = new G4PenelopeIonisationXSHandler(nBins);    
    }
  
  const G4PenelopeCrossSection* theXS = 
    theCrossSectionHandler->GetCrossSectionTableForCouple(theParticle,
                              material,
                              cutEnergy);
  if (!theXS)
    {
      //If we are here, it means that Initialize() was inkoved, but the MaterialTable was 
      //not filled up. This can happen in a UnitTest or via G4EmCalculator
      if (verboseLevel > 0)
    {
      //Issue a G4Exception (warning) only in verbose mode
      G4ExceptionDescription ed;
      ed << "Unable to retrieve the cross section table for " << theParticle->GetParticleName() << 
        " in " << material->GetName() << ", cut = " << cutEnergy/keV << " keV " << G4endl;
      ed << "This can happen only in Unit Tests or via G4EmCalculator" << G4endl;
      G4Exception("G4PenelopeIonisationModel::CrossSectionPerVolume()",
              "em2038",JustWarning,ed);
    }
      //protect file reading via autolock
      G4AutoLock lock(&PenelopeIonisationModelMutex);
      theCrossSectionHandler->BuildXSTable(material,cutEnergy,theParticle);
      lock.unlock();
      //now it should be ok
      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;
}

Here is the call graph for this function:

GetVerbosityLevel()

G4int G4PenelopeIonisationModel::GetVerbosityLevel ( )

inline

Definition at line 113 of file G4PenelopeIonisationModel.hh.

{return verboseLevel;};

Initialise()

void G4PenelopeIonisationModel::Initialise ( const G4ParticleDefinition *  particle, const G4DataVector &  theCuts  )

virtual

Implements G4VEmModel.

Definition at line 119 of file G4PenelopeIonisationModel.cc.

{
  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;
    }

  if (fAtomDeexcitation)
    fPIXEflag = fAtomDeexcitation->IsPIXEActive();

  //If the PIXE flag is active, the PIXE interface will take care of the 
  //atomic de-excitation. The model does not need to do that.
  //Issue warnings here
  if (fPIXEflag && IsMaster() && particle==G4Electron::Electron())
    {
      G4String theModel = G4EmParameters::Instance()->PIXEElectronCrossSectionModel();
      G4cout << "======================================================================" << G4endl;
      G4cout << "The G4PenelopeIonisationModel is being used with the PIXE flag ON." << G4endl;
      G4cout << "Atomic de-excitation will be produced statistically by the PIXE " << G4endl;
      G4cout << "interface by using the shell cross section --> " << theModel << G4endl;
      G4cout << "The built-in model procedure for atomic de-excitation is disabled. " << G4endl;
      G4cout << "*Please be sure this is intended*, or disable PIXE by" << G4endl;
      G4cout << "/process/em/pixe false" << G4endl;
      /*
      if (theModel != "Penelope")
    {
      ed << "To use the PIXE interface with the Penelope-based shell cross sections" << G4endl;
      ed << "/process/em/pixeElecXSmodel Penelope" << G4endl;
      
    }
      */
      G4cout << "======================================================================" << G4endl;

    }
 
   

  SetParticle(particle);

  //Only the master model creates/manages the tables. All workers get the 
  //pointer to the table, and use it as readonly
  if (IsMaster() && particle == fParticle)
    {

      //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);

      //Build tables for all materials
      G4ProductionCutsTable* theCoupleTable = 
    G4ProductionCutsTable::GetProductionCutsTable();      
      for (size_t i=0;i<theCoupleTable->GetTableSize();i++)
    {
      const G4Material* theMat = 
        theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
      //Forces the building of the cross section tables
      theCrossSectionHandler->BuildXSTable(theMat,theCuts.at(i),particle,
                           IsMaster());
     
    }

      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;

  return;
}

Here is the call graph for this function:

InitialiseLocal()

void G4PenelopeIonisationModel::InitialiseLocal ( const G4ParticleDefinition *  part, G4VEmModel *  masterModel  )

virtual

Reimplemented from G4VEmModel.

Definition at line 219 of file G4PenelopeIonisationModel.cc.

{
  if (verboseLevel > 3)
    G4cout << "Calling  G4PenelopeIonisationModel::InitialiseLocal()" << G4endl;
 
  //
  //Check that particle matches: one might have multiple master models (e.g. 
  //for e+ and e-).
  //
  if (part == fParticle)
    {
      //Get the const table pointers from the master to the workers
      const G4PenelopeIonisationModel* theModel = 
    static_cast<G4PenelopeIonisationModel*> (masterModel);
      
      //Copy pointers to the data tables
      theCrossSectionHandler = theModel->theCrossSectionHandler;
      
      //copy data
      nBins = theModel->nBins;
      
      //Same verbosity for all workers, as the master
      verboseLevel = theModel->verboseLevel;
    }

  return;
}

MinEnergyCut()

G4double G4PenelopeIonisationModel::MinEnergyCut ( const G4ParticleDefinition *  , const G4MaterialCutsCouple *    )

virtual

Reimplemented from G4VEmModel.

Definition at line 460 of file G4PenelopeIonisationModel.cc.

{
  return fIntrinsicLowEnergyLimit;
}

operator=()

G4PenelopeIonisationModel& G4PenelopeIonisationModel::operator= ( const G4PenelopeIonisationModel &  right )

private

SampleFinalStateElectron()

void G4PenelopeIonisationModel::SampleFinalStateElectron ( const G4Material *  mat, G4double  cutEnergy, G4double  kineticEnergy  )

private

Definition at line 710 of file G4PenelopeIonisationModel.cc.

{
  // This method sets the final ionisation parameters
  // kineticEnergy1, cosThetaPrimary (= updates of the primary e-)
  // energySecondary, cosThetaSecondary (= info of delta-ray)
  // targetOscillator (= ionised oscillator)
  //
  // The method implements SUBROUTINE EINa of Penelope
  //

  G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableIonisation(mat);
  size_t numberOfOscillators = theTable->size();
  const G4PenelopeCrossSection* theXS = 
    theCrossSectionHandler->GetCrossSectionTableForCouple(G4Electron::Electron(),mat,
                              cutEnergy);
  G4double delta = theCrossSectionHandler->GetDensityCorrection(mat,kineticEnergy);
 
  // Selection of the active oscillator
  G4double TST = G4UniformRand();
  targetOscillator = numberOfOscillators-1; //initialization, last oscillator
  G4double XSsum = 0.;

  for (size_t i=0;i<numberOfOscillators-1;i++)
    {
      /* testing purposes
      G4cout << "sampling: " << i << " " << XSsum << " " << TST << " " << 
    theXS->GetShellCrossSection(i,kineticEnergy) << " " << 
    theXS->GetNormalizedShellCrossSection(i,kineticEnergy) << " " << 
    mat->GetName() << G4endl;
      */
      XSsum += theXS->GetNormalizedShellCrossSection(i,kineticEnergy); 
        
      if (XSsum > TST)
    {
      targetOscillator = (G4int) i;
      break;
    }
    }

  
  if (verboseLevel > 3)
    {
      G4cout << "SampleFinalStateElectron: sampled oscillator #" << targetOscillator << "." << G4endl;
      G4cout << "Ionisation energy: " << (*theTable)[targetOscillator]->GetIonisationEnergy()/eV << 
    " eV " << G4endl;
      G4cout << "Resonance energy: : " << (*theTable)[targetOscillator]->GetResonanceEnergy()/eV << " eV "
         << G4endl;
    }
  //Constants
  G4double rb = kineticEnergy + 2.0*electron_mass_c2;
  G4double gam = 1.0+kineticEnergy/electron_mass_c2;
  G4double gam2 = gam*gam;
  G4double beta2 = (gam2-1.0)/gam2;
  G4double amol = ((gam-1.0)/gam)*((gam-1.0)/gam);

  //Partial cross section of the active oscillator
  G4double resEne = (*theTable)[targetOscillator]->GetResonanceEnergy();
  G4double invResEne = 1.0/resEne;
  G4double ionEne = (*theTable)[targetOscillator]->GetIonisationEnergy();
  G4double cutoffEne = (*theTable)[targetOscillator]->GetCutoffRecoilResonantEnergy();
  G4double XHDL = 0.;
  G4double XHDT = 0.;
  G4double QM = 0.;
  G4double cps = 0.;
  G4double cp = 0.;

  //Distant excitations
  if (resEne > cutEnergy && resEne < kineticEnergy)
    {
      cps = kineticEnergy*rb;
      cp = std::sqrt(cps);
      G4double XHDT0 = std::max(std::log(gam2)-beta2-delta,0.);
      if (resEne > 1.0e-6*kineticEnergy)
    {
      G4double cpp = std::sqrt((kineticEnergy-resEne)*(kineticEnergy-resEne+2.0*electron_mass_c2));
      QM = std::sqrt((cp-cpp)*(cp-cpp)+electron_mass_c2*electron_mass_c2)-electron_mass_c2;
    }
      else
    {
      QM = resEne*resEne/(beta2*2.0*electron_mass_c2);
      QM *= (1.0-QM*0.5/electron_mass_c2);
    }
      if (QM < cutoffEne)
    {
      XHDL = std::log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2)))
        *invResEne;
      XHDT = XHDT0*invResEne;     
    }
      else
    {
      QM = cutoffEne;
      XHDL = 0.;
      XHDT = 0.;
    }
    }
  else
    {
      QM = cutoffEne;
      cps = 0.;
      cp = 0.;
      XHDL = 0.;
      XHDT = 0.;
    }

  //Close collisions
  G4double EE = kineticEnergy + ionEne;
  G4double wmaxc = 0.5*EE;
  G4double wcl = std::max(cutEnergy,cutoffEne);
  G4double rcl = wcl/EE;
  G4double XHC = 0.;
  if (wcl < wmaxc)
    {
      G4double rl1 = 1.0-rcl;
      G4double rrl1 = 1.0/rl1;
      XHC = (amol*(0.5-rcl)+1.0/rcl-rrl1+
         (1.0-amol)*std::log(rcl*rrl1))/EE;
    }

  //Total cross section per molecule for the active shell, in cm2
  G4double XHTOT = XHC + XHDL + XHDT;

  //very small cross section, do nothing
  if (XHTOT < 1.e-14*barn)
    {
      kineticEnergy1=kineticEnergy;
      cosThetaPrimary=1.0;
      energySecondary=0.0;
      cosThetaSecondary=1.0;
      targetOscillator = numberOfOscillators-1;
      return;
    }
  
  //decide which kind of interaction we'll have
  TST = XHTOT*G4UniformRand();
  

  // Hard close collision
  G4double TS1 = XHC;
  
  if (TST < TS1)
    {
      G4double A = 5.0*amol;
      G4double ARCL = A*0.5*rcl;
      G4double rk=0.;
      G4bool loopAgain = false;
      do
    {
      loopAgain = false;
      G4double fb = (1.0+ARCL)*G4UniformRand();
      if (fb < 1)
        rk = rcl/(1.0-fb*(1.0-(rcl+rcl)));       
      else
        rk = rcl + (fb-1.0)*(0.5-rcl)/ARCL;
      G4double rk2 = rk*rk;
      G4double rkf = rk/(1.0-rk);
      G4double phi = 1.0+rkf*rkf-rkf+amol*(rk2+rkf);
      if (G4UniformRand()*(1.0+A*rk2) > phi)
        loopAgain = true;
    }while(loopAgain);
      //energy and scattering angle (primary electron)
      G4double deltaE = rk*EE;
      
      kineticEnergy1 = kineticEnergy - deltaE;
      cosThetaPrimary = std::sqrt(kineticEnergy1*rb/(kineticEnergy*(rb-deltaE)));
      //energy and scattering angle of the delta ray
      energySecondary = deltaE - ionEne; //subtract ionisation energy
      cosThetaSecondary= std::sqrt(deltaE*rb/(kineticEnergy*(deltaE+2.0*electron_mass_c2)));
      if (verboseLevel > 3)
    G4cout << "SampleFinalStateElectron: sampled close collision " << G4endl;
      return;                    
    }

  //Hard distant longitudinal collisions
  TS1 += XHDL;
  G4double deltaE = resEne;
  kineticEnergy1 = kineticEnergy - deltaE;

  if (TST < TS1)
    {
      G4double QS = QM/(1.0+QM*0.5/electron_mass_c2);
      G4double Q = QS/(std::pow((QS/cutoffEne)*(1.0+cutoffEne*0.5/electron_mass_c2),G4UniformRand())
               - (QS*0.5/electron_mass_c2));
      G4double QTREV = Q*(Q+2.0*electron_mass_c2);
      G4double cpps = kineticEnergy1*(kineticEnergy1+2.0*electron_mass_c2);
      cosThetaPrimary = (cpps+cps-QTREV)/(2.0*cp*std::sqrt(cpps));
      if (cosThetaPrimary > 1.) 
    cosThetaPrimary = 1.0;
      //energy and emission angle of the delta ray
      energySecondary = deltaE - ionEne;
      cosThetaSecondary = 0.5*(deltaE*(kineticEnergy+rb-deltaE)+QTREV)/std::sqrt(cps*QTREV);
      if (cosThetaSecondary > 1.0)
    cosThetaSecondary = 1.0;
      if (verboseLevel > 3)
    G4cout << "SampleFinalStateElectron: sampled distant longitudinal collision " << G4endl;
      return;      
    }

  //Hard distant transverse collisions
  cosThetaPrimary = 1.0;
  //energy and emission angle of the delta ray
  energySecondary = deltaE - ionEne;
  cosThetaSecondary = 0.5;
  if (verboseLevel > 3)
    G4cout << "SampleFinalStateElectron: sampled distant transverse collision " << G4endl;

  return;
}

Here is the call graph for this function:

Here is the caller graph for this function:

SampleFinalStatePositron()

void G4PenelopeIonisationModel::SampleFinalStatePositron ( const G4Material *  mat, G4double  cutEnergy, G4double  kineticEnergy  )

private

Definition at line 923 of file G4PenelopeIonisationModel.cc.

{
  // This method sets the final ionisation parameters
  // kineticEnergy1, cosThetaPrimary (= updates of the primary e-)
  // energySecondary, cosThetaSecondary (= info of delta-ray)
  // targetOscillator (= ionised oscillator)
  //
  // The method implements SUBROUTINE PINa of Penelope
  // 
 
  G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableIonisation(mat);
  size_t numberOfOscillators = theTable->size();
  const G4PenelopeCrossSection* theXS = 
    theCrossSectionHandler->GetCrossSectionTableForCouple(G4Positron::Positron(),mat,
                              cutEnergy);
  G4double delta = theCrossSectionHandler->GetDensityCorrection(mat,kineticEnergy);

  // Selection of the active oscillator
  G4double TST = G4UniformRand();
  targetOscillator = numberOfOscillators-1; //initialization, last oscillator
  G4double XSsum = 0.;
  for (size_t i=0;i<numberOfOscillators-1;i++)
    {
      XSsum += theXS->GetNormalizedShellCrossSection(i,kineticEnergy);      
      if (XSsum > TST)
    {
      targetOscillator = (G4int) i;
      break;
    }
    }

  if (verboseLevel > 3)
    {
      G4cout << "SampleFinalStatePositron: sampled oscillator #" << targetOscillator << "." << G4endl;
      G4cout << "Ionisation energy: " << (*theTable)[targetOscillator]->GetIonisationEnergy()/eV 
         << " eV " << G4endl; 
      G4cout << "Resonance energy: : " << (*theTable)[targetOscillator]->GetResonanceEnergy()/eV 
         << " eV " << G4endl;
    }

  //Constants
  G4double rb = kineticEnergy + 2.0*electron_mass_c2;
  G4double gam = 1.0+kineticEnergy/electron_mass_c2;
  G4double gam2 = gam*gam;
  G4double beta2 = (gam2-1.0)/gam2;
  G4double g12 = (gam+1.0)*(gam+1.0);
  G4double amol = ((gam-1.0)/gam)*((gam-1.0)/gam);
  //Bhabha coefficients
  G4double bha1 = amol*(2.0*g12-1.0)/(gam2-1.0);
  G4double bha2 = amol*(3.0+1.0/g12);
  G4double bha3 = amol*2.0*gam*(gam-1.0)/g12;
  G4double bha4 = amol*(gam-1.0)*(gam-1.0)/g12;

  //
  //Partial cross section of the active oscillator
  //
  G4double resEne = (*theTable)[targetOscillator]->GetResonanceEnergy();
  G4double invResEne = 1.0/resEne;
  G4double ionEne = (*theTable)[targetOscillator]->GetIonisationEnergy();
  G4double cutoffEne = (*theTable)[targetOscillator]->GetCutoffRecoilResonantEnergy();

  G4double XHDL = 0.;
  G4double XHDT = 0.;
  G4double QM = 0.;
  G4double cps = 0.;
  G4double cp = 0.;

  //Distant excitations XS (same as for electrons)
  if (resEne > cutEnergy && resEne < kineticEnergy)
    {
      cps = kineticEnergy*rb;
      cp = std::sqrt(cps);
      G4double XHDT0 = std::max(std::log(gam2)-beta2-delta,0.);
      if (resEne > 1.0e-6*kineticEnergy)
    {
      G4double cpp = std::sqrt((kineticEnergy-resEne)*(kineticEnergy-resEne+2.0*electron_mass_c2));
      QM = std::sqrt((cp-cpp)*(cp-cpp)+electron_mass_c2*electron_mass_c2)-electron_mass_c2;
    }
      else
    {
      QM = resEne*resEne/(beta2*2.0*electron_mass_c2);
      QM *= (1.0-QM*0.5/electron_mass_c2);
    }
      if (QM < cutoffEne)
    {
      XHDL = std::log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2)))
        *invResEne;
      XHDT = XHDT0*invResEne;     
    }
      else
    {
      QM = cutoffEne;
      XHDL = 0.;
      XHDT = 0.;
    }
    }
  else
    {
      QM = cutoffEne;
      cps = 0.;
      cp = 0.;
      XHDL = 0.;
      XHDT = 0.;
    }
  //Close collisions (Bhabha)
  G4double wmaxc = kineticEnergy;
  G4double wcl = std::max(cutEnergy,cutoffEne);
  G4double rcl = wcl/kineticEnergy;
  G4double XHC = 0.;
  if (wcl < wmaxc)
    {
      G4double rl1 = 1.0-rcl;
      XHC = ((1.0/rcl-1.0)+bha1*std::log(rcl)+bha2*rl1
         + (bha3/2.0)*(rcl*rcl-1.0) 
         + (bha4/3.0)*(1.0-rcl*rcl*rcl))/kineticEnergy;
    }

  //Total cross section per molecule for the active shell, in cm2
  G4double XHTOT = XHC + XHDL + XHDT;

  //very small cross section, do nothing
  if (XHTOT < 1.e-14*barn)
    {
      kineticEnergy1=kineticEnergy;
      cosThetaPrimary=1.0;
      energySecondary=0.0;
      cosThetaSecondary=1.0;
      targetOscillator = numberOfOscillators-1;
      return;
    }

  //decide which kind of interaction we'll have
  TST = XHTOT*G4UniformRand();
  
  // Hard close collision
  G4double TS1 = XHC;
  if (TST < TS1)
    {
      G4double rl1 = 1.0-rcl;
      G4double rk=0.;
      G4bool loopAgain = false;
      do
    {
      loopAgain = false;
      rk = rcl/(1.0-G4UniformRand()*rl1);
      G4double phi = 1.0-rk*(bha1-rk*(bha2-rk*(bha3-bha4*rk)));
      if (G4UniformRand() > phi)
        loopAgain = true;
    }while(loopAgain);
      //energy and scattering angle (primary electron)
      G4double deltaE = rk*kineticEnergy;      
      kineticEnergy1 = kineticEnergy - deltaE;
      cosThetaPrimary = std::sqrt(kineticEnergy1*rb/(kineticEnergy*(rb-deltaE)));
      //energy and scattering angle of the delta ray
      energySecondary = deltaE - ionEne; //subtract ionisation energy
      cosThetaSecondary= std::sqrt(deltaE*rb/(kineticEnergy*(deltaE+2.0*electron_mass_c2)));
      if (verboseLevel > 3)
    G4cout << "SampleFinalStatePositron: sampled close collision " << G4endl;
      return;                    
    }

  //Hard distant longitudinal collisions
  TS1 += XHDL;
  G4double deltaE = resEne;
  kineticEnergy1 = kineticEnergy - deltaE;
  if (TST < TS1)
    {
      G4double QS = QM/(1.0+QM*0.5/electron_mass_c2);
      G4double Q = QS/(std::pow((QS/cutoffEne)*(1.0+cutoffEne*0.5/electron_mass_c2),G4UniformRand())
               - (QS*0.5/electron_mass_c2));
      G4double QTREV = Q*(Q+2.0*electron_mass_c2);
      G4double cpps = kineticEnergy1*(kineticEnergy1+2.0*electron_mass_c2);
      cosThetaPrimary = (cpps+cps-QTREV)/(2.0*cp*std::sqrt(cpps));
      if (cosThetaPrimary > 1.) 
    cosThetaPrimary = 1.0;
      //energy and emission angle of the delta ray
      energySecondary = deltaE - ionEne;
      cosThetaSecondary = 0.5*(deltaE*(kineticEnergy+rb-deltaE)+QTREV)/std::sqrt(cps*QTREV);
      if (cosThetaSecondary > 1.0)
    cosThetaSecondary = 1.0;
      if (verboseLevel > 3)
    G4cout << "SampleFinalStatePositron: sampled distant longitudinal collision " << G4endl;
      return;      
    }

  //Hard distant transverse collisions
  cosThetaPrimary = 1.0;
  //energy and emission angle of the delta ray
  energySecondary = deltaE - ionEne;
  cosThetaSecondary = 0.5;

  if (verboseLevel > 3)    
    G4cout << "SampleFinalStatePositron: sampled distant transverse collision " << G4endl;
    
  return;
}

Here is the call graph for this function:

Here is the caller graph for this function:

SampleSecondaries()

void G4PenelopeIonisationModel::SampleSecondaries ( std::vector< G4DynamicParticle *> *  fvect, const G4MaterialCutsCouple *  couple, const G4DynamicParticle *  aDynamicParticle, G4double  tmin, G4double  maxEnergy  )

virtual

Implements G4VEmModel.

Definition at line 468 of file G4PenelopeIonisationModel.cc.

{
  // 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();
  const 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
  //Disable the built-in de-excitation of the PIXE flag is active. In this 
  //case, the PIXE interface takes care (statistically) of producing the 
  //de-excitation.
  //Now, take care of fluorescence, if required
  if (fAtomDeexcitation && !fPIXEflag && 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;      
    }
    
}

Here is the call graph for this function:

SetParticle()

void G4PenelopeIonisationModel::SetParticle ( const G4ParticleDefinition *  p )

private

Definition at line 1124 of file G4PenelopeIonisationModel.cc.

{
  if(!fParticle) {
    fParticle = p;  
  }
}

Here is the caller graph for this function:

SetVerbosityLevel()

void G4PenelopeIonisationModel::SetVerbosityLevel ( G4int  lev )

inline

Definition at line 112 of file G4PenelopeIonisationModel.hh.

{verboseLevel = lev;};

Member Data Documentation

cosThetaPrimary

G4double G4PenelopeIonisationModel::cosThetaPrimary

private

Definition at line 145 of file G4PenelopeIonisationModel.hh.

cosThetaSecondary

G4double G4PenelopeIonisationModel::cosThetaSecondary

private

Definition at line 147 of file G4PenelopeIonisationModel.hh.

energySecondary

G4double G4PenelopeIonisationModel::energySecondary

private

Definition at line 146 of file G4PenelopeIonisationModel.hh.

fAtomDeexcitation

G4VAtomDeexcitation* G4PenelopeIonisationModel::fAtomDeexcitation

private

Definition at line 140 of file G4PenelopeIonisationModel.hh.

fIntrinsicHighEnergyLimit

G4double G4PenelopeIonisationModel::fIntrinsicHighEnergyLimit

private

Definition at line 135 of file G4PenelopeIonisationModel.hh.

fIntrinsicLowEnergyLimit

G4double G4PenelopeIonisationModel::fIntrinsicLowEnergyLimit

private

Definition at line 134 of file G4PenelopeIonisationModel.hh.

fLocalTable

G4bool G4PenelopeIonisationModel::fLocalTable

private

Definition at line 156 of file G4PenelopeIonisationModel.hh.

fParticle

const G4ParticleDefinition* G4PenelopeIonisationModel::fParticle

protected

Definition at line 117 of file G4PenelopeIonisationModel.hh.

fParticleChange

G4ParticleChangeForLoss* G4PenelopeIonisationModel::fParticleChange

protected

Definition at line 113 of file G4PenelopeIonisationModel.hh.

fPIXEflag

G4bool G4PenelopeIonisationModel::fPIXEflag

private

Definition at line 141 of file G4PenelopeIonisationModel.hh.

isInitialised

G4bool G4PenelopeIonisationModel::isInitialised

private

Definition at line 139 of file G4PenelopeIonisationModel.hh.

kineticEnergy1

G4double G4PenelopeIonisationModel::kineticEnergy1

private

Definition at line 144 of file G4PenelopeIonisationModel.hh.

nBins

size_t G4PenelopeIonisationModel::nBins

private

Definition at line 153 of file G4PenelopeIonisationModel.hh.

oscManager

G4PenelopeOscillatorManager* G4PenelopeIonisationModel::oscManager

private

Definition at line 150 of file G4PenelopeIonisationModel.hh.

targetOscillator

G4int G4PenelopeIonisationModel::targetOscillator

private

Definition at line 148 of file G4PenelopeIonisationModel.hh.

theCrossSectionHandler

G4PenelopeIonisationXSHandler* G4PenelopeIonisationModel::theCrossSectionHandler

private

Definition at line 151 of file G4PenelopeIonisationModel.hh.

verboseLevel

G4int G4PenelopeIonisationModel::verboseLevel

private

Definition at line 137 of file G4PenelopeIonisationModel.hh.

---

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

• G4PenelopeIonisationModel.hh
• G4PenelopeIonisationModel.cc