G4InitXscPAI Class Reference

#include <G4InitXscPAI.hh>

Collaboration diagram for G4InitXscPAI:

Public Member Functions
	G4InitXscPAI (const G4MaterialCutsCouple *matCC)
virtual	~G4InitXscPAI ()
void	KillCloseIntervals ()
void	Normalisation ()
G4double	RutherfordIntegral (G4int intervalNumber, G4double limitLow, G4double limitHigh)
G4double	IntegralTerm (G4double omega)
G4double	ImPartDielectricConst (G4int intervalNumber, G4double energy)
G4double	RePartDielectricConst (G4double energy)
G4double	ModuleSqDielectricConst (G4int intervalNumber, G4double energy)
G4double	DifPAIxSection (G4double omega)
G4double	DifPAIdEdx (G4double omega)
G4double	PAIdNdxCherenkov (G4double omega)
G4double	PAIdNdxPlasmon (G4double omega)
void	IntegralPAIxSection (G4double bg2, G4double Tmax)
void	IntegralCherenkov (G4double bg2, G4double Tmax)
void	IntegralPlasmon (G4double bg2, G4double Tmax)
void	IntegralPAIdEdx (G4double bg2, G4double Tmax)
G4double	GetPhotonLambda (G4double omega)
G4double	GetStepEnergyLoss (G4double step)
G4double	GetStepCerenkovLoss (G4double step)
G4double	GetStepPlasmonLoss (G4double step)
G4int	GetIntervalNumber () const
G4int	GetBinPAI () const
G4double	GetNormalizationCof () const
G4double	GetMatSandiaMatrix (G4int i, G4int j) const
G4PhysicsLogVector *	GetPAIxscVector () const
G4PhysicsLogVector *	GetPAIdEdxVector () const
G4PhysicsLogVector *	GetPAIphotonVector () const
G4PhysicsLogVector *	GetPAIelectronVector () const
G4PhysicsLogVector *	GetChCosSqVector () const
G4PhysicsLogVector *	GetChWidthVector () const

Private Member Functions
G4InitXscPAI &	operator= (const G4InitXscPAI &right)
	G4InitXscPAI (const G4InitXscPAI &)

Private Attributes
G4int	fIntervalNumber
G4double	fNormalizationCof
G4int	fCurrentInterval
G4int	fIntervalTmax
G4double	fBetaGammaSq
G4double	fTmax
G4double	fDensity
G4double	fElectronDensity
G4OrderedTable *	fMatSandiaMatrix
G4SandiaTable *	fSandia
G4PhysicsLogVector *	fPAIxscVector
G4PhysicsLogVector *	fPAIdEdxVector
G4PhysicsLogVector *	fPAIphotonVector
G4PhysicsLogVector *	fPAIelectronVector
G4PhysicsLogVector *	fChCosSqVector
G4PhysicsLogVector *	fChWidthVector

Static Private Attributes
static const G4double	fDelta = 0.005
static const G4int	fPAIbin = 100
static const G4double	fSolidDensity = 0.05*g/cm3

Detailed Description

Definition at line 48 of file G4InitXscPAI.hh.

Constructor & Destructor Documentation

G4InitXscPAI() [1/2]

G4InitXscPAI::G4InitXscPAI

(

const G4MaterialCutsCouple *

matCC

)

