G4LivermorePolarizedComptonModel Class Reference

#include <G4LivermorePolarizedComptonModel.hh>

Inheritance diagram for G4LivermorePolarizedComptonModel:

Public Member Functions
	G4LivermorePolarizedComptonModel (const G4ParticleDefinition *p=0, const G4String &nam="LivermorePolarizedCompton")
virtual	~G4LivermorePolarizedComptonModel ()
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)
Protected Attributes
G4ParticleChangeForGamma *	fParticleChange

---

Detailed Description

Definition at line 45 of file G4LivermorePolarizedComptonModel.hh.

---

Constructor & Destructor Documentation

G4LivermorePolarizedComptonModel::G4LivermorePolarizedComptonModel	(	const G4ParticleDefinition *	p = 0,
		const G4String &	nam = "LivermorePolarizedCompton"
	)

Definition at line 51 of file G4LivermorePolarizedComptonModel.cc.

References G4cout, G4endl, and G4VEmModel::SetHighEnergyLimit().

  :G4VEmModel(nam),fParticleChange(0),isInitialised(false),
   meanFreePathTable(0),scatterFunctionData(0),crossSectionHandler(0)
{
  lowEnergyLimit = 250 * eV; 
  highEnergyLimit = 100 * GeV;
  //SetLowEnergyLimit(lowEnergyLimit);
  SetHighEnergyLimit(highEnergyLimit);
  
  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

  if( verboseLevel>0 ) { 
  G4cout << "Livermore Polarized Compton is constructed " << G4endl
         << "Energy range: "
         << lowEnergyLimit / eV << " eV - "
         << highEnergyLimit / GeV << " GeV"
         << G4endl;
  }
}

G4LivermorePolarizedComptonModel::~G4LivermorePolarizedComptonModel

(

)

[virtual]

Definition at line 80 of file G4LivermorePolarizedComptonModel.cc.

{  
  if (meanFreePathTable)   delete meanFreePathTable;
  if (crossSectionHandler) delete crossSectionHandler;
  if (scatterFunctionData) delete scatterFunctionData;
}

---

Member Function Documentation

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

Reimplemented from G4VEmModel.

Definition at line 143 of file G4LivermorePolarizedComptonModel.cc.

References G4VCrossSectionHandler::FindValue(), G4cout, and G4endl.

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

  if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) return 0.0;

  G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
  return cs;
}

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

Implements G4VEmModel.

Definition at line 89 of file G4LivermorePolarizedComptonModel.cc.

References G4VCrossSectionHandler::BuildMeanFreePathForMaterials(), G4VCrossSectionHandler::Clear(), fParticleChange, G4cout, G4endl, G4VEmModel::GetParticleChangeForGamma(), G4VEmModel::HighEnergyLimit(), G4VEmModel::InitialiseElementSelectors(), G4ShellData::LoadData(), G4VEMDataSet::LoadData(), G4VCrossSectionHandler::LoadData(), G4VEmModel::LowEnergyLimit(), and G4ShellData::SetOccupancyData().

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

  if (crossSectionHandler)
  {
    crossSectionHandler->Clear();
    delete crossSectionHandler;
  }

  // Reading of data files - all materials are read
  
  crossSectionHandler = new G4CrossSectionHandler;
  crossSectionHandler->Clear();
  G4String crossSectionFile = "comp/ce-cs-";
  crossSectionHandler->LoadData(crossSectionFile);

  meanFreePathTable = 0;
  meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();

  G4VDataSetAlgorithm* scatterInterpolation = new G4LogLogInterpolation;
  G4String scatterFile = "comp/ce-sf-";
  scatterFunctionData = new G4CompositeEMDataSet(scatterInterpolation, 1., 1.);
  scatterFunctionData->LoadData(scatterFile);

  // For Doppler broadening
  shellData.SetOccupancyData();
  G4String file = "/doppler/shell-doppler";
  shellData.LoadData(file);

  if (verboseLevel > 2) 
    G4cout << "Loaded cross section files for Livermore Polarized Compton model" << G4endl;

  InitialiseElementSelectors(particle,cuts);

  if(  verboseLevel>0 ) { 
    G4cout << "Livermore Polarized Compton model is initialized " << G4endl
         << "Energy range: "
         << LowEnergyLimit() / eV << " eV - "
         << HighEnergyLimit() / GeV << " GeV"
         << G4endl;
  }
  
  //
    
  if(isInitialised) return;
  fParticleChange = GetParticleChangeForGamma();
  isInitialised = true;
}

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

Implements G4VEmModel.

Definition at line 160 of file G4LivermorePolarizedComptonModel.cc.

References G4ShellData::BindingEnergy(), G4Electron::Electron(), G4VEMDataSet::FindValue(), fParticleChange, fStopAndKill, G4cout, G4endl, G4UniformRand, G4DynamicParticle::GetDefinition(), G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMomentumDirection(), G4DynamicParticle::GetPolarization(), G4Element::GetZ(), G4VParticleChange::ProposeLocalEnergyDeposit(), G4ParticleChangeForGamma::ProposeMomentumDirection(), G4ParticleChangeForGamma::ProposePolarization(), G4VParticleChange::ProposeTrackStatus(), G4DopplerProfile::RandomSelectMomentum(), G4VEmModel::SelectRandomAtom(), G4ShellData::SelectRandomShell(), and G4ParticleChangeForGamma::SetProposedKineticEnergy().

