#include <G4PenelopeComptonModel.hh>

Inheritance diagram for G4PenelopeComptonModel:

Public Member Functions
	G4PenelopeComptonModel (const G4ParticleDefinition *p=0, const G4String &processName="PenCompton")

virtual	~G4PenelopeComptonModel ()

virtual void	Initialise (const G4ParticleDefinition *, const G4DataVector &)

virtual void	InitialiseLocal (const G4ParticleDefinition , G4VEmModel masterModel)

virtual G4double	CrossSectionPerVolume (const G4Material , const G4ParticleDefinition , G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)

virtual G4double	ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double, G4double, G4double, G4double, G4double)

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

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	ComputeDEDXPerVolume (const G4Material , const G4ParticleDefinition , G4double kineticEnergy, G4double cutEnergy=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 G4double	MinEnergyCut (const G4ParticleDefinition , const G4MaterialCutsCouple )

virtual void	SetupForMaterial (const G4ParticleDefinition , const G4Material , G4double kineticEnergy)

virtual void	DefineForRegion (const G4Region *)

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

Protected Attributes
G4ParticleChangeForGamma *	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

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 *)

Detailed Description

Definition at line 63 of file G4PenelopeComptonModel.hh.

Constructor & Destructor Documentation

G4PenelopeComptonModel::G4PenelopeComptonModel	(	const G4ParticleDefinition *	p = `0`,
		const G4String &	processName = `"PenCompton"`
	)

Definition at line 62 of file G4PenelopeComptonModel.cc.

References python.hepunit::eV, G4PenelopeOscillatorManager::GetOscillatorManager(), python.hepunit::GeV, G4AtomicTransitionManager::Instance(), G4VEmModel::SetDeexcitationFlag(), and G4VEmModel::SetHighEnergyLimit().

   :G4VEmModel(nam),fParticleChange(0),fParticle(0),
    isInitialised(false),fAtomDeexcitation(0),
    oscManager(0)
 {
   fIntrinsicLowEnergyLimit = 100.0*eV;
   fIntrinsicHighEnergyLimit = 100.0*GeV;
   //  SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
   SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
   //
   oscManager = G4PenelopeOscillatorManager::GetOscillatorManager();
 
   if (part)
     SetParticle(part);
  
   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();
 }

G4PenelopeComptonModel::~G4PenelopeComptonModel ( )

virtual

Definition at line 94 of file G4PenelopeComptonModel.cc.

95 {;}

Member Function Documentation

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

virtual

Reimplemented from G4VEmModel.

Definition at line 234 of file G4PenelopeComptonModel.cc.

References G4cout, and G4endl.

 {
   G4cout << "*** G4PenelopeComptonModel -- WARNING ***" << G4endl;
   G4cout << "Penelope Compton 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 G4PenelopeComptonModel::CrossSectionPerVolume	(	const G4Material *	material,
		const G4ParticleDefinition *	p,
		G4double	kineticEnergy,
		G4double	cutEnergy = `0.0`,
		G4double	maxEnergy = `DBL_MAX`
	)

virtual

Reimplemented from G4VEmModel.

Definition at line 164 of file G4PenelopeComptonModel.cc.

References python.hepunit::classic_electr_radius, G4cout, G4endl, G4PenelopeOscillatorManager::GetAtomsPerMolecule(), G4Material::GetName(), G4PenelopeOscillatorManager::GetOscillatorTableCompton(), G4Material::GetTotNbOfAtomsPerVolume(), python.hepunit::keV, python.hepunit::MeV, python.hepunit::mm, python.hepunit::pi, and G4VEmModel::SetupForMaterial().

 {
   // Penelope model v2008 to calculate the Compton scattering cross section:
   // D. Brusa et al., Nucl. Instrum. Meth. A 379 (1996) 167
   // 
   // The cross section for Compton scattering is calculated according to the Klein-Nishina 
   // formula for energy > 5 MeV.
   // For E < 5 MeV it is used a parametrization for the differential cross-section dSigma/dOmega,
   // which is integrated numerically in cos(theta), G4PenelopeComptonModel::DifferentialCrossSection().
   // The parametrization includes the J(p) 
   // distribution profiles for the atomic shells, that are tabulated from Hartree-Fock calculations 
   // from F. Biggs et al., At. Data Nucl. Data Tables 16 (1975) 201
   //
   if (verboseLevel > 3)
     G4cout << "Calling CrossSectionPerVolume() of G4PenelopeComptonModel" << G4endl;
   SetupForMaterial(p, material, energy);
 
   //Retrieve the oscillator table for this material
   G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableCompton(material);
 
   G4double cs = 0;
 
   if (energy < 5*MeV) //explicit calculation for E < 5 MeV
     {
       size_t numberOfOscillators = theTable->size();
       for (size_t i=0;i<numberOfOscillators;i++)
     {
       G4PenelopeOscillator* theOsc = (*theTable)[i];
       //sum contributions coming from each oscillator
       cs += OscillatorTotalCrossSection(energy,theOsc);
     }
     }
   else //use Klein-Nishina for E>5 MeV
     cs = KleinNishinaCrossSection(energy,material);
        
   //cross sections are in units of pi*classic_electr_radius^2
   cs *= pi*classic_electr_radius*classic_electr_radius;
 
   //Now, cs is the cross section *per molecule*, let's calculate the 
   //cross section per volume
 
   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 csvolume = cs*moleculeDensity;
   
   if (verboseLevel > 2)
     G4cout << "Compton mean free path at " << energy/keV << " keV for material " << 
             material->GetName() << " = " << (1./csvolume)/mm << " mm" << G4endl;
   return csvolume;
 }