Definition at line 70 of file G4InitXscPAI.cc.

  : fPAIxscVector(NULL),
    fPAIdEdxVector(NULL),
    fPAIphotonVector(NULL),
    fPAIelectronVector(NULL),
    fChCosSqVector(NULL),
    fChWidthVector(NULL)
{
  G4int i, j, matIndex;
 
  fDensity         = matCC->GetMaterial()->GetDensity();
  fElectronDensity = matCC->GetMaterial()->GetElectronDensity();
  matIndex         = matCC->GetMaterial()->GetIndex();

  fSandia          = new G4SandiaTable(matIndex);
  fIntervalNumber  = fSandia->GetMaxInterval()-1;

  fMatSandiaMatrix = new G4OrderedTable();
 
  for (i = 0; i < fIntervalNumber; i++)
  {
    fMatSandiaMatrix->push_back(new G4DataVector(5,0.));
  }             
  for (i = 0; i < fIntervalNumber; i++)
  {
    (*(*fMatSandiaMatrix)[i])[0] = fSandia->GetSandiaMatTable(i,0);

    for(j = 1; j < 5 ; j++)
    {
      (*(*fMatSandiaMatrix)[i])[j] = fSandia->GetSandiaMatTable(i,j)*fDensity;
    }     
  }
  KillCloseIntervals();
  Normalisation();
  fBetaGammaSq = fTmax = 0.0;
  fIntervalTmax = fCurrentInterval = 0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

~G4InitXscPAI()

G4InitXscPAI::~G4InitXscPAI ( )

virtual

Definition at line 115 of file G4InitXscPAI.cc.

{
  if(fPAIxscVector)      delete fPAIxscVector;  
  if(fPAIdEdxVector)     delete fPAIdEdxVector;  
  if(fPAIphotonVector)   delete fPAIphotonVector;  
  if(fPAIelectronVector) delete fPAIelectronVector;  
  if(fChCosSqVector)     delete fChCosSqVector;  
  if(fChWidthVector)     delete fChWidthVector;  
  delete fSandia;
  delete fMatSandiaMatrix;
}

G4InitXscPAI() [2/2]

G4InitXscPAI::G4InitXscPAI ( const G4InitXscPAI &  )

private

Member Function Documentation

DifPAIdEdx()

G4double G4InitXscPAI::DifPAIdEdx

(

G4double

omega

)

Definition at line 477 of file G4InitXscPAI.cc.

{
  G4double dEdx = omega*DifPAIxSection(omega);
  return dEdx;
}

Here is the call graph for this function:

Here is the caller graph for this function:

DifPAIxSection()

G4double G4InitXscPAI::DifPAIxSection

(

G4double

omega

)

Definition at line 414 of file G4InitXscPAI.cc.

{
  G4int i = fCurrentInterval;
  G4double  betaGammaSq = fBetaGammaSq;       
  G4double integralTerm = IntegralTerm(omega);
  G4double be2,cof,x1,x2,x3,x4,x5,x6,x7,x8,result ;
  G4double epsilonRe = RePartDielectricConst(omega);
  G4double epsilonIm = ImPartDielectricConst(i,omega);
  G4double be4 ;
  static const G4double betaBohr2 = fine_structure_const*fine_structure_const ;
  static const G4double betaBohr4 = betaBohr2*betaBohr2*4.0 ;
  be2 = betaGammaSq/(1 + betaGammaSq) ;
  be4 = be2*be2 ;
 
   cof = 1 ;
   x1 = log(2*electron_mass_c2/omega) ;

   if( betaGammaSq < 0.01 ) x2 = log(be2) ;
   else
   {
     x2 = -log( (1/betaGammaSq - epsilonRe)*
            (1/betaGammaSq - epsilonRe) + 
            epsilonIm*epsilonIm )/2 ;
   }
   if( epsilonIm == 0.0 || betaGammaSq < 0.01 )
   {
     x6=0 ;
   }
   else
   {
     x3 = -epsilonRe + 1/betaGammaSq ;
     x5 = -1 - epsilonRe + be2*((1 +epsilonRe)*(1 + epsilonRe) +
       epsilonIm*epsilonIm) ;

     x7 = atan2(epsilonIm,x3) ;
     x6 = x5 * x7 ;
   }
    // if(fImPartDielectricConst[i] == 0) x6 = 0 ;
   
   x4 = ((x1 + x2)*epsilonIm + x6)/hbarc ;
   //   if( x4 < 0.0 ) x4 = 0.0 ;
   x8 = (1 + epsilonRe)*(1 + epsilonRe) + 
        epsilonIm*epsilonIm;

   result = (x4 + cof*integralTerm/omega/omega) ;
   if(result < 1.0e-8) result = 1.0e-8 ;
   result *= fine_structure_const/be2/pi ;
   //   result *= (1-exp(-beta/betaBohr))*(1-exp(-beta/betaBohr)) ;
   //  result *= (1-exp(-be2/betaBohr2)) ;
   result *= (1-exp(-be4/betaBohr4)) ;
   if(fDensity >= fSolidDensity)
   { 
      result /= x8 ;
   }
   return result ;

} // end of DifPAIxSection 

Here is the call graph for this function:

Here is the caller graph for this function:

GetBinPAI()

G4int G4InitXscPAI::GetBinPAI ( ) const

inline

Definition at line 106 of file G4InitXscPAI.hh.

{ return fPAIbin ; }

GetChCosSqVector()

G4PhysicsLogVector* G4InitXscPAI::GetChCosSqVector ( ) const

inline

Definition at line 117 of file G4InitXscPAI.hh.

{ return fChCosSqVector;}

GetChWidthVector()

G4PhysicsLogVector* G4InitXscPAI::GetChWidthVector ( ) const

inline

Definition at line 118 of file G4InitXscPAI.hh.

{ return fChWidthVector;}

Here is the call graph for this function:

GetIntervalNumber()

G4int G4InitXscPAI::GetIntervalNumber ( ) const

inline

Definition at line 105 of file G4InitXscPAI.hh.

{ return fIntervalNumber ; }

GetMatSandiaMatrix()

G4double G4InitXscPAI::GetMatSandiaMatrix ( G4int  i, G4int  j  ) const

inline

Definition at line 110 of file G4InitXscPAI.hh.

          { return (*(*fMatSandiaMatrix)[i])[j]; }

GetNormalizationCof()

G4double G4InitXscPAI::GetNormalizationCof ( ) const

inline

Definition at line 108 of file G4InitXscPAI.hh.

{ return fNormalizationCof ; }

GetPAIdEdxVector()

G4PhysicsLogVector* G4InitXscPAI::GetPAIdEdxVector ( ) const

inline

Definition at line 114 of file G4InitXscPAI.hh.

{ return fPAIdEdxVector;}

GetPAIelectronVector()

G4PhysicsLogVector* G4InitXscPAI::GetPAIelectronVector ( ) const

inline

Definition at line 116 of file G4InitXscPAI.hh.

{ return fPAIelectronVector;}

GetPAIphotonVector()

G4PhysicsLogVector* G4InitXscPAI::GetPAIphotonVector ( ) const

inline

Definition at line 115 of file G4InitXscPAI.hh.

{ return fPAIphotonVector;}

GetPAIxscVector()

G4PhysicsLogVector* G4InitXscPAI::GetPAIxscVector ( ) const

inline

Definition at line 113 of file G4InitXscPAI.hh.

{ return fPAIxscVector;}

GetPhotonLambda()

G4double G4InitXscPAI::GetPhotonLambda

(

G4double

omega

)

Definition at line 933 of file G4InitXscPAI.cc.

{  
  G4int i ;
  G4double omega2, omega3, omega4, a1, a2, a3, a4, lambda ;

  omega2 = omega*omega ;
  omega3 = omega2*omega ;
  omega4 = omega2*omega2 ;

  for(i = 0; i < fIntervalNumber;i++)
  {
    if( omega < (*(*fMatSandiaMatrix)[i])[0] ) break ;
  }
  if( i == 0 )
  {
    G4cout<<"Warning: energy in G4InitXscPAI::GetPhotonLambda < I1"<<G4endl;
  }
  else i-- ;

  a1 = (*(*fMatSandiaMatrix)[i])[1]; 
  a2 = (*(*fMatSandiaMatrix)[i])[2]; 
  a3 = (*(*fMatSandiaMatrix)[i])[3]; 
  a4 = (*(*fMatSandiaMatrix)[i])[4]; 

  lambda = 1./(a1/omega + a2/omega2 + a3/omega3 + a4/omega4);

  return lambda ;
}

GetStepCerenkovLoss()

G4double G4InitXscPAI::GetStepCerenkovLoss

(

G4double

step

)

Definition at line 982 of file G4InitXscPAI.cc.

{  
  G4double loss = 0.0 ;
  loss *= step;

  return loss ;
}

GetStepEnergyLoss()

G4double G4InitXscPAI::GetStepEnergyLoss

(

G4double

step

)

Definition at line 970 of file G4InitXscPAI.cc.

{  
  G4double loss = 0.0 ;
  loss *= step;

  return loss ;
}

GetStepPlasmonLoss()

G4double G4InitXscPAI::GetStepPlasmonLoss

(

G4double

step

)

Definition at line 994 of file G4InitXscPAI.cc.

{  


  G4double loss = 0.0 ;
  loss *= step;
  return loss ;
}

ImPartDielectricConst()

G4double G4InitXscPAI::ImPartDielectricConst	(	G4int	intervalNumber,
		G4double	energy
	)

Definition at line 294 of file G4InitXscPAI.cc.

{
   G4double energy2,energy3,energy4,a1,a2,a3,a4,result;

   a1 = (*(*fMatSandiaMatrix)[k])[1]; 
   a2 = (*(*fMatSandiaMatrix)[k])[2]; 
   a3 = (*(*fMatSandiaMatrix)[k])[3]; 
   a4 = (*(*fMatSandiaMatrix)[k])[4]; 

   energy2 = energy1*energy1;
   energy3 = energy2*energy1;
   energy4 = energy3*energy1;
   
   result  = a1/energy1+a2/energy2+a3/energy3+a4/energy4 ;  
   result *= hbarc/energy1 ;
   
   return result ;

}  // end of ImPartDielectricConst 

Here is the caller graph for this function:

IntegralCherenkov()

void G4InitXscPAI::IntegralCherenkov	(	G4double	bg2,
		G4double	Tmax
	)

Definition at line 758 of file G4InitXscPAI.cc.

{
  G4int i,k,i1,i2;
  G4double energy1, energy2, beta2, module2, cos2, width, result = 0.;
 
  fBetaGammaSq = bg2;
  fTmax        = Tmax;
  beta2        = bg2/(1+bg2);

  if(fPAIphotonVector) delete fPAIphotonVector;  
  if(fChCosSqVector)   delete fChCosSqVector;  
  if(fChWidthVector)   delete fChWidthVector;  
  
  fPAIphotonVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
  fChCosSqVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
  fChWidthVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);

  fPAIphotonVector->PutValue(fPAIbin-1,result);
  fChCosSqVector->PutValue(fPAIbin-1,1.);
  fChWidthVector->PutValue(fPAIbin-1,1e-7);
                
  for( i = fIntervalNumber - 1; i >= 0; i-- )
  {
    if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
  }
  if (i < 0) i = 0; // Tmax should be more than 
                    // first ionisation potential
  fIntervalTmax = i;

  G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;

  for( k = fPAIbin - 2; k >= 0; k-- )
  {
    energy1 = fPAIphotonVector->GetLowEdgeEnergy(k);
    energy2 = fPAIphotonVector->GetLowEdgeEnergy(k+1);

    for( i = fIntervalTmax; i >= 0; i-- ) 
    {
      if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
    }
    if(i < 0) i = 0;
    i2 = i;

    for( i = fIntervalTmax; i >= 0; i-- ) 
    {
      if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
    }
    if(i < 0) i = 0;
    i1 = i;

    module2 = ModuleSqDielectricConst(i1,energy1);
    cos2    = RePartDielectricConst(energy1)/module2/beta2;      
    width   = ImPartDielectricConst(i1,energy1)/module2/beta2;
      
    fChCosSqVector->PutValue(k,cos2);
    fChWidthVector->PutValue(k,width);

    if( i1 == i2 )
    {
      fCurrentInterval = i1;
      result += integral.Legendre10(this,&G4InitXscPAI::PAIdNdxCherenkov,
                                    energy1,energy2);
      fPAIphotonVector->PutValue(k,result);

    }
    else
    {
      for( i = i2; i >= i1; i-- ) 
      {
        fCurrentInterval = i;

        if( i==i2 )        result += integral.Legendre10(this,
                           &G4InitXscPAI::PAIdNdxCherenkov,
                           (*(*fMatSandiaMatrix)[i])[0] ,energy2);

    else if( i == i1 ) result += integral.Legendre10(this,
                           &G4InitXscPAI::PAIdNdxCherenkov,energy1,
                           (*(*fMatSandiaMatrix)[i+1])[0]);

        else               result += integral.Legendre10(this,
                           &G4InitXscPAI::PAIdNdxCherenkov,
                       (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
      }
      fPAIphotonVector->PutValue(k,result);
    }
    // G4cout<<k<<"\t"<<result<<G4endl; 
  }
  return;
}   // end of IntegralCerenkov 

Here is the call graph for this function:

IntegralPAIdEdx()

void G4InitXscPAI::IntegralPAIdEdx	(	G4double	bg2,
		G4double	Tmax
	)

Definition at line 678 of file G4InitXscPAI.cc.

{
  G4int i,k,i1,i2;
  G4double energy1, energy2, result = 0.;
 
  fBetaGammaSq = bg2;
  fTmax        = Tmax;

  if(fPAIdEdxVector) delete fPAIdEdxVector;  
  
  fPAIdEdxVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
  fPAIdEdxVector->PutValue(fPAIbin-1,result);
                
  for( i = fIntervalNumber - 1; i >= 0; i-- )
  {
    if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
  }
  if (i < 0) i = 0; // Tmax should be more than 
                    // first ionisation potential
  fIntervalTmax = i;

  G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;

  for( k = fPAIbin - 2; k >= 0; k-- )
  {
    energy1 = fPAIdEdxVector->GetLowEdgeEnergy(k);
    energy2 = fPAIdEdxVector->GetLowEdgeEnergy(k+1);

    for( i = fIntervalTmax; i >= 0; i-- ) 
    {
      if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
    }
    if(i < 0) i = 0;
    i2 = i;

    for( i = fIntervalTmax; i >= 0; i-- ) 
    {
      if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
    }
    if(i < 0) i = 0;
    i1 = i;

    if( i1 == i2 )
    {
      fCurrentInterval = i1;
      result += integral.Legendre10(this,&G4InitXscPAI::DifPAIdEdx,
                                    energy1,energy2);
      fPAIdEdxVector->PutValue(k,result);
    }
    else
    {
      for( i = i2; i >= i1; i-- ) 
      {
        fCurrentInterval = i;

        if( i==i2 )        result += integral.Legendre10(this,
                           &G4InitXscPAI::DifPAIdEdx,
                           (*(*fMatSandiaMatrix)[i])[0] ,energy2);

    else if( i == i1 ) result += integral.Legendre10(this,
                           &G4InitXscPAI::DifPAIdEdx,energy1,
                           (*(*fMatSandiaMatrix)[i+1])[0]);

        else               result += integral.Legendre10(this,
                           &G4InitXscPAI::DifPAIdEdx,
                       (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
      }
      fPAIdEdxVector->PutValue(k,result);
    }
    // G4cout<<k<<"\t"<<result<<G4endl; 
  }
  return ;
}

Here is the call graph for this function:

IntegralPAIxSection()

void G4InitXscPAI::IntegralPAIxSection	(	G4double	bg2,
		G4double	Tmax
	)

Definition at line 597 of file G4InitXscPAI.cc.

{
  G4int i,k,i1,i2;
  G4double energy1, energy2, result = 0.;
 
  fBetaGammaSq = bg2;
  fTmax        = Tmax;

  if(fPAIxscVector) delete fPAIxscVector;  
  
  fPAIxscVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
  fPAIxscVector->PutValue(fPAIbin-1,result);
                
  for( i = fIntervalNumber - 1; i >= 0; i-- )
  {
    if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
  }
  if (i < 0) i = 0; // Tmax should be more than 
                    // first ionisation potential
  fIntervalTmax = i;

  G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;

  for( k = fPAIbin - 2; k >= 0; k-- )
  {
    energy1 = fPAIxscVector->GetLowEdgeEnergy(k);
    energy2 = fPAIxscVector->GetLowEdgeEnergy(k+1);

    for( i = fIntervalTmax; i >= 0; i-- ) 
    {
      if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
    }
    if(i < 0) i = 0;
    i2 = i;

    for( i = fIntervalTmax; i >= 0; i-- ) 
    {
      if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
    }
    if(i < 0) i = 0;
    i1 = i;

    if( i1 == i2 )
    {
      fCurrentInterval = i1;
      result += integral.Legendre10(this,&G4InitXscPAI::DifPAIxSection,
                                    energy1,energy2);
      fPAIxscVector->PutValue(k,result);
    }
    else
    {
      for( i = i2; i >= i1; i-- ) 
      {
        fCurrentInterval = i;

        if( i==i2 )        result += integral.Legendre10(this,
                           &G4InitXscPAI::DifPAIxSection,
                           (*(*fMatSandiaMatrix)[i])[0] ,energy2);

    else if( i == i1 ) result += integral.Legendre10(this,
                           &G4InitXscPAI::DifPAIxSection,energy1,
                           (*(*fMatSandiaMatrix)[i+1])[0]);

        else               result += integral.Legendre10(this,
                           &G4InitXscPAI::DifPAIxSection,
                       (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
      }
      fPAIxscVector->PutValue(k,result);
    }
    // G4cout<<k<<"\t"<<result<<G4endl; 
  }
  return ;
}

Here is the call graph for this function:

IntegralPlasmon()

void G4InitXscPAI::IntegralPlasmon	(	G4double	bg2,
		G4double	Tmax
	)

Definition at line 854 of file G4InitXscPAI.cc.

{
  G4int i,k,i1,i2;
  G4double energy1, energy2, result = 0.;
 
  fBetaGammaSq = bg2;
  fTmax        = Tmax;

  if(fPAIelectronVector) delete fPAIelectronVector;  
  
  fPAIelectronVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
  fPAIelectronVector->PutValue(fPAIbin-1,result);
                
  for( i = fIntervalNumber - 1; i >= 0; i-- )
  {
    if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
  }
  if (i < 0) i = 0; // Tmax should be more than 
                    // first ionisation potential
  fIntervalTmax = i;

  G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;

  for( k = fPAIbin - 2; k >= 0; k-- )
  {
    energy1 = fPAIelectronVector->GetLowEdgeEnergy(k);
    energy2 = fPAIelectronVector->GetLowEdgeEnergy(k+1);

    for( i = fIntervalTmax; i >= 0; i-- ) 
    {
      if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
    }
    if(i < 0) i = 0;
    i2 = i;

    for( i = fIntervalTmax; i >= 0; i-- ) 
    {
      if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
    }
    if(i < 0) i = 0;
    i1 = i;

    if( i1 == i2 )
    {
      fCurrentInterval = i1;
      result += integral.Legendre10(this,&G4InitXscPAI::PAIdNdxPlasmon,
                                    energy1,energy2);
      fPAIelectronVector->PutValue(k,result);
    }
    else
    {
      for( i = i2; i >= i1; i-- ) 
      {
        fCurrentInterval = i;

        if( i==i2 )        result += integral.Legendre10(this,
                           &G4InitXscPAI::PAIdNdxPlasmon,
                           (*(*fMatSandiaMatrix)[i])[0] ,energy2);

    else if( i == i1 ) result += integral.Legendre10(this,
                           &G4InitXscPAI::PAIdNdxPlasmon,energy1,
                           (*(*fMatSandiaMatrix)[i+1])[0]);

        else               result += integral.Legendre10(this,
                           &G4InitXscPAI::PAIdNdxPlasmon,
                       (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
      }
      fPAIelectronVector->PutValue(k,result);
    }
    // G4cout<<k<<"\t"<<result<<G4endl; 
  }
  return;
}   // end of IntegralPlasmon

Here is the call graph for this function:

IntegralTerm()

G4double G4InitXscPAI::IntegralTerm

(

G4double

omega

)

Definition at line 256 of file G4InitXscPAI.cc.

{
  G4int i;
  G4double energy1, energy2, result = 0.; 
                
  for( i = 0; i <= fIntervalTmax; i++ )
  {
    if(i == fIntervalTmax) 
    {
      energy1 = (*(*fMatSandiaMatrix)[i])[0];
      result += RutherfordIntegral(i,energy1,omega);
    }
    else 
    {
      if( omega <= (*(*fMatSandiaMatrix)[i+1])[0])
      {
        energy1 = (*(*fMatSandiaMatrix)[i])[0];
        result += RutherfordIntegral(i,energy1,omega);
        break;
      }
      else
      {
        energy1 = (*(*fMatSandiaMatrix)[i])[0];
        energy2 = (*(*fMatSandiaMatrix)[i+1])[0];
        result += RutherfordIntegral(i,energy1,energy2);
      }
    }
    // G4cout<<"IntegralTerm<<"("<<omega<<")"<<" = "<<result<<G4endl;
  }
  return result;
}

Here is the call graph for this function:

Here is the caller graph for this function:

KillCloseIntervals()

void G4InitXscPAI::KillCloseIntervals

(

)

Definition at line 131 of file G4InitXscPAI.cc.

{
  G4int i, j, k;
  G4double energy1, energy2; 
                
  for( i = 0 ; i < fIntervalNumber - 1 ; i++ )
  {
    energy1 = (*(*fMatSandiaMatrix)[i])[0];
    energy2 = (*(*fMatSandiaMatrix)[i+1])[0];

    if( energy2 - energy1 > 1.5*fDelta*(energy1 + energy2) )  continue ;
    else
    {
      for(j = i; j < fIntervalNumber-1; j++)
      {
        for( k = 0; k < 5; k++ )
        {
          (*(*fMatSandiaMatrix)[j])[k] = (*(*fMatSandiaMatrix)[j+1])[k];
    }
      }
      fIntervalNumber-- ;
      i-- ;
    }
  }
  
}

Here is the caller graph for this function:

ModuleSqDielectricConst()

G4double G4InitXscPAI::ModuleSqDielectricConst	(	G4int	intervalNumber,
		G4double	energy
	)

Definition at line 320 of file G4InitXscPAI.cc.

{
   G4double eIm2, eRe2, result;

   result = ImPartDielectricConst(k,omega);
   eIm2   = result*result;

   result = RePartDielectricConst(omega);
   eRe2   = result*result;

   result = eIm2 + eRe2;
  
   return result ;
}  

Here is the call graph for this function:

Here is the caller graph for this function:

Normalisation()

void G4InitXscPAI::Normalisation

(

)

Definition at line 162 of file G4InitXscPAI.cc.

{
  G4int i, j;
  G4double energy1, energy2, /*delta,*/ cof; // , shift;

  energy1 = (*(*fMatSandiaMatrix)[fIntervalNumber-1])[0];
  energy2 = 2.*(*(*fMatSandiaMatrix)[fIntervalNumber-1])[0];

 
  cof = RutherfordIntegral(fIntervalNumber-1,energy1,energy2);
                
  for( i = fIntervalNumber-2; i >= 0; i-- )
  {
    energy1 = (*(*fMatSandiaMatrix)[i])[0];
    energy2 = (*(*fMatSandiaMatrix)[i+1])[0];

    cof += RutherfordIntegral(i,energy1,energy2);
    // G4cout<<"norm. cof = "<<cof<<G4endl;
  }
  fNormalizationCof  = 2*pi*pi*hbarc*hbarc*fine_structure_const/electron_mass_c2 ;
  fNormalizationCof *= fElectronDensity;
  //delta = fNormalizationCof - cof;
  fNormalizationCof /= cof;
  //  G4cout<<"G4InitXscPAI::fNormalizationCof/cof = "<<fNormalizationCof
  //    <<";  at delta ="<<delta<<G4endl ;

  for (i = 0; i < fIntervalNumber; i++) // renormalisation on QM sum rule
  {
    for(j = 1; j < 5 ; j++)
    {
      (*(*fMatSandiaMatrix)[i])[j] *= fNormalizationCof;
    }     
  }
  /* 
  if(delta > 0) // shift the first energy interval
  {
    for(i=1;i<100;i++)
    {
      energy1 = (1.-i/100.)*(*(*fMatSandiaMatrix)[0])[0];
      energy2 = (*(*fMatSandiaMatrix)[0])[0];
      shift   = RutherfordIntegral(0,energy1,energy2);
      G4cout<<shift<<"\t";
      if(shift >= delta) break;
    }
    (*(*fMatSandiaMatrix)[0])[0] = energy1;
    cof += shift;
  }
  else if(delta < 0)
  {
    for(i=1;i<100;i++)
    {
      energy1 = (*(*fMatSandiaMatrix)[0])[0];
      energy2 = (*(*fMatSandiaMatrix)[0])[0] + 
              ( (*(*fMatSandiaMatrix)[0])[0] - (*(*fMatSandiaMatrix)[0])[0] )*i/100.;
      shift   = RutherfordIntegral(0,energy1,energy2);
      if( shift >= std::abs(delta) ) break;
    }
    (*(*fMatSandiaMatrix)[0])[0] = energy2;
    cof -= shift;
  }
  G4cout<<G4cout<<"G4InitXscPAI::fNormalizationCof/cof = "<<fNormalizationCof/cof
        <<";  at delta ="<<delta<<"  and i = "<<i<<G4endl ;
 */ 
}

Here is the call graph for this function:

Here is the caller graph for this function:

operator=()

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

private

Here is the caller graph for this function:

PAIdNdxCherenkov()

G4double G4InitXscPAI::PAIdNdxCherenkov

(

G4double

omega

)

Definition at line 487 of file G4InitXscPAI.cc.

{        
  G4int i = fCurrentInterval;
  G4double  betaGammaSq = fBetaGammaSq;       
  G4double epsilonRe = RePartDielectricConst(omega);
  G4double epsilonIm = ImPartDielectricConst(i,omega);

  G4double /*cof,*/ logarithm, x3, x5, argument, modul2, dNdxC ; 
  G4double be2, be4;

  //cof         = 1.0 ;
  static const G4double cofBetaBohr = 4.0 ;
  static const G4double betaBohr2   = fine_structure_const*fine_structure_const ;
  static const G4double betaBohr4   = betaBohr2*betaBohr2*cofBetaBohr ;

   be2 = betaGammaSq/(1 + betaGammaSq) ;
   be4 = be2*be2 ;

   if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq) ; // 0.0 ;
   else
   {
     logarithm  = -log( (1/betaGammaSq - epsilonRe)*
                    (1/betaGammaSq - epsilonRe) + 
                    epsilonIm*epsilonIm )*0.5 ;
     logarithm += log(1+1.0/betaGammaSq) ;
   }

   if( epsilonIm == 0.0 || betaGammaSq < 0.01 )
   {
     argument = 0.0 ;
   }
   else
   {
     x3 = -epsilonRe + 1.0/betaGammaSq ;
     x5 = -1.0 - epsilonRe +
          be2*((1.0 +epsilonRe)*(1.0 + epsilonRe) +
      epsilonIm*epsilonIm) ;
     if( x3 == 0.0 ) argument = 0.5*pi;
     else            argument = atan2(epsilonIm,x3) ;
     argument *= x5  ;
   }   
   dNdxC = ( logarithm*epsilonIm + argument )/hbarc ;
  
   if(dNdxC < 1.0e-8) dNdxC = 1.0e-8 ;

   dNdxC *= fine_structure_const/be2/pi ;

   dNdxC *= (1-exp(-be4/betaBohr4)) ;

   if(fDensity >= fSolidDensity)
   { 
      modul2 = (1.0 + epsilonRe)*(1.0 + epsilonRe) + 
                    epsilonIm*epsilonIm;
      dNdxC /= modul2 ;
   }
   return dNdxC ;

} // end of PAIdNdxCerenkov 

Here is the call graph for this function:

Here is the caller graph for this function:

PAIdNdxPlasmon()

G4double G4InitXscPAI::PAIdNdxPlasmon

(

G4double

omega

)

Definition at line 551 of file G4InitXscPAI.cc.

{        
  G4int i = fCurrentInterval;
  G4double  betaGammaSq = fBetaGammaSq;       
  G4double integralTerm = IntegralTerm(omega);
  G4double epsilonRe = RePartDielectricConst(omega);
  G4double epsilonIm = ImPartDielectricConst(i,omega);

   G4double cof, resonance, modul2, dNdxP ;
   G4double be2, be4;

   cof = 1 ;
   static const G4double cofBetaBohr = 4.0 ;
   static const G4double betaBohr2   = fine_structure_const*fine_structure_const ;
   static const G4double betaBohr4   = betaBohr2*betaBohr2*cofBetaBohr ;

   be2 = betaGammaSq/(1 + betaGammaSq) ;
   be4 = be2*be2 ;
 
   resonance  = log(2*electron_mass_c2*be2/omega) ;  
   resonance *= epsilonIm/hbarc ;


   dNdxP = ( resonance + cof*integralTerm/omega/omega ) ;

   if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8 ;

   dNdxP *= fine_structure_const/be2/pi ;
   dNdxP *= (1-exp(-be4/betaBohr4)) ;

   if( fDensity >= fSolidDensity )
   { 
     modul2 = (1 + epsilonRe)*(1 + epsilonRe) + 
               epsilonIm*epsilonIm;
     dNdxP /= modul2 ;
   }
   return dNdxP ;

} // end of PAIdNdxPlasmon 

Here is the call graph for this function:

Here is the caller graph for this function:

RePartDielectricConst()

G4double G4InitXscPAI::RePartDielectricConst

(

G4double

energy

)

Definition at line 343 of file G4InitXscPAI.cc.

{
  G4int i;       
   G4double x0, x02, x03, x04, x05, x1, x2, a1,a2,a3,a4,xx1 ,xx2 , xx12,
            c1, c2, c3, cof1, cof2, xln1, xln2, xln3, result ;

   x0 = enb ;
   result = 0 ;
   
   for( i = 0; i < fIntervalNumber-1; i++)
   {
      x1 = (*(*fMatSandiaMatrix)[i])[0];
      x2 = (*(*fMatSandiaMatrix)[i+1])[0] ;

      a1 = (*(*fMatSandiaMatrix)[i])[1]; 
      a2 = (*(*fMatSandiaMatrix)[i])[2]; 
      a3 = (*(*fMatSandiaMatrix)[i])[3]; 
      a4 = (*(*fMatSandiaMatrix)[i])[4];
 
      if( std::abs(x0-x1) < 0.5*(x0+x1)*fDelta ) 
      {
        if(x0 >= x1) x0 = x1*(1+fDelta);
        else         x0 = x1*(1-fDelta);
      } 
      if( std::abs(x0-x2) < 0.5*(x0+x2)*fDelta ) 
      {
        if(x0 >= x2) x0 = x2*(1+fDelta);
        else         x0 = x2*(1-fDelta);
      }
      xx1 = x1 - x0 ;
      xx2 = x2 - x0 ;
      xx12 = xx2/xx1 ;
      
      if( xx12 < 0 ) xx12 = -xx12;
      
      xln1 = log(x2/x1) ;
      xln2 = log(xx12) ;
      xln3 = log((x2 + x0)/(x1 + x0)) ;

      x02 = x0*x0 ;
      x03 = x02*x0 ;
      x04 = x03*x0 ;
      x05 = x04*x0;

      c1  = (x2 - x1)/x1/x2 ;
      c2  = (x2 - x1)*(x2 +x1)/x1/x1/x2/x2 ;
      c3  = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2 ;

      result -= (a1/x02 + a3/x04)*xln1 ;
      result -= (a2/x02 + a4/x04)*c1 ;
      result -= a3*c2/2/x02 ;
      result -= a4*c3/3/x02 ;

      cof1 = a1/x02 + a3/x04 ;
      cof2 = a2/x03 + a4/x05 ;

      result += 0.5*(cof1 +cof2)*xln2 ;
      result += 0.5*(cof1 - cof2)*xln3 ;
   } 
   result *= 2*hbarc/pi ;
   
   return result ;

}   // end of RePartDielectricConst 

Here is the caller graph for this function:

RutherfordIntegral()

G4double G4InitXscPAI::RutherfordIntegral	(	G4int	intervalNumber,
		G4double	limitLow,
		G4double	limitHigh
	)

Definition at line 232 of file G4InitXscPAI.cc.

{
   G4double  c1, c2, c3, a1, a2, a3, a4 ;

   a1 = (*(*fMatSandiaMatrix)[k])[1]; 
   a2 = (*(*fMatSandiaMatrix)[k])[2]; 
   a3 = (*(*fMatSandiaMatrix)[k])[3]; 
   a4 = (*(*fMatSandiaMatrix)[k])[4]; 
   // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;   
   c1 = (x2 - x1)/x1/x2 ;
   c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2 ;
   c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2 ;
   // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = "<<c2<<"; "<<"c3 = "<<c3<<G4endl;   
   
   return  a1*log(x2/x1) + a2*c1 + a3*c2/2 + a4*c3/3 ;

}   // end of RutherfordIntegral 

Here is the caller graph for this function:

Member Data Documentation

fBetaGammaSq

G4double G4InitXscPAI::fBetaGammaSq

private

Definition at line 137 of file G4InitXscPAI.hh.

fChCosSqVector

G4PhysicsLogVector* G4InitXscPAI::fChCosSqVector

private

Definition at line 153 of file G4InitXscPAI.hh.

fChWidthVector

G4PhysicsLogVector* G4InitXscPAI::fChWidthVector

private

Definition at line 154 of file G4InitXscPAI.hh.

fCurrentInterval

G4int G4InitXscPAI::fCurrentInterval

private

Definition at line 135 of file G4InitXscPAI.hh.

fDelta

const G4double G4InitXscPAI::fDelta = 0.005

staticprivate

Definition at line 129 of file G4InitXscPAI.hh.

fDensity

G4double G4InitXscPAI::fDensity

private

Definition at line 139 of file G4InitXscPAI.hh.

fElectronDensity

G4double G4InitXscPAI::fElectronDensity

private

Definition at line 140 of file G4InitXscPAI.hh.

fIntervalNumber

G4int G4InitXscPAI::fIntervalNumber

private

Definition at line 133 of file G4InitXscPAI.hh.

fIntervalTmax

G4int G4InitXscPAI::fIntervalTmax

private

Definition at line 136 of file G4InitXscPAI.hh.

fMatSandiaMatrix

G4OrderedTable* G4InitXscPAI::fMatSandiaMatrix

private

Definition at line 144 of file G4InitXscPAI.hh.

fNormalizationCof

G4double G4InitXscPAI::fNormalizationCof

private

Definition at line 134 of file G4InitXscPAI.hh.

fPAIbin

const G4int G4InitXscPAI::fPAIbin = 100

staticprivate

Definition at line 130 of file G4InitXscPAI.hh.

fPAIdEdxVector

G4PhysicsLogVector* G4InitXscPAI::fPAIdEdxVector

private

Definition at line 150 of file G4InitXscPAI.hh.

fPAIelectronVector

G4PhysicsLogVector* G4InitXscPAI::fPAIelectronVector

private

Definition at line 152 of file G4InitXscPAI.hh.

fPAIphotonVector

G4PhysicsLogVector* G4InitXscPAI::fPAIphotonVector

private

Definition at line 151 of file G4InitXscPAI.hh.

fPAIxscVector

G4PhysicsLogVector* G4InitXscPAI::fPAIxscVector

private

Definition at line 149 of file G4InitXscPAI.hh.

fSandia

G4SandiaTable* G4InitXscPAI::fSandia

private

Definition at line 145 of file G4InitXscPAI.hh.

fSolidDensity

const G4double G4InitXscPAI::fSolidDensity = 0.05*g/cm3

staticprivate

Definition at line 131 of file G4InitXscPAI.hh.

fTmax

G4double G4InitXscPAI::fTmax

private

Definition at line 138 of file G4InitXscPAI.hh.

---

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

• G4InitXscPAI.hh
• G4InitXscPAI.cc