G4InitXscPAI Class Reference

#include <G4InitXscPAI.hh>

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

Detailed Description

Definition at line 48 of file G4InitXscPAI.hh.

Constructor & Destructor Documentation

G4InitXscPAI::G4InitXscPAI

(

const G4MaterialCutsCouple *

matCC

)

Definition at line 70 of file G4InitXscPAI.cc.

References G4Material::GetDensity(), G4Material::GetElectronDensity(), G4Material::GetIndex(), G4MaterialCutsCouple::GetMaterial(), G4SandiaTable::GetMaxInterval(), G4SandiaTable::GetSandiaMatTable(), KillCloseIntervals(), and Normalisation().

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

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

Member Function Documentation

G4double G4InitXscPAI::DifPAIdEdx

(

G4double

omega

)

Definition at line 477 of file G4InitXscPAI.cc.

References DifPAIxSection().

Referenced by IntegralPAIdEdx().

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

G4double G4InitXscPAI::DifPAIxSection

(

G4double

omega

)

Definition at line 414 of file G4InitXscPAI.cc.

References python.hepunit::electron_mass_c2, python.hepunit::fine_structure_const, python.hepunit::hbarc, ImPartDielectricConst(), IntegralTerm(), python.hepunit::pi, and RePartDielectricConst().

Referenced by DifPAIdEdx(), and IntegralPAIxSection().

