G4PenelopeRayleighModel Class Reference

#include <G4PenelopeRayleighModel.hh>

Inheritance diagram for G4PenelopeRayleighModel:

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

---

Detailed Description

Definition at line 57 of file G4PenelopeRayleighModel.hh.

---

Constructor & Destructor Documentation

G4PenelopeRayleighModel::G4PenelopeRayleighModel	(	const G4ParticleDefinition *	p = 0,
		const G4String &	processName = "PenRayleigh"
	)

Definition at line 53 of file G4PenelopeRayleighModel.cc.

References G4VEmModel::SetHighEnergyLimit().

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

  //build the energy grid. It is the same for all materials
  G4double logenergy = std::log(fIntrinsicLowEnergyLimit/2.);
  G4double logmaxenergy = std::log(1.5*fIntrinsicHighEnergyLimit);
  //finer grid below 160 keV
  G4double logtransitionenergy = std::log(160*keV); 
  G4double logfactor1 = std::log(10.)/250.;
  G4double logfactor2 = logfactor1*10;
  logEnergyGridPMax.push_back(logenergy);
  do{
    if (logenergy < logtransitionenergy)
      logenergy += logfactor1;
    else
      logenergy += logfactor2;
    logEnergyGridPMax.push_back(logenergy);      
  }while (logenergy < logmaxenergy);
}

G4PenelopeRayleighModel::~G4PenelopeRayleighModel

(

)

[virtual]

Definition at line 90 of file G4PenelopeRayleighModel.cc.

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

   if (atomicFormFactor)
     {
       for (i=atomicFormFactor->begin();i != atomicFormFactor->end();i++)
         if (i->second) delete i->second;
       delete atomicFormFactor;
     }

  ClearTables();
}

---

Member Function Documentation

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

Reimplemented from G4VEmModel.

Definition at line 186 of file G4PenelopeRayleighModel.cc.

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

{
  // Cross section of Rayleigh scattering in Penelope v2008 is calculated by the EPDL97 
  // tabulation, Cuellen et al. (1997), with non-relativistic form factors from Hubbel 
  // et al. J. Phys. Chem. Ref. Data 4 (1975) 471; Erratum ibid. 6 (1977) 615.

   if (verboseLevel > 3)
    G4cout << "Calling CrossSectionPerAtom() of G4PenelopeRayleighModel" << G4endl;
 
   G4int iZ = (G4int) Z;

   //read data files
   if (!logAtomicCrossSection->count(iZ))
     ReadDataFile(iZ);
   //now it should be ok
   if (!logAtomicCrossSection->count(iZ))
     {
       G4Exception("G4PenelopeRayleighModel::ComputeCrossSectionPerAtom()",
                   "em2040",FatalException,"Unable to load the cross section table");
     }

   G4double cross = 0;

   G4PhysicsFreeVector* atom = logAtomicCrossSection->find(iZ)->second;
   if (!atom)
     {
       G4ExceptionDescription ed;
       ed << "Unable to find Z=" << iZ << " in the atomic cross section table" << G4endl;
       G4Exception("G4PenelopeRayleighModel::ComputeCrossSectionPerAtom()",
                   "em2041",FatalException,ed);
       return 0;
     }
   G4double logene = std::log(energy);
   G4double logXS = atom->Value(logene);
   cross = std::exp(logXS);

   if (verboseLevel > 2)
    G4cout << "Rayleigh cross section at " << energy/keV << " keV for Z=" << Z << 
      " = " << cross/barn << " barn" << G4endl;
    return cross;
}

void G4PenelopeRayleighModel::DumpFormFactorTable

(

const G4Material *

)

Definition at line 1160 of file G4PenelopeRayleighModel.cc.

References G4cout, G4PhysicsVector::GetLowEdgeEnergy(), G4Material::GetName(), and G4PhysicsVector::GetVectorLength().

{
  G4cout << "*****************************************************************" << G4endl;
  G4cout << "G4PenelopeRayleighModel: Form Factor Table for " << mat->GetName() << G4endl;
  //try to use the same format as Penelope-Fortran, namely Q (/m_e*c) and F
  G4cout <<  "Q/(m_e*c)                 F(Q)     " << G4endl;
  G4cout << "*****************************************************************" << G4endl;
  if (!logFormFactorTable->count(mat))
    BuildFormFactorTable(mat);
  
  G4PhysicsFreeVector* theVec = logFormFactorTable->find(mat)->second;
  for (size_t i=0;i<theVec->GetVectorLength();i++)
    {
      G4double logQ2 = theVec->GetLowEdgeEnergy(i);
      G4double Q = std::exp(0.5*logQ2);
      G4double logF2 = (*theVec)[i];
      G4double F = std::exp(0.5*logF2);
      G4cout << Q << "              " << F << G4endl;
    }
  //DONE
  return;
}

G4int G4PenelopeRayleighModel::GetVerbosityLevel

(

)

[inline]

Definition at line 82 of file G4PenelopeRayleighModel.hh.

{return verboseLevel;};

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

Implements G4VEmModel.

Definition at line 145 of file G4PenelopeRayleighModel.cc.

References fParticleChange, G4cout, G4endl, G4VEmModel::GetParticleChangeForGamma(), G4VEmModel::HighEnergyLimit(), and G4VEmModel::LowEnergyLimit().

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

  //clear tables depending on materials, not the atomic ones
  ClearTables();
  
  //create new tables
  //
  // logAtomicCrossSection and atomicFormFactor are created only once,
  // since they are never cleared
  if (!logAtomicCrossSection)
    logAtomicCrossSection = new std::map<G4int,G4PhysicsFreeVector*>;
  if (!atomicFormFactor)
    atomicFormFactor = new std::map<G4int,G4PhysicsFreeVector*>;

  if (!logFormFactorTable)
    logFormFactorTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
  if (!pMaxTable)
    pMaxTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
  if (!samplingTable)
    samplingTable = new std::map<const G4Material*,G4PenelopeSamplingData*>;


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

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

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