G4int G4PenelopeComptonModel::GetVerbosityLevel ( )

inline

Definition at line 100 of file G4PenelopeComptonModel.hh.

100 {return verboseLevel;};

void G4PenelopeComptonModel::Initialise	(	const G4ParticleDefinition *	part,
		const G4DataVector &
	)

virtual

Implements G4VEmModel.

Definition at line 99 of file G4PenelopeComptonModel.cc.

References G4LossTableManager::AtomDeexcitation(), fParticle, fParticleChange, G4cout, G4endl, G4VEmModel::GetParticleChangeForGamma(), python.hepunit::GeV, G4VEmModel::HighEnergyLimit(), G4LossTableManager::Instance(), G4VEmModel::IsMaster(), python.hepunit::keV, and G4VEmModel::LowEnergyLimit().

 {
   if (verboseLevel > 3)
     G4cout << "Calling G4PenelopeComptonModel::Initialise()" << G4endl;
 
   fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
   //Issue warning if the AtomicDeexcitation has not been declared
   if (!fAtomDeexcitation)
     {
       G4cout << G4endl;
       G4cout << "WARNING from G4PenelopeComptonModel " << 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;
     }
 
   SetParticle(part);
 
   if (IsMaster() && part == fParticle) 
     {
 
       if (verboseLevel > 0) 
     {
       G4cout << "Penelope Compton model v2008 is initialized " << G4endl
          << "Energy range: "
          << LowEnergyLimit() / keV << " keV - "
          << HighEnergyLimit() / GeV << " GeV";  
     }
     }      
 
   if(isInitialised) return;
   fParticleChange = GetParticleChangeForGamma();
   isInitialised = true; 
 
 }

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

virtual

Reimplemented from G4VEmModel.

Definition at line 138 of file G4PenelopeComptonModel.cc.

References fParticle, G4cout, and G4endl.

 {
   if (verboseLevel > 3)
     G4cout << "Calling  G4PenelopeComptonModel::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 G4PenelopeComptonModel* theModel = 
         static_cast<G4PenelopeComptonModel*> (masterModel);
       
       //Same verbosity for all workers, as the master
       verboseLevel = theModel->verboseLevel;
     }
 
   return;
 }

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

virtual

Implements G4VEmModel.

Definition at line 250 of file G4PenelopeComptonModel.cc.