{
  // The scattered gamma energy is sampled according to Klein - Nishina formula.
  // The random number techniques of Butcher & Messel are used (Nuc Phys 20(1960),15).
  // GEANT4 internal units
  //
  // Note : Effects due to binding of atomic electrons are negliged.

  if (verboseLevel > 3)
    G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedComptonModel" << G4endl;

  G4double gammaEnergy0 = aDynamicGamma->GetKineticEnergy();
  G4ThreeVector gammaPolarization0 = aDynamicGamma->GetPolarization();

  // Protection: a polarisation parallel to the
  // direction causes problems;
  // in that case find a random polarization

  G4ThreeVector gammaDirection0 = aDynamicGamma->GetMomentumDirection();

  // Make sure that the polarization vector is perpendicular to the
  // gamma direction. If not

  if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0))
    { // only for testing now
      gammaPolarization0 = GetRandomPolarization(gammaDirection0);
    }
  else
    {
      if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0)
        {
          gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0);
        }
    }

  // End of Protection

  // Within energy limit?

  if(gammaEnergy0 <= lowEnergyLimit)
    {
      fParticleChange->ProposeTrackStatus(fStopAndKill);
      fParticleChange->SetProposedKineticEnergy(0.);
      fParticleChange->ProposeLocalEnergyDeposit(gammaEnergy0);
      return;
    }

  G4double E0_m = gammaEnergy0 / electron_mass_c2 ;

  // Select randomly one element in the current material
  //G4int Z = crossSectionHandler->SelectRandomAtom(couple,gammaEnergy0);
  const G4ParticleDefinition* particle =  aDynamicGamma->GetDefinition();
  const G4Element* elm = SelectRandomAtom(couple,particle,gammaEnergy0);
  G4int Z = (G4int)elm->GetZ();

  // Sample the energy and the polarization of the scattered photon

  G4double epsilon, epsilonSq, onecost, sinThetaSqr, greject ;

  G4double epsilon0Local = 1./(1. + 2*E0_m);
  G4double epsilon0Sq = epsilon0Local*epsilon0Local;
  G4double alpha1   = - std::log(epsilon0Local);
  G4double alpha2 = 0.5*(1.- epsilon0Sq);

  G4double wlGamma = h_Planck*c_light/gammaEnergy0;
  G4double gammaEnergy1;
  G4ThreeVector gammaDirection1;

  do {
    if ( alpha1/(alpha1+alpha2) > G4UniformRand() )
      {
        epsilon   = std::exp(-alpha1*G4UniformRand());  
        epsilonSq = epsilon*epsilon; 
      }
    else 
      {
        epsilonSq = epsilon0Sq + (1.- epsilon0Sq)*G4UniformRand();
        epsilon   = std::sqrt(epsilonSq);
      }

    onecost = (1.- epsilon)/(epsilon*E0_m);
    sinThetaSqr   = onecost*(2.-onecost);

    // Protection
    if (sinThetaSqr > 1.)
      {
        G4cout
          << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
          << "sin(theta)**2 = "
          << sinThetaSqr
          << "; set to 1"
          << G4endl;
        sinThetaSqr = 1.;
      }
    if (sinThetaSqr < 0.)
      {
        G4cout
          << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
          << "sin(theta)**2 = "
          << sinThetaSqr
          << "; set to 0"
          << G4endl;
        sinThetaSqr = 0.;
      }
    // End protection

    G4double x =  std::sqrt(onecost/2.) / (wlGamma/cm);;
    G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1);
    greject = (1. - epsilon*sinThetaSqr/(1.+ epsilonSq))*scatteringFunction;

  } while(greject < G4UniformRand()*Z);


  // ****************************************************
  //            Phi determination
  // ****************************************************

  G4double phi = SetPhi(epsilon,sinThetaSqr);

  //
  // scattered gamma angles. ( Z - axis along the parent gamma)
  //

  G4double cosTheta = 1. - onecost;

  // Protection

  if (cosTheta > 1.)
    {
      G4cout
        << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
        << "cosTheta = "
        << cosTheta
        << "; set to 1"
        << G4endl;
      cosTheta = 1.;
    }
  if (cosTheta < -1.)
    {
      G4cout 
        << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
        << "cosTheta = " 
        << cosTheta
        << "; set to -1"
        << G4endl;
      cosTheta = -1.;
    }
  // End protection      
  
  
  G4double sinTheta = std::sqrt (sinThetaSqr);
  
  // Protection
  if (sinTheta > 1.)
    {
      G4cout 
        << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
        << "sinTheta = " 
        << sinTheta
        << "; set to 1"
        << G4endl;
      sinTheta = 1.;
    }
  if (sinTheta < -1.)
    {
      G4cout 
        << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
        << "sinTheta = " 
        << sinTheta
        << "; set to -1" 
        << G4endl;
      sinTheta = -1.;
    }
  // End protection
  
      
  G4double dirx = sinTheta*std::cos(phi);
  G4double diry = sinTheta*std::sin(phi);
  G4double dirz = cosTheta ;
  

  // oneCosT , eom

  // Doppler broadening -  Method based on:
  // Y. Namito, S. Ban and H. Hirayama, 
  // "Implementation of the Doppler Broadening of a Compton-Scattered Photon Into the EGS4 Code" 
  // NIM A 349, pp. 489-494, 1994
  
  // Maximum number of sampling iterations

  G4int maxDopplerIterations = 1000;
  G4double bindingE = 0.;
  G4double photonEoriginal = epsilon * gammaEnergy0;
  G4double photonE = -1.;
  G4int iteration = 0;
  G4double eMax = gammaEnergy0;

  do
    {
      iteration++;
      // Select shell based on shell occupancy
      G4int shell = shellData.SelectRandomShell(Z);
      bindingE = shellData.BindingEnergy(Z,shell);
      
      eMax = gammaEnergy0 - bindingE;
     
      // Randomly sample bound electron momentum (memento: the data set is in Atomic Units)
      G4double pSample = profileData.RandomSelectMomentum(Z,shell);
      // Rescale from atomic units
      G4double pDoppler = pSample * fine_structure_const;
      G4double pDoppler2 = pDoppler * pDoppler;
      G4double var2 = 1. + onecost * E0_m;
      G4double var3 = var2*var2 - pDoppler2;
      G4double var4 = var2 - pDoppler2 * cosTheta;
      G4double var = var4*var4 - var3 + pDoppler2 * var3;
      if (var > 0.)
        {
          G4double varSqrt = std::sqrt(var);        
          G4double scale = gammaEnergy0 / var3;  
          // Random select either root
          if (G4UniformRand() < 0.5) photonE = (var4 - varSqrt) * scale;               
          else photonE = (var4 + varSqrt) * scale;
        } 
      else
        {
          photonE = -1.;
        }
   } while ( iteration <= maxDopplerIterations && 
             (photonE < 0. || photonE > eMax || photonE < eMax*G4UniformRand()) );
 
  // End of recalculation of photon energy with Doppler broadening
  // Revert to original if maximum number of iterations threshold has been reached
  if (iteration >= maxDopplerIterations)
    {
      photonE = photonEoriginal;
      bindingE = 0.;
    }

  gammaEnergy1 = photonE;
 
  //
  // update G4VParticleChange for the scattered photon 
  //

  //  gammaEnergy1 = epsilon*gammaEnergy0;


  // New polarization

  G4ThreeVector gammaPolarization1 = SetNewPolarization(epsilon,
                                                        sinThetaSqr,
                                                        phi,
                                                        cosTheta);

  // Set new direction
  G4ThreeVector tmpDirection1( dirx,diry,dirz );
  gammaDirection1 = tmpDirection1;

  // Change reference frame.

  SystemOfRefChange(gammaDirection0,gammaDirection1,
                    gammaPolarization0,gammaPolarization1);

  if (gammaEnergy1 > 0.)
    {
      fParticleChange->SetProposedKineticEnergy( gammaEnergy1 ) ;
      fParticleChange->ProposeMomentumDirection( gammaDirection1 );
      fParticleChange->ProposePolarization( gammaPolarization1 );
    }
  else
    {
      gammaEnergy1 = 0.;
      fParticleChange->SetProposedKineticEnergy(0.) ;
      fParticleChange->ProposeTrackStatus(fStopAndKill);
    }

  //
  // kinematic of the scattered electron
  //

  G4double ElecKineEnergy = gammaEnergy0 - gammaEnergy1 -bindingE;

  // SI -protection against negative final energy: no e- is created
  // like in G4LivermoreComptonModel.cc
  if(ElecKineEnergy < 0.0) {
    fParticleChange->ProposeLocalEnergyDeposit(gammaEnergy0 - gammaEnergy1);
    return;
  }
 
  // SI - Removed range test
  
  G4double ElecMomentum = std::sqrt(ElecKineEnergy*(ElecKineEnergy+2.*electron_mass_c2));

  G4ThreeVector ElecDirection((gammaEnergy0 * gammaDirection0 -
                                   gammaEnergy1 * gammaDirection1) * (1./ElecMomentum));

  fParticleChange->ProposeLocalEnergyDeposit(bindingE);
  
  G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),ElecDirection.unit(),ElecKineEnergy) ;
  fvect->push_back(dp);

}

---

Field Documentation

G4ParticleChangeForGamma* G4LivermorePolarizedComptonModel::fParticleChange [protected]

Definition at line 73 of file G4LivermorePolarizedComptonModel.hh.

Referenced by Initialise(), and SampleSecondaries().

---

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

• G4LivermorePolarizedComptonModel.hh
• G4LivermorePolarizedComptonModel.cc

---