{
  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 

G4int G4InitXscPAI::GetBinPAI ( ) const

inline

Definition at line 106 of file G4InitXscPAI.hh.

{ return fPAIbin ; }

G4PhysicsLogVector* G4InitXscPAI::GetChCosSqVector ( ) const

inline

Definition at line 117 of file G4InitXscPAI.hh.

{ return fChCosSqVector;}

G4PhysicsLogVector* G4InitXscPAI::GetChWidthVector ( ) const

inline

Definition at line 118 of file G4InitXscPAI.hh.

{ return fChWidthVector;}

G4int G4InitXscPAI::GetIntervalNumber ( ) const

inline

Definition at line 105 of file G4InitXscPAI.hh.

{ return fIntervalNumber ; }

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

inline

Definition at line 110 of file G4InitXscPAI.hh.

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

G4double G4InitXscPAI::GetNormalizationCof ( ) const

inline

Definition at line 108 of file G4InitXscPAI.hh.

{ return fNormalizationCof ; }

G4PhysicsLogVector* G4InitXscPAI::GetPAIdEdxVector ( ) const

inline

Definition at line 114 of file G4InitXscPAI.hh.

{ return fPAIdEdxVector;}

G4PhysicsLogVector* G4InitXscPAI::GetPAIelectronVector ( ) const

inline

Definition at line 116 of file G4InitXscPAI.hh.

{ return fPAIelectronVector;}

G4PhysicsLogVector* G4InitXscPAI::GetPAIphotonVector ( ) const

inline

Definition at line 115 of file G4InitXscPAI.hh.

{ return fPAIphotonVector;}

G4PhysicsLogVector* G4InitXscPAI::GetPAIxscVector ( ) const

inline

Definition at line 113 of file G4InitXscPAI.hh.

{ return fPAIxscVector;}

G4double G4InitXscPAI::GetPhotonLambda

(

G4double

omega

)

Definition at line 933 of file G4InitXscPAI.cc.

References G4cout, G4endl, and G4InuclParticleNames::lambda.

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

G4double G4InitXscPAI::GetStepCerenkovLoss

(

G4double

step

)

Definition at line 982 of file G4InitXscPAI.cc.

{  
  G4double loss = 0.0 ;
  loss *= step;

  return loss ;
}

G4double G4InitXscPAI::GetStepEnergyLoss

(

G4double

step

)

Definition at line 970 of file G4InitXscPAI.cc.

{  
  G4double loss = 0.0 ;
  loss *= step;

  return loss ;
}

G4double G4InitXscPAI::GetStepPlasmonLoss

(

G4double

step

)

Definition at line 994 of file G4InitXscPAI.cc.

{  


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

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

Definition at line 294 of file G4InitXscPAI.cc.

References python.hepunit::hbarc.

Referenced by DifPAIxSection(), IntegralCherenkov(), ModuleSqDielectricConst(), PAIdNdxCherenkov(), and PAIdNdxPlasmon().

{
   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 

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

Definition at line 758 of file G4InitXscPAI.cc.

References G4PhysicsVector::GetLowEdgeEnergy(), ImPartDielectricConst(), G4Integrator< T, F >::Legendre10(), ModuleSqDielectricConst(), PAIdNdxCherenkov(), G4PhysicsVector::PutValue(), RePartDielectricConst(), and width.

{
  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 

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

Definition at line 678 of file G4InitXscPAI.cc.

References DifPAIdEdx(), G4PhysicsVector::GetLowEdgeEnergy(), G4Integrator< T, F >::Legendre10(), and G4PhysicsVector::PutValue().

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

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

Definition at line 597 of file G4InitXscPAI.cc.

References DifPAIxSection(), G4PhysicsVector::GetLowEdgeEnergy(), G4Integrator< T, F >::Legendre10(), and G4PhysicsVector::PutValue().

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

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

Definition at line 854 of file G4InitXscPAI.cc.

References G4PhysicsVector::GetLowEdgeEnergy(), G4Integrator< T, F >::Legendre10(), PAIdNdxPlasmon(), and G4PhysicsVector::PutValue().

{
  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

G4double G4InitXscPAI::IntegralTerm

(

G4double

omega

)

Definition at line 256 of file G4InitXscPAI.cc.

References RutherfordIntegral().

Referenced by DifPAIxSection(), and PAIdNdxPlasmon().

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

void G4InitXscPAI::KillCloseIntervals

(

)

Definition at line 131 of file G4InitXscPAI.cc.

Referenced by G4InitXscPAI().

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

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

Definition at line 320 of file G4InitXscPAI.cc.

References ImPartDielectricConst(), and RePartDielectricConst().

Referenced by IntegralCherenkov().

{
   G4double eIm2, eRe2, result;

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

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

   result = eIm2 + eRe2;
  
   return result ;
}  

void G4InitXscPAI::Normalisation

(

)

Definition at line 162 of file G4InitXscPAI.cc.

References python.hepunit::electron_mass_c2, python.hepunit::fine_structure_const, python.hepunit::hbarc, python.hepunit::pi, and RutherfordIntegral().

Referenced by G4InitXscPAI().

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

G4double G4InitXscPAI::PAIdNdxCherenkov

(

G4double

omega

)

Definition at line 487 of file G4InitXscPAI.cc.

References python.hepunit::fine_structure_const, python.hepunit::hbarc, ImPartDielectricConst(), python.hepunit::pi, and RePartDielectricConst().

Referenced by IntegralCherenkov().

{        
  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 

G4double G4InitXscPAI::PAIdNdxPlasmon

(

G4double

omega

)

Definition at line 551 of file G4InitXscPAI.cc.

References python.hepunit::electron_mass_c2, python.hepunit::fine_structure_const, python.hepunit::hbarc, ImPartDielectricConst(), IntegralTerm(), python.hepunit::pi, and RePartDielectricConst().

Referenced by IntegralPlasmon().

{        
  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 

G4double G4InitXscPAI::RePartDielectricConst

(

G4double

energy

)

Definition at line 343 of file G4InitXscPAI.cc.

References plottest35::c1, python.hepunit::hbarc, and python.hepunit::pi.

Referenced by DifPAIxSection(), IntegralCherenkov(), ModuleSqDielectricConst(), PAIdNdxCherenkov(), and PAIdNdxPlasmon().

{
  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 

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

Definition at line 232 of file G4InitXscPAI.cc.

References plottest35::c1.

Referenced by IntegralTerm(), and Normalisation().

{
   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 

---

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

• G4InitXscPAI.hh
• G4InitXscPAI.cc