References G4InuclSpecialFunctions::bindingEnergy(), G4AtomicShell::BindingEnergy(), G4VAtomDeexcitation::CheckDeexcitationActiveRegion(), G4Electron::Definition(), G4Gamma::Definition(), G4InuclParticleNames::electron, G4Electron::Electron(), python.hepunit::electron_mass_c2, python.hepunit::eV, fParticleChange, fStopAndKill, G4cout, G4endl, G4UniformRand, G4VAtomDeexcitation::GenerateParticles(), G4MaterialCutsCouple::GetIndex(), G4DynamicParticle::GetKineticEnergy(), G4MaterialCutsCouple::GetMaterial(), G4DynamicParticle::GetMomentumDirection(), G4PenelopeOscillatorManager::GetOscillatorTableCompton(), G4PenelopeOscillatorManager::GetTotalZ(), python.hepunit::keV, eplot::material, G4INCL::Math::max(), python.hepunit::MeV, G4INCL::Math::min(), nmax, python.hepunit::pi, G4VParticleChange::ProposeLocalEnergyDeposit(), G4ParticleChangeForGamma::ProposeMomentumDirection(), G4VParticleChange::ProposeTrackStatus(), CLHEP::Hep3Vector::rotateUz(), G4InuclParticleNames::s0, G4ParticleChangeForGamma::SetProposedKineticEnergy(), G4AtomicTransitionManager::Shell(), and python.hepunit::twopi.

 {
   
   // Penelope model v2008 to sample the Compton scattering final state.
   // D. Brusa et al., Nucl. Instrum. Meth. A 379 (1996) 167
   // The model determines also the original shell from which the electron is expelled, 
   // in order to produce fluorescence de-excitation (from G4DeexcitationManager)
   // 
   // The final state for Compton scattering is calculated according to the Klein-Nishina 
   // formula for energy > 5 MeV. In this case, the Doppler broadening is negligible and 
   // one can assume that the target electron is at rest.
   // For E < 5 MeV it is used the parametrization for the differential cross-section dSigma/dOmega,
   // to sample the scattering angle and the energy of the emerging electron, which is  
   // G4PenelopeComptonModel::DifferentialCrossSection(). The rejection method is 
   // used to sample cos(theta). The efficiency increases monotonically with photon energy and is 
   // nearly independent on the Z; typical values are 35%, 80% and 95% for 1 keV, 1 MeV and 10 MeV, 
   // respectively.
   // The parametrization includes the J(p) distribution profiles for the atomic shells, that are 
   // tabulated 
   // from Hartree-Fock calculations from F. Biggs et al., At. Data Nucl. Data Tables 16 (1975) 201. 
   // Doppler broadening is included.
   //
 
   if (verboseLevel > 3)
     G4cout << "Calling SampleSecondaries() of G4PenelopeComptonModel" << G4endl;
   
   G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
 
   if (photonEnergy0 <= fIntrinsicLowEnergyLimit)
     {
       fParticleChange->ProposeTrackStatus(fStopAndKill);
       fParticleChange->SetProposedKineticEnergy(0.);
       fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
       return ;
     }
 
   G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
   const G4Material* material = couple->GetMaterial();
 
   G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableCompton(material); 
 
   const G4int nmax = 64;
   G4double rn[nmax]={0.0};
   G4double pac[nmax]={0.0};
   
   G4double S=0.0;
   G4double epsilon = 0.0;
   G4double cosTheta = 1.0;
   G4double hartreeFunc = 0.0;
   G4double oscStren = 0.0;
   size_t numberOfOscillators = theTable->size();
   size_t targetOscillator = 0;
   G4double ionEnergy = 0.0*eV;
 
   G4double ek = photonEnergy0/electron_mass_c2;
   G4double ek2 = 2.*ek+1.0;
   G4double eks = ek*ek;
   G4double ek1 = eks-ek2-1.0;
 
   G4double taumin = 1.0/ek2;
   G4double a1 = std::log(ek2);
   G4double a2 = a1+2.0*ek*(1.0+ek)/(ek2*ek2);
 
   G4double TST = 0;
   G4double tau = 0.;
  
   //If the incoming photon is above 5 MeV, the quicker approach based on the 
   //pure Klein-Nishina formula is used
   if (photonEnergy0 > 5*MeV)
     {
       do{
     do{
       if ((a2*G4UniformRand()) < a1)
         tau = std::pow(taumin,G4UniformRand());     
       else
         tau = std::sqrt(1.0+G4UniformRand()*(taumin*taumin-1.0));       
       //rejection function
       TST = (1.0+tau*(ek1+tau*(ek2+tau*eks)))/(eks*tau*(1.0+tau*tau));
     }while (G4UniformRand()> TST);
     epsilon=tau;
     cosTheta = 1.0 - (1.0-tau)/(ek*tau);
 
     //Target shell electrons    
     TST = oscManager->GetTotalZ(material)*G4UniformRand();
     targetOscillator = numberOfOscillators-1; //last level
     S=0.0;
     G4bool levelFound = false;
     for (size_t j=0;j<numberOfOscillators && !levelFound; j++)
       {
         S += (*theTable)[j]->GetOscillatorStrength();       
         if (S > TST) 
           {
         targetOscillator = j;
         levelFound = true;
           }
       }
     //check whether the level is valid
     ionEnergy = (*theTable)[targetOscillator]->GetIonisationEnergy();
       }while((epsilon*photonEnergy0-photonEnergy0+ionEnergy) >0);
     }
   else //photonEnergy0 < 5 MeV
     {
       //Incoherent scattering function for theta=PI
       G4double s0=0.0;
       G4double pzomc=0.0;
       G4double rni=0.0;
       G4double aux=0.0;
       for (size_t i=0;i<numberOfOscillators;i++)
     {
       ionEnergy = (*theTable)[i]->GetIonisationEnergy();
       if (photonEnergy0 > ionEnergy)
         {
           G4double aux2 = photonEnergy0*(photonEnergy0-ionEnergy)*2.0;
           hartreeFunc = (*theTable)[i]->GetHartreeFactor(); 
           oscStren = (*theTable)[i]->GetOscillatorStrength();
           pzomc = hartreeFunc*(aux2-electron_mass_c2*ionEnergy)/
         (electron_mass_c2*std::sqrt(2.0*aux2+ionEnergy*ionEnergy));
           if (pzomc > 0)    
         rni = 1.0-0.5*std::exp(0.5-(std::sqrt(0.5)+std::sqrt(2.0)*pzomc)*
                        (std::sqrt(0.5)+std::sqrt(2.0)*pzomc));      
           else        
         rni = 0.5*std::exp(0.5-(std::sqrt(0.5)-std::sqrt(2.0)*pzomc)*
                    (std::sqrt(0.5)-std::sqrt(2.0)*pzomc));      
           s0 += oscStren*rni;
         }
     }      
       //Sampling tau
       G4double cdt1 = 0.;
       do
     {
       if ((G4UniformRand()*a2) < a1)        
         tau = std::pow(taumin,G4UniformRand());     
       else      
         tau = std::sqrt(1.0+G4UniformRand()*(taumin*taumin-1.0));       
       cdt1 = (1.0-tau)/(ek*tau);
       //Incoherent scattering function
       S = 0.;
       for (size_t i=0;i<numberOfOscillators;i++)
         {
           ionEnergy = (*theTable)[i]->GetIonisationEnergy();
           if (photonEnergy0 > ionEnergy) //sum only on excitable levels
         {
           aux = photonEnergy0*(photonEnergy0-ionEnergy)*cdt1;
           hartreeFunc = (*theTable)[i]->GetHartreeFactor(); 
           oscStren = (*theTable)[i]->GetOscillatorStrength();
           pzomc = hartreeFunc*(aux-electron_mass_c2*ionEnergy)/
             (electron_mass_c2*std::sqrt(2.0*aux+ionEnergy*ionEnergy));
           if (pzomc > 0) 
             rn[i] = 1.0-0.5*std::exp(0.5-(std::sqrt(0.5)+std::sqrt(2.0)*pzomc)*
                          (std::sqrt(0.5)+std::sqrt(2.0)*pzomc));            
           else          
             rn[i] = 0.5*std::exp(0.5-(std::sqrt(0.5)-std::sqrt(2.0)*pzomc)*
                      (std::sqrt(0.5)-std::sqrt(2.0)*pzomc));            
           S += oscStren*rn[i];
           pac[i] = S;
         }
           else
         pac[i] = S-1e-6;        
         }
       //Rejection function
       TST = S*(1.0+tau*(ek1+tau*(ek2+tau*eks)))/(eks*tau*(1.0+tau*tau));  
     }while ((G4UniformRand()*s0) > TST);
 
       cosTheta = 1.0 - cdt1;
       G4double fpzmax=0.0,fpz=0.0;
       G4double A=0.0;
       //Target electron shell
       do
     {
       do
         {
           TST = S*G4UniformRand();
           targetOscillator = numberOfOscillators-1; //last level
           G4bool levelFound = false;
           for (size_t i=0;i<numberOfOscillators && !levelFound;i++)
         {
           if (pac[i]>TST) 
             {            
               targetOscillator = i;
               levelFound = true;
             }
         }
           A = G4UniformRand()*rn[targetOscillator];
           hartreeFunc = (*theTable)[targetOscillator]->GetHartreeFactor(); 
           oscStren = (*theTable)[targetOscillator]->GetOscillatorStrength();
           if (A < 0.5) 
         pzomc = (std::sqrt(0.5)-std::sqrt(0.5-std::log(2.0*A)))/
           (std::sqrt(2.0)*hartreeFunc);       
           else      
         pzomc = (std::sqrt(0.5-std::log(2.0-2.0*A))-std::sqrt(0.5))/
           (std::sqrt(2.0)*hartreeFunc); 
         } while (pzomc < -1);
 
       // F(EP) rejection
       G4double XQC = 1.0+tau*(tau-2.0*cosTheta);
       G4double AF = std::sqrt(XQC)*(1.0+tau*(tau-cosTheta)/XQC);
       if (AF > 0) 
         fpzmax = 1.0+AF*0.2;
       else
         fpzmax = 1.0-AF*0.2;        
       fpz = 1.0+AF*std::max(std::min(pzomc,0.2),-0.2);
     }while ((fpzmax*G4UniformRand())>fpz);
   
       //Energy of the scattered photon
       G4double T = pzomc*pzomc;
       G4double b1 = 1.0-T*tau*tau;
       G4double b2 = 1.0-T*tau*cosTheta;
       if (pzomc > 0.0)  
     epsilon = (tau/b1)*(b2+std::sqrt(std::abs(b2*b2-b1*(1.0-T))));  
       else  
     epsilon = (tau/b1)*(b2-std::sqrt(std::abs(b2*b2-b1*(1.0-T))));  
     } //energy < 5 MeV
   
   //Ok, the kinematics has been calculated.
   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 ;
 
   // Update G4VParticleChange for the scattered photon
   G4ThreeVector photonDirection1(dirx,diry,dirz);
   photonDirection1.rotateUz(photonDirection0);
   fParticleChange->ProposeMomentumDirection(photonDirection1) ;
 
   G4double photonEnergy1 = epsilon * photonEnergy0;
 
   if (photonEnergy1 > 0.)  
     fParticleChange->SetProposedKineticEnergy(photonEnergy1) ;  
   else
   {
     fParticleChange->SetProposedKineticEnergy(0.) ;
     fParticleChange->ProposeTrackStatus(fStopAndKill);
   }
   
   //Create scattered electron
   G4double diffEnergy = photonEnergy0*(1-epsilon);
   ionEnergy = (*theTable)[targetOscillator]->GetIonisationEnergy();
 
   G4double Q2 = 
     photonEnergy0*photonEnergy0+photonEnergy1*(photonEnergy1-2.0*photonEnergy0*cosTheta);
   G4double cosThetaE = 0.; //scattering angle for the electron
 
   if (Q2 > 1.0e-12)    
     cosThetaE = (photonEnergy0-photonEnergy1*cosTheta)/std::sqrt(Q2);    
   else    
     cosThetaE = 1.0;    
   G4double sinThetaE = std::sqrt(1-cosThetaE*cosThetaE);
 
   //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-
   G4double bindingEnergy = 0.*eV;
   const G4AtomicShell* shell = 0;
 
   //Real level
   if (Z > 0 && shFlag<30)
     {
       shell = fTransitionManager->Shell(Z,shFlag-1);
       bindingEnergy = shell->BindingEnergy();    
     }
 
   G4double ionEnergyInPenelopeDatabase = ionEnergy;
   //protection against energy non-conservation
   ionEnergy = std::max(bindingEnergy,ionEnergyInPenelopeDatabase);  
 
   //subtract the excitation energy. If not emitted by fluorescence
   //the ionization energy is deposited as local energy deposition
   G4double eKineticEnergy = diffEnergy - ionEnergy; 
   G4double localEnergyDeposit = ionEnergy; 
   G4double energyInFluorescence = 0.; //testing purposes only
   G4double energyInAuger = 0; //testing purposes
 
   if (eKineticEnergy < 0) 
     {
       //It means that there was some problem/mismatch between the two databases. 
       //Try to make it work
       //In this case available Energy (diffEnergy) < ionEnergy
       //Full residual energy is deposited locally
       localEnergyDeposit = diffEnergy;
       eKineticEnergy = 0.0;
     }
 
   //the local energy deposit is what remains: part of this may be spent for fluorescence.
   //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;
           
         }
         }
 
     }
     }
 
 
   /*
   if(DeexcitationFlag() && Z > 5 && shellId>0) {
 
     const G4ProductionCutsTable* theCoupleTable=
       G4ProductionCutsTable::GetProductionCutsTable();
 
     size_t index = couple->GetIndex();
     G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[index];
     G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[index];
 
     // Generation of fluorescence
     // Data in EADL are available only for Z > 5
     // Protection to avoid generating photons in the unphysical case of
     // shell binding energy > photon energy
     if (localEnergyDeposit > cutg || localEnergyDeposit > cute)
       { 
     G4DynamicParticle* aPhoton;
     deexcitationManager.SetCutForSecondaryPhotons(cutg);
     deexcitationManager.SetCutForAugerElectrons(cute);
 
     photonVector = deexcitationManager.GenerateParticles(Z,shellId);
     if(photonVector) 
       {
         size_t nPhotons = photonVector->size();
         for (size_t k=0; k<nPhotons; k++)
           {
         aPhoton = (*photonVector)[k];
         if (aPhoton)
           {
             G4double itsEnergy = aPhoton->GetKineticEnergy();
             G4bool keepIt = false;
             if (itsEnergy <= localEnergyDeposit)
               {
             //check if good! 
             if(aPhoton->GetDefinition() == G4Gamma::Gamma()
                && itsEnergy >= cutg)
               {
                 keepIt = true;
                 energyInFluorescence += itsEnergy;            
               }
             if (aPhoton->GetDefinition() == G4Electron::Electron() && 
                 itsEnergy >= cute)
               {
                 energyInAuger += itsEnergy;
                 keepIt = true;
               }
               }
             //good secondary, register it
             if (keepIt)
               {
             localEnergyDeposit -= itsEnergy;
             fvect->push_back(aPhoton);
               }         
             else
               {
             delete aPhoton;
             (*photonVector)[k] = 0;
               }
           }
           }
         delete photonVector;
       }
       }
   }
   */
 
 
   //Always produce explicitely the electron 
   G4DynamicParticle* electron = 0;
 
   G4double xEl = sinThetaE * std::cos(phi+pi); 
   G4double yEl = sinThetaE * std::sin(phi+pi);
   G4double zEl = cosThetaE;
   G4ThreeVector eDirection(xEl,yEl,zEl); //electron direction
   eDirection.rotateUz(photonDirection0);
   electron = new G4DynamicParticle (G4Electron::Electron(),
                     eDirection,eKineticEnergy) ;
   fvect->push_back(electron);
     
 
   if (localEnergyDeposit < 0)
     {
       G4cout << "WARNING-" 
          << "G4PenelopeComptonModel::SampleSecondaries - Negative energy deposit"
          << G4endl;
       localEnergyDeposit=0.;
     }
   fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
   
   G4double electronEnergy = 0.;
   if (electron)
     electronEnergy = eKineticEnergy;
   if (verboseLevel > 1)
     {
       G4cout << "-----------------------------------------------------------" << G4endl;
       G4cout << "Energy balance from G4PenelopeCompton" << G4endl;
       G4cout << "Incoming photon energy: " << photonEnergy0/keV << " keV" << G4endl;
       G4cout << "-----------------------------------------------------------" << G4endl;
       G4cout << "Scattered photon: " << photonEnergy1/keV << " keV" << G4endl;
       G4cout << "Scattered electron " << electronEnergy/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: " << (photonEnergy1+electronEnergy+energyInFluorescence+
                       localEnergyDeposit+energyInAuger)/keV << 
     " keV" << G4endl;
       G4cout << "-----------------------------------------------------------" << G4endl;
     }
   if (verboseLevel > 0)
     {
       G4double energyDiff = std::fabs(photonEnergy1+
                       electronEnergy+energyInFluorescence+
                       localEnergyDeposit+energyInAuger-photonEnergy0);
       if (energyDiff > 0.05*keV)
     G4cout << "Warning from G4PenelopeCompton: problem with energy conservation: " << 
       (photonEnergy1+electronEnergy+energyInFluorescence+energyInAuger+localEnergyDeposit)/keV << 
       " keV (final) vs. " << 
       photonEnergy0/keV << " keV (initial)" << G4endl;
     }    
 }

void G4PenelopeComptonModel::SetVerbosityLevel ( G4int lev )

inline

Definition at line 99 of file G4PenelopeComptonModel.hh.

99 {verboseLevel = lev;};

Field Documentation

const G4ParticleDefinition* G4PenelopeComptonModel::fParticle

protected

Definition at line 105 of file G4PenelopeComptonModel.hh.

Referenced by Initialise(), and InitialiseLocal().

G4ParticleChangeForGamma* G4PenelopeComptonModel::fParticleChange

protected

Definition at line 100 of file G4PenelopeComptonModel.hh.

Referenced by Initialise(), and SampleSecondaries().

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

Public Member Functions

Protected Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Field Documentation