Implements G4VEmModel.

Definition at line 306 of file G4PenelopeRayleighModel.cc.

References FatalException, fParticleChange, fStopAndKill, G4cout, G4endl, G4Exception(), G4UniformRand, G4DynamicParticle::GetKineticEnergy(), G4MaterialCutsCouple::GetMaterial(), G4DynamicParticle::GetMomentumDirection(), G4PenelopeSamplingData::GetNumberOfStoredPoints(), G4PenelopeSamplingData::GetX(), G4VParticleChange::ProposeLocalEnergyDeposit(), G4ParticleChangeForGamma::ProposeMomentumDirection(), G4VParticleChange::ProposeTrackStatus(), G4PenelopeSamplingData::SampleValue(), G4ParticleChangeForGamma::SetProposedKineticEnergy(), and G4PhysicsVector::Value().

{
  // Sampling of the Rayleigh final state (namely, scattering angle of the photon) 
  // from the Penelope2008 model. The scattering angle is sampled from the atomic 
  // cross section dOmega/d(cosTheta) from Born ("Atomic Phyisics", 1969), disregarding 
  // anomalous scattering effects. The Form Factor F(Q) function which appears in the 
  // analytical cross section is retrieved via the method GetFSquared(); atomic data 
  // are tabulated for F(Q). Form factor for compounds is calculated according to 
  // the additivity rule. The sampling from the F(Q) is made via a Rational Inverse 
  // Transform with Aliasing (RITA) algorithm; RITA parameters are calculated once 
  // for each material and managed by G4PenelopeSamplingData objects.
  // The sampling algorithm (rejection method) has efficiency 67% at low energy, and 
  // increases with energy. For E=100 keV the efficiency is 100% and 86% for 
  // hydrogen and uranium, respectively.

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

  G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
 
  if (photonEnergy0 <= fIntrinsicLowEnergyLimit)
    {
      fParticleChange->ProposeTrackStatus(fStopAndKill);
      fParticleChange->SetProposedKineticEnergy(0.);
      fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
      return ;
    }

  G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
  
  const G4Material* theMat = couple->GetMaterial();
  
  //1) Verify if tables are ready
  if (!pMaxTable || !samplingTable)
    {
      G4Exception("G4PenelopeRayleighModel::SampleSecondaries()",
                  "em2043",FatalException,"Invalid model initialization");    
      return;
    }
  
  //2) retrieve or build the sampling table
  if (!(samplingTable->count(theMat)))
    InitializeSamplingAlgorithm(theMat);
  G4PenelopeSamplingData* theDataTable = samplingTable->find(theMat)->second;
  
  //3) retrieve or build the pMax data
  if (!pMaxTable->count(theMat))
    GetPMaxTable(theMat);
  G4PhysicsFreeVector* thePMax = pMaxTable->find(theMat)->second;

  G4double cosTheta = 1.0;
  
  //OK, ready to go!
  G4double qmax = 2.0*photonEnergy0/electron_mass_c2; //this is non-dimensional now

  if (qmax < 1e-10) //very low momentum transfer
    {
      G4bool loopAgain=false;
      do
        {
          loopAgain = false;
          cosTheta = 1.0-2.0*G4UniformRand();
          G4double G = 0.5*(1+cosTheta*cosTheta);
          if (G4UniformRand()>G)
            loopAgain = true;
        }while(loopAgain);
    }
  else //larger momentum transfer
    {
      size_t nData = theDataTable->GetNumberOfStoredPoints();
      G4double LastQ2inTheTable = theDataTable->GetX(nData-1);
      G4double q2max = std::min(qmax*qmax,LastQ2inTheTable);

      G4bool loopAgain = false;
      G4double MaxPValue = thePMax->Value(photonEnergy0);
      G4double xx=0;
      
      //Sampling by rejection method. The rejection function is 
      //G = 0.5*(1+cos^2(theta))
      //
      do{
        loopAgain = false;
        G4double RandomMax = G4UniformRand()*MaxPValue;
        xx = theDataTable->SampleValue(RandomMax);
        //xx is a random value of q^2 in (0,q2max),sampled according to 
        //F(Q^2) via the RITA algorithm
        if (xx > q2max)
          loopAgain = true;
        cosTheta = 1.0-2.0*xx/q2max;
        G4double G = 0.5*(1+cosTheta*cosTheta);
        if (G4UniformRand()>G)
          loopAgain = true;
      }while(loopAgain);
    }
  
  G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
 
  // Scattered photon angles. ( Z - axis along the parent photon)
  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) ;
  fParticleChange->SetProposedKineticEnergy(photonEnergy0) ;
 
  return;
}

void G4PenelopeRayleighModel::SetVerbosityLevel

(

G4int

lev

)

[inline]

Definition at line 81 of file G4PenelopeRayleighModel.hh.

{verboseLevel = lev;};

---

Field Documentation

G4ParticleChangeForGamma* G4PenelopeRayleighModel::fParticleChange [protected]

Definition at line 88 of file G4PenelopeRayleighModel.hh.

Referenced by Initialise(), and SampleSecondaries().

---

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

• G4PenelopeRayleighModel.hh
• G4PenelopeRayleighModel.cc

---