Geant4: G4InitXscPAI Class Reference

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().

00071   : fPAIxscVector(NULL),
00072     fPAIdEdxVector(NULL),
00073     fPAIphotonVector(NULL),
00074     fPAIelectronVector(NULL),
00075     fChCosSqVector(NULL),
00076     fChWidthVector(NULL)
00077 {
00078   G4int i, j, matIndex;
00079  
00080   fDensity         = matCC->GetMaterial()->GetDensity();
00081   fElectronDensity = matCC->GetMaterial()->GetElectronDensity();
00082   matIndex         = matCC->GetMaterial()->GetIndex();
00083 
00084   fSandia          = new G4SandiaTable(matIndex);
00085   fIntervalNumber  = fSandia->GetMaxInterval()-1;
00086 
00087   fMatSandiaMatrix = new G4OrderedTable();
00088  
00089   for (i = 0; i < fIntervalNumber; i++)
00090   {
00091     fMatSandiaMatrix->push_back(new G4DataVector(5,0.));
00092   }                     
00093   for (i = 0; i < fIntervalNumber; i++)
00094   {
00095     (*(*fMatSandiaMatrix)[i])[0] = fSandia->GetSandiaMatTable(i,0);
00096 
00097     for(j = 1; j < 5 ; j++)
00098     {
00099       (*(*fMatSandiaMatrix)[i])[j] = fSandia->GetSandiaMatTable(i,j)*fDensity;
00100     }     
00101   }
00102   KillCloseIntervals();
00103   Normalisation();
00104   fBetaGammaSq = fTmax = 0.0;
00105   fIntervalTmax = fCurrentInterval = 0;
00106 }

G4InitXscPAI::~G4InitXscPAI ( ) [virtual]

Definition at line 115 of file G4InitXscPAI.cc.

00116 {
00117   if(fPAIxscVector)      delete fPAIxscVector;  
00118   if(fPAIdEdxVector)     delete fPAIdEdxVector;  
00119   if(fPAIphotonVector)   delete fPAIphotonVector;  
00120   if(fPAIelectronVector) delete fPAIelectronVector;  
00121   if(fChCosSqVector)     delete fChCosSqVector;  
00122   if(fChWidthVector)     delete fChWidthVector;  
00123   delete fSandia;
00124   delete fMatSandiaMatrix;
00125 }

G4double G4InitXscPAI::DifPAIdEdx ( G4double omega )

Definition at line 477 of file G4InitXscPAI.cc.

References DifPAIxSection().

Referenced by IntegralPAIdEdx().

00478 {
00479   G4double dEdx = omega*DifPAIxSection(omega);
00480   return dEdx;
00481 }

G4double G4InitXscPAI::DifPAIxSection ( G4double omega )

Definition at line 414 of file G4InitXscPAI.cc.

References ImPartDielectricConst(), IntegralTerm(), G4INCL::Math::pi, and RePartDielectricConst().

Referenced by DifPAIdEdx(), and IntegralPAIxSection().

00415 {
00416   G4int i = fCurrentInterval;
00417   G4double  betaGammaSq = fBetaGammaSq;       
00418   G4double integralTerm = IntegralTerm(omega);
00419   G4double be2,cof,x1,x2,x3,x4,x5,x6,x7,x8,result ;
00420   G4double epsilonRe = RePartDielectricConst(omega);
00421   G4double epsilonIm = ImPartDielectricConst(i,omega);
00422   G4double be4 ;
00423   G4double betaBohr2 = fine_structure_const*fine_structure_const ;
00424   G4double betaBohr4 = betaBohr2*betaBohr2*4.0 ;
00425   be2 = betaGammaSq/(1 + betaGammaSq) ;
00426   be4 = be2*be2 ;
00427  
00428    cof = 1 ;
00429    x1 = log(2*electron_mass_c2/omega) ;
00430 
00431    if( betaGammaSq < 0.01 ) x2 = log(be2) ;
00432    else
00433    {
00434      x2 = -log( (1/betaGammaSq - epsilonRe)*
00435                 (1/betaGammaSq - epsilonRe) + 
00436                 epsilonIm*epsilonIm )/2 ;
00437    }
00438    if( epsilonIm == 0.0 || betaGammaSq < 0.01 )
00439    {
00440      x6=0 ;
00441    }
00442    else
00443    {
00444      x3 = -epsilonRe + 1/betaGammaSq ;
00445      x5 = -1 - epsilonRe + be2*((1 +epsilonRe)*(1 + epsilonRe) +
00446            epsilonIm*epsilonIm) ;
00447 
00448      x7 = atan2(epsilonIm,x3) ;
00449      x6 = x5 * x7 ;
00450    }
00451     // if(fImPartDielectricConst[i] == 0) x6 = 0 ;
00452    
00453    x4 = ((x1 + x2)*epsilonIm + x6)/hbarc ;
00454    //   if( x4 < 0.0 ) x4 = 0.0 ;
00455    x8 = (1 + epsilonRe)*(1 + epsilonRe) + 
00456         epsilonIm*epsilonIm;
00457 
00458    result = (x4 + cof*integralTerm/omega/omega) ;
00459    if(result < 1.0e-8) result = 1.0e-8 ;
00460    result *= fine_structure_const/be2/pi ;
00461    //   result *= (1-exp(-beta/betaBohr))*(1-exp(-beta/betaBohr)) ;
00462    //  result *= (1-exp(-be2/betaBohr2)) ;
00463    result *= (1-exp(-be4/betaBohr4)) ;
00464    if(fDensity >= fSolidDensity)
00465    { 
00466       result /= x8 ;
00467    }
00468    return result ;
00469 
00470 } // end of DifPAIxSection

G4int G4InitXscPAI::GetBinPAI ( ) const [inline]

Definition at line 106 of file G4InitXscPAI.hh.

00106 { return fPAIbin ; }

G4PhysicsLogVector* G4InitXscPAI::GetChCosSqVector ( ) const [inline]

Definition at line 117 of file G4InitXscPAI.hh.

00117 { return fChCosSqVector;}

G4PhysicsLogVector* G4InitXscPAI::GetChWidthVector ( ) const [inline]

Definition at line 118 of file G4InitXscPAI.hh.

00118 { return fChWidthVector;}

G4int G4InitXscPAI::GetIntervalNumber ( ) const [inline]

Definition at line 105 of file G4InitXscPAI.hh.

00105 { return fIntervalNumber ; }

G4double G4InitXscPAI::GetMatSandiaMatrix	(	G4int	i,
		G4int	j
	)			const `[inline]`

Definition at line 110 of file G4InitXscPAI.hh.

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

G4double G4InitXscPAI::GetNormalizationCof ( ) const [inline]

Definition at line 108 of file G4InitXscPAI.hh.

00108 { return fNormalizationCof ; }

G4PhysicsLogVector* G4InitXscPAI::GetPAIdEdxVector ( ) const [inline]

Definition at line 114 of file G4InitXscPAI.hh.

00114 { return fPAIdEdxVector;}

G4PhysicsLogVector* G4InitXscPAI::GetPAIelectronVector ( ) const [inline]

Definition at line 116 of file G4InitXscPAI.hh.

00116 { return fPAIelectronVector;}

G4PhysicsLogVector* G4InitXscPAI::GetPAIphotonVector ( ) const [inline]

Definition at line 115 of file G4InitXscPAI.hh.

00115 { return fPAIphotonVector;}

G4PhysicsLogVector* G4InitXscPAI::GetPAIxscVector ( ) const [inline]

Definition at line 113 of file G4InitXscPAI.hh.

00113 { return fPAIxscVector;}

G4double G4InitXscPAI::GetPhotonLambda ( G4double omega )

Definition at line 933 of file G4InitXscPAI.cc.

References G4cout, G4endl, and G4InuclParticleNames::lambda.

00934 {  
00935   G4int i ;
00936   G4double omega2, omega3, omega4, a1, a2, a3, a4, lambda ;
00937 
00938   omega2 = omega*omega ;
00939   omega3 = omega2*omega ;
00940   omega4 = omega2*omega2 ;
00941 
00942   for(i = 0; i < fIntervalNumber;i++)
00943   {
00944     if( omega < (*(*fMatSandiaMatrix)[i])[0] ) break ;
00945   }
00946   if( i == 0 )
00947   {
00948     G4cout<<"Warning: energy in G4InitXscPAI::GetPhotonLambda < I1"<<G4endl;
00949   }
00950   else i-- ;
00951 
00952   a1 = (*(*fMatSandiaMatrix)[i])[1]; 
00953   a2 = (*(*fMatSandiaMatrix)[i])[2]; 
00954   a3 = (*(*fMatSandiaMatrix)[i])[3]; 
00955   a4 = (*(*fMatSandiaMatrix)[i])[4]; 
00956 
00957   lambda = 1./(a1/omega + a2/omega2 + a3/omega3 + a4/omega4);
00958 
00959   return lambda ;
00960 }

G4double G4InitXscPAI::GetStepCerenkovLoss ( G4double step )

Definition at line 982 of file G4InitXscPAI.cc.

00983 {  
00984   G4double loss = 0.0 ;
00985   loss *= step;
00986 
00987   return loss ;
00988 }

G4double G4InitXscPAI::GetStepEnergyLoss ( G4double step )

Definition at line 970 of file G4InitXscPAI.cc.

00971 {  
00972   G4double loss = 0.0 ;
00973   loss *= step;
00974 
00975   return loss ;
00976 }

G4double G4InitXscPAI::GetStepPlasmonLoss ( G4double step )

Definition at line 994 of file G4InitXscPAI.cc.

00995 {  
00996 
00997 
00998   G4double loss = 0.0 ;
00999   loss *= step;
01000   return loss ;
01001 }

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

Definition at line 294 of file G4InitXscPAI.cc.

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

00296 {
00297    G4double energy2,energy3,energy4,a1,a2,a3,a4,result;
00298 
00299    a1 = (*(*fMatSandiaMatrix)[k])[1]; 
00300    a2 = (*(*fMatSandiaMatrix)[k])[2]; 
00301    a3 = (*(*fMatSandiaMatrix)[k])[3]; 
00302    a4 = (*(*fMatSandiaMatrix)[k])[4]; 
00303 
00304    energy2 = energy1*energy1;
00305    energy3 = energy2*energy1;
00306    energy4 = energy3*energy1;
00307    
00308    result  = a1/energy1+a2/energy2+a3/energy3+a4/energy4 ;  
00309    result *= hbarc/energy1 ;
00310    
00311    return result ;
00312 
00313 }  // 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(), and RePartDielectricConst().

00759 {
00760   G4int i,k,i1,i2;
00761   G4double energy1, energy2, beta2, module2, cos2, width, result = 0.;
00762  
00763   fBetaGammaSq = bg2;
00764   fTmax        = Tmax;
00765   beta2        = bg2/(1+bg2);
00766 
00767   if(fPAIphotonVector) delete fPAIphotonVector;  
00768   if(fChCosSqVector)   delete fChCosSqVector;  
00769   if(fChWidthVector)   delete fChWidthVector;  
00770   
00771   fPAIphotonVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
00772   fChCosSqVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
00773   fChWidthVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
00774 
00775   fPAIphotonVector->PutValue(fPAIbin-1,result);
00776   fChCosSqVector->PutValue(fPAIbin-1,1.);
00777   fChWidthVector->PutValue(fPAIbin-1,1e-7);
00778                         
00779   for( i = fIntervalNumber - 1; i >= 0; i-- )
00780   {
00781     if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
00782   }
00783   if (i < 0) i = 0; // Tmax should be more than 
00784                     // first ionisation potential
00785   fIntervalTmax = i;
00786 
00787   G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;
00788 
00789   for( k = fPAIbin - 2; k >= 0; k-- )
00790   {
00791     energy1 = fPAIphotonVector->GetLowEdgeEnergy(k);
00792     energy2 = fPAIphotonVector->GetLowEdgeEnergy(k+1);
00793 
00794     for( i = fIntervalTmax; i >= 0; i-- ) 
00795     {
00796       if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
00797     }
00798     if(i < 0) i = 0;
00799     i2 = i;
00800 
00801     for( i = fIntervalTmax; i >= 0; i-- ) 
00802     {
00803       if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
00804     }
00805     if(i < 0) i = 0;
00806     i1 = i;
00807 
00808     module2 = ModuleSqDielectricConst(i1,energy1);
00809     cos2    = RePartDielectricConst(energy1)/module2/beta2;      
00810     width   = ImPartDielectricConst(i1,energy1)/module2/beta2;
00811       
00812     fChCosSqVector->PutValue(k,cos2);
00813     fChWidthVector->PutValue(k,width);
00814 
00815     if( i1 == i2 )
00816     {
00817       fCurrentInterval = i1;
00818       result += integral.Legendre10(this,&G4InitXscPAI::PAIdNdxCherenkov,
00819                                     energy1,energy2);
00820       fPAIphotonVector->PutValue(k,result);
00821 
00822     }
00823     else
00824     {
00825       for( i = i2; i >= i1; i-- ) 
00826       {
00827         fCurrentInterval = i;
00828 
00829         if( i==i2 )        result += integral.Legendre10(this,
00830                            &G4InitXscPAI::PAIdNdxCherenkov,
00831                            (*(*fMatSandiaMatrix)[i])[0] ,energy2);
00832 
00833         else if( i == i1 ) result += integral.Legendre10(this,
00834                            &G4InitXscPAI::PAIdNdxCherenkov,energy1,
00835                            (*(*fMatSandiaMatrix)[i+1])[0]);
00836 
00837         else               result += integral.Legendre10(this,
00838                            &G4InitXscPAI::PAIdNdxCherenkov,
00839                        (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
00840       }
00841       fPAIphotonVector->PutValue(k,result);
00842     }
00843     // G4cout<<k<<"\t"<<result<<G4endl; 
00844   }
00845   return;
00846 }   // 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().

00679 {
00680   G4int i,k,i1,i2;
00681   G4double energy1, energy2, result = 0.;
00682  
00683   fBetaGammaSq = bg2;
00684   fTmax        = Tmax;
00685 
00686   if(fPAIdEdxVector) delete fPAIdEdxVector;  
00687   
00688   fPAIdEdxVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
00689   fPAIdEdxVector->PutValue(fPAIbin-1,result);
00690                         
00691   for( i = fIntervalNumber - 1; i >= 0; i-- )
00692   {
00693     if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
00694   }
00695   if (i < 0) i = 0; // Tmax should be more than 
00696                     // first ionisation potential
00697   fIntervalTmax = i;
00698 
00699   G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;
00700 
00701   for( k = fPAIbin - 2; k >= 0; k-- )
00702   {
00703     energy1 = fPAIdEdxVector->GetLowEdgeEnergy(k);
00704     energy2 = fPAIdEdxVector->GetLowEdgeEnergy(k+1);
00705 
00706     for( i = fIntervalTmax; i >= 0; i-- ) 
00707     {
00708       if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
00709     }
00710     if(i < 0) i = 0;
00711     i2 = i;
00712 
00713     for( i = fIntervalTmax; i >= 0; i-- ) 
00714     {
00715       if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
00716     }
00717     if(i < 0) i = 0;
00718     i1 = i;
00719 
00720     if( i1 == i2 )
00721     {
00722       fCurrentInterval = i1;
00723       result += integral.Legendre10(this,&G4InitXscPAI::DifPAIdEdx,
00724                                     energy1,energy2);
00725       fPAIdEdxVector->PutValue(k,result);
00726     }
00727     else
00728     {
00729       for( i = i2; i >= i1; i-- ) 
00730       {
00731         fCurrentInterval = i;
00732 
00733         if( i==i2 )        result += integral.Legendre10(this,
00734                            &G4InitXscPAI::DifPAIdEdx,
00735                            (*(*fMatSandiaMatrix)[i])[0] ,energy2);
00736 
00737         else if( i == i1 ) result += integral.Legendre10(this,
00738                            &G4InitXscPAI::DifPAIdEdx,energy1,
00739                            (*(*fMatSandiaMatrix)[i+1])[0]);
00740 
00741         else               result += integral.Legendre10(this,
00742                            &G4InitXscPAI::DifPAIdEdx,
00743                        (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
00744       }
00745       fPAIdEdxVector->PutValue(k,result);
00746     }
00747     // G4cout<<k<<"\t"<<result<<G4endl; 
00748   }
00749   return ;
00750 }

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().

00598 {
00599   G4int i,k,i1,i2;
00600   G4double energy1, energy2, result = 0.;
00601  
00602   fBetaGammaSq = bg2;
00603   fTmax        = Tmax;
00604 
00605   if(fPAIxscVector) delete fPAIxscVector;  
00606   
00607   fPAIxscVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
00608   fPAIxscVector->PutValue(fPAIbin-1,result);
00609                         
00610   for( i = fIntervalNumber - 1; i >= 0; i-- )
00611   {
00612     if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
00613   }
00614   if (i < 0) i = 0; // Tmax should be more than 
00615                     // first ionisation potential
00616   fIntervalTmax = i;
00617 
00618   G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;
00619 
00620   for( k = fPAIbin - 2; k >= 0; k-- )
00621   {
00622     energy1 = fPAIxscVector->GetLowEdgeEnergy(k);
00623     energy2 = fPAIxscVector->GetLowEdgeEnergy(k+1);
00624 
00625     for( i = fIntervalTmax; i >= 0; i-- ) 
00626     {
00627       if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
00628     }
00629     if(i < 0) i = 0;
00630     i2 = i;
00631 
00632     for( i = fIntervalTmax; i >= 0; i-- ) 
00633     {
00634       if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
00635     }
00636     if(i < 0) i = 0;
00637     i1 = i;
00638 
00639     if( i1 == i2 )
00640     {
00641       fCurrentInterval = i1;
00642       result += integral.Legendre10(this,&G4InitXscPAI::DifPAIxSection,
00643                                     energy1,energy2);
00644       fPAIxscVector->PutValue(k,result);
00645     }
00646     else
00647     {
00648       for( i = i2; i >= i1; i-- ) 
00649       {
00650         fCurrentInterval = i;
00651 
00652         if( i==i2 )        result += integral.Legendre10(this,
00653                            &G4InitXscPAI::DifPAIxSection,
00654                            (*(*fMatSandiaMatrix)[i])[0] ,energy2);
00655 
00656         else if( i == i1 ) result += integral.Legendre10(this,
00657                            &G4InitXscPAI::DifPAIxSection,energy1,
00658                            (*(*fMatSandiaMatrix)[i+1])[0]);
00659 
00660         else               result += integral.Legendre10(this,
00661                            &G4InitXscPAI::DifPAIxSection,
00662                        (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
00663       }
00664       fPAIxscVector->PutValue(k,result);
00665     }
00666     // G4cout<<k<<"\t"<<result<<G4endl; 
00667   }
00668   return ;
00669 }

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().

00855 {
00856   G4int i,k,i1,i2;
00857   G4double energy1, energy2, result = 0.;
00858  
00859   fBetaGammaSq = bg2;
00860   fTmax        = Tmax;
00861 
00862   if(fPAIelectronVector) delete fPAIelectronVector;  
00863   
00864   fPAIelectronVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
00865   fPAIelectronVector->PutValue(fPAIbin-1,result);
00866                         
00867   for( i = fIntervalNumber - 1; i >= 0; i-- )
00868   {
00869     if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
00870   }
00871   if (i < 0) i = 0; // Tmax should be more than 
00872                     // first ionisation potential
00873   fIntervalTmax = i;
00874 
00875   G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;
00876 
00877   for( k = fPAIbin - 2; k >= 0; k-- )
00878   {
00879     energy1 = fPAIelectronVector->GetLowEdgeEnergy(k);
00880     energy2 = fPAIelectronVector->GetLowEdgeEnergy(k+1);
00881 
00882     for( i = fIntervalTmax; i >= 0; i-- ) 
00883     {
00884       if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
00885     }
00886     if(i < 0) i = 0;
00887     i2 = i;
00888 
00889     for( i = fIntervalTmax; i >= 0; i-- ) 
00890     {
00891       if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
00892     }
00893     if(i < 0) i = 0;
00894     i1 = i;
00895 
00896     if( i1 == i2 )
00897     {
00898       fCurrentInterval = i1;
00899       result += integral.Legendre10(this,&G4InitXscPAI::PAIdNdxPlasmon,
00900                                     energy1,energy2);
00901       fPAIelectronVector->PutValue(k,result);
00902     }
00903     else
00904     {
00905       for( i = i2; i >= i1; i-- ) 
00906       {
00907         fCurrentInterval = i;
00908 
00909         if( i==i2 )        result += integral.Legendre10(this,
00910                            &G4InitXscPAI::PAIdNdxPlasmon,
00911                            (*(*fMatSandiaMatrix)[i])[0] ,energy2);
00912 
00913         else if( i == i1 ) result += integral.Legendre10(this,
00914                            &G4InitXscPAI::PAIdNdxPlasmon,energy1,
00915                            (*(*fMatSandiaMatrix)[i+1])[0]);
00916 
00917         else               result += integral.Legendre10(this,
00918                            &G4InitXscPAI::PAIdNdxPlasmon,
00919                        (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
00920       }
00921       fPAIelectronVector->PutValue(k,result);
00922     }
00923     // G4cout<<k<<"\t"<<result<<G4endl; 
00924   }
00925   return;
00926 }   // end of IntegralPlasmon

G4double G4InitXscPAI::IntegralTerm ( G4double omega )

Definition at line 256 of file G4InitXscPAI.cc.

References RutherfordIntegral().

Referenced by DifPAIxSection(), and PAIdNdxPlasmon().

00257 {
00258   G4int i;
00259   G4double energy1, energy2, result = 0.; 
00260                         
00261   for( i = 0; i <= fIntervalTmax; i++ )
00262   {
00263     if(i == fIntervalTmax) 
00264     {
00265       energy1 = (*(*fMatSandiaMatrix)[i])[0];
00266       result += RutherfordIntegral(i,energy1,omega);
00267     }
00268     else 
00269     {
00270       if( omega <= (*(*fMatSandiaMatrix)[i+1])[0])
00271       {
00272         energy1 = (*(*fMatSandiaMatrix)[i])[0];
00273         result += RutherfordIntegral(i,energy1,omega);
00274         break;
00275       }
00276       else
00277       {
00278         energy1 = (*(*fMatSandiaMatrix)[i])[0];
00279         energy2 = (*(*fMatSandiaMatrix)[i+1])[0];
00280         result += RutherfordIntegral(i,energy1,energy2);
00281       }
00282     }
00283     // G4cout<<"IntegralTerm<<"("<<omega<<")"<<" = "<<result<<G4endl;
00284   }
00285   return result;
00286 }

void G4InitXscPAI::KillCloseIntervals ( )

Definition at line 131 of file G4InitXscPAI.cc.

Referenced by G4InitXscPAI().

00132 {
00133   G4int i, j, k;
00134   G4double energy1, energy2; 
00135                         
00136   for( i = 0 ; i < fIntervalNumber - 1 ; i++ )
00137   {
00138     energy1 = (*(*fMatSandiaMatrix)[i])[0];
00139     energy2 = (*(*fMatSandiaMatrix)[i+1])[0];
00140 
00141     if( energy2 - energy1 > 1.5*fDelta*(energy1 + energy2) )  continue ;
00142     else
00143     {
00144       for(j = i; j < fIntervalNumber-1; j++)
00145       {
00146         for( k = 0; k < 5; k++ )
00147         {
00148           (*(*fMatSandiaMatrix)[j])[k] = (*(*fMatSandiaMatrix)[j+1])[k];
00149         }
00150       }
00151       fIntervalNumber-- ;
00152       i-- ;
00153     }
00154   }
00155   
00156 }

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

Definition at line 320 of file G4InitXscPAI.cc.

References ImPartDielectricConst(), and RePartDielectricConst().

Referenced by IntegralCherenkov().

00322 {
00323    G4double eIm2, eRe2, result;
00324 
00325    result = ImPartDielectricConst(k,omega);
00326    eIm2   = result*result;
00327 
00328    result = RePartDielectricConst(omega);
00329    eRe2   = result*result;
00330 
00331    result = eIm2 + eRe2;
00332   
00333    return result ;
00334 }

void G4InitXscPAI::Normalisation ( )

Definition at line 162 of file G4InitXscPAI.cc.

References G4INCL::Math::pi, and RutherfordIntegral().

Referenced by G4InitXscPAI().

00163 {
00164   G4int i, j;
00165   G4double energy1, energy2, /*delta,*/ cof; // , shift;
00166 
00167   energy1 = (*(*fMatSandiaMatrix)[fIntervalNumber-1])[0];
00168   energy2 = 2.*(*(*fMatSandiaMatrix)[fIntervalNumber-1])[0];
00169 
00170  
00171   cof = RutherfordIntegral(fIntervalNumber-1,energy1,energy2);
00172                         
00173   for( i = fIntervalNumber-2; i >= 0; i-- )
00174   {
00175     energy1 = (*(*fMatSandiaMatrix)[i])[0];
00176     energy2 = (*(*fMatSandiaMatrix)[i+1])[0];
00177 
00178     cof += RutherfordIntegral(i,energy1,energy2);
00179     // G4cout<<"norm. cof = "<<cof<<G4endl;
00180   }
00181   fNormalizationCof  = 2*pi*pi*hbarc*hbarc*fine_structure_const/electron_mass_c2 ;
00182   fNormalizationCof *= fElectronDensity;
00183   //delta = fNormalizationCof - cof;
00184   fNormalizationCof /= cof;
00185   //  G4cout<<"G4InitXscPAI::fNormalizationCof/cof = "<<fNormalizationCof
00186   //    <<";  at delta ="<<delta<<G4endl ;
00187 
00188   for (i = 0; i < fIntervalNumber; i++) // renormalisation on QM sum rule
00189   {
00190     for(j = 1; j < 5 ; j++)
00191     {
00192       (*(*fMatSandiaMatrix)[i])[j] *= fNormalizationCof;
00193     }     
00194   }
00195   /* 
00196   if(delta > 0) // shift the first energy interval
00197   {
00198     for(i=1;i<100;i++)
00199     {
00200       energy1 = (1.-i/100.)*(*(*fMatSandiaMatrix)[0])[0];
00201       energy2 = (*(*fMatSandiaMatrix)[0])[0];
00202       shift   = RutherfordIntegral(0,energy1,energy2);
00203       G4cout<<shift<<"\t";
00204       if(shift >= delta) break;
00205     }
00206     (*(*fMatSandiaMatrix)[0])[0] = energy1;
00207     cof += shift;
00208   }
00209   else if(delta < 0)
00210   {
00211     for(i=1;i<100;i++)
00212     {
00213       energy1 = (*(*fMatSandiaMatrix)[0])[0];
00214       energy2 = (*(*fMatSandiaMatrix)[0])[0] + 
00215               ( (*(*fMatSandiaMatrix)[0])[0] - (*(*fMatSandiaMatrix)[0])[0] )*i/100.;
00216       shift   = RutherfordIntegral(0,energy1,energy2);
00217       if( shift >= std::abs(delta) ) break;
00218     }
00219     (*(*fMatSandiaMatrix)[0])[0] = energy2;
00220     cof -= shift;
00221   }
00222   G4cout<<G4cout<<"G4InitXscPAI::fNormalizationCof/cof = "<<fNormalizationCof/cof
00223         <<";  at delta ="<<delta<<"  and i = "<<i<<G4endl ;
00224  */ 
00225 }

G4double G4InitXscPAI::PAIdNdxCherenkov ( G4double omega )

Definition at line 487 of file G4InitXscPAI.cc.

References ImPartDielectricConst(), G4INCL::Math::pi, and RePartDielectricConst().

Referenced by IntegralCherenkov().

00488 {        
00489   G4int i = fCurrentInterval;
00490   G4double  betaGammaSq = fBetaGammaSq;       
00491   G4double epsilonRe = RePartDielectricConst(omega);
00492   G4double epsilonIm = ImPartDielectricConst(i,omega);
00493 
00494   G4double /*cof,*/ logarithm, x3, x5, argument, modul2, dNdxC ; 
00495   G4double be2, be4, betaBohr2,betaBohr4,cofBetaBohr ;
00496 
00497   //cof         = 1.0 ;
00498    cofBetaBohr = 4.0 ;
00499    betaBohr2   = fine_structure_const*fine_structure_const ;
00500    betaBohr4   = betaBohr2*betaBohr2*cofBetaBohr ;
00501 
00502    be2 = betaGammaSq/(1 + betaGammaSq) ;
00503    be4 = be2*be2 ;
00504 
00505    if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq) ; // 0.0 ;
00506    else
00507    {
00508      logarithm  = -log( (1/betaGammaSq - epsilonRe)*
00509                         (1/betaGammaSq - epsilonRe) + 
00510                         epsilonIm*epsilonIm )*0.5 ;
00511      logarithm += log(1+1.0/betaGammaSq) ;
00512    }
00513 
00514    if( epsilonIm == 0.0 || betaGammaSq < 0.01 )
00515    {
00516      argument = 0.0 ;
00517    }
00518    else
00519    {
00520      x3 = -epsilonRe + 1.0/betaGammaSq ;
00521      x5 = -1.0 - epsilonRe +
00522           be2*((1.0 +epsilonRe)*(1.0 + epsilonRe) +
00523           epsilonIm*epsilonIm) ;
00524      if( x3 == 0.0 ) argument = 0.5*pi;
00525      else            argument = atan2(epsilonIm,x3) ;
00526      argument *= x5  ;
00527    }   
00528    dNdxC = ( logarithm*epsilonIm + argument )/hbarc ;
00529   
00530    if(dNdxC < 1.0e-8) dNdxC = 1.0e-8 ;
00531 
00532    dNdxC *= fine_structure_const/be2/pi ;
00533 
00534    dNdxC *= (1-exp(-be4/betaBohr4)) ;
00535 
00536    if(fDensity >= fSolidDensity)
00537    { 
00538       modul2 = (1.0 + epsilonRe)*(1.0 + epsilonRe) + 
00539                     epsilonIm*epsilonIm;
00540       dNdxC /= modul2 ;
00541    }
00542    return dNdxC ;
00543 
00544 } // end of PAIdNdxCerenkov

G4double G4InitXscPAI::PAIdNdxPlasmon ( G4double omega )

Definition at line 551 of file G4InitXscPAI.cc.

References ImPartDielectricConst(), IntegralTerm(), G4INCL::Math::pi, and RePartDielectricConst().

Referenced by IntegralPlasmon().

00552 {        
00553   G4int i = fCurrentInterval;
00554   G4double  betaGammaSq = fBetaGammaSq;       
00555   G4double integralTerm = IntegralTerm(omega);
00556   G4double epsilonRe = RePartDielectricConst(omega);
00557   G4double epsilonIm = ImPartDielectricConst(i,omega);
00558 
00559    G4double cof, resonance, modul2, dNdxP ;
00560    G4double be2, be4, betaBohr2, betaBohr4, cofBetaBohr ;
00561 
00562    cof = 1 ;
00563    cofBetaBohr = 4.0 ;
00564    betaBohr2   = fine_structure_const*fine_structure_const ;
00565    betaBohr4   = betaBohr2*betaBohr2*cofBetaBohr ;
00566 
00567    be2 = betaGammaSq/(1 + betaGammaSq) ;
00568    be4 = be2*be2 ;
00569  
00570    resonance  = log(2*electron_mass_c2*be2/omega) ;  
00571    resonance *= epsilonIm/hbarc ;
00572 
00573 
00574    dNdxP = ( resonance + cof*integralTerm/omega/omega ) ;
00575 
00576    if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8 ;
00577 
00578    dNdxP *= fine_structure_const/be2/pi ;
00579    dNdxP *= (1-exp(-be4/betaBohr4)) ;
00580 
00581    if( fDensity >= fSolidDensity )
00582    { 
00583      modul2 = (1 + epsilonRe)*(1 + epsilonRe) + 
00584                epsilonIm*epsilonIm;
00585      dNdxP /= modul2 ;
00586    }
00587    return dNdxP ;
00588 
00589 } // end of PAIdNdxPlasmon

G4double G4InitXscPAI::RePartDielectricConst ( G4double energy )

Definition at line 343 of file G4InitXscPAI.cc.

References G4INCL::Math::pi.

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

00344 {
00345   G4int i;       
00346    G4double x0, x02, x03, x04, x05, x1, x2, a1,a2,a3,a4,xx1 ,xx2 , xx12,
00347             c1, c2, c3, cof1, cof2, xln1, xln2, xln3, result ;
00348 
00349    x0 = enb ;
00350    result = 0 ;
00351    
00352    for( i = 0; i < fIntervalNumber-1; i++)
00353    {
00354       x1 = (*(*fMatSandiaMatrix)[i])[0];
00355       x2 = (*(*fMatSandiaMatrix)[i+1])[0] ;
00356 
00357       a1 = (*(*fMatSandiaMatrix)[i])[1]; 
00358       a2 = (*(*fMatSandiaMatrix)[i])[2]; 
00359       a3 = (*(*fMatSandiaMatrix)[i])[3]; 
00360       a4 = (*(*fMatSandiaMatrix)[i])[4];
00361  
00362       if( std::abs(x0-x1) < 0.5*(x0+x1)*fDelta ) 
00363       {
00364         if(x0 >= x1) x0 = x1*(1+fDelta);
00365         else         x0 = x1*(1-fDelta);
00366       } 
00367       if( std::abs(x0-x2) < 0.5*(x0+x2)*fDelta ) 
00368       {
00369         if(x0 >= x2) x0 = x2*(1+fDelta);
00370         else         x0 = x2*(1-fDelta);
00371       }
00372       xx1 = x1 - x0 ;
00373       xx2 = x2 - x0 ;
00374       xx12 = xx2/xx1 ;
00375       
00376       if( xx12 < 0 ) xx12 = -xx12;
00377       
00378       xln1 = log(x2/x1) ;
00379       xln2 = log(xx12) ;
00380       xln3 = log((x2 + x0)/(x1 + x0)) ;
00381 
00382       x02 = x0*x0 ;
00383       x03 = x02*x0 ;
00384       x04 = x03*x0 ;
00385       x05 = x04*x0;
00386 
00387       c1  = (x2 - x1)/x1/x2 ;
00388       c2  = (x2 - x1)*(x2 +x1)/x1/x1/x2/x2 ;
00389       c3  = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2 ;
00390 
00391       result -= (a1/x02 + a3/x04)*xln1 ;
00392       result -= (a2/x02 + a4/x04)*c1 ;
00393       result -= a3*c2/2/x02 ;
00394       result -= a4*c3/3/x02 ;
00395 
00396       cof1 = a1/x02 + a3/x04 ;
00397       cof2 = a2/x03 + a4/x05 ;
00398 
00399       result += 0.5*(cof1 +cof2)*xln2 ;
00400       result += 0.5*(cof1 - cof2)*xln3 ;
00401    } 
00402    result *= 2*hbarc/pi ;
00403    
00404    return result ;
00405 
00406 }   // end of RePartDielectricConst

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

Definition at line 232 of file G4InitXscPAI.cc.

Referenced by IntegralTerm(), and Normalisation().

00235 {
00236    G4double  c1, c2, c3, a1, a2, a3, a4 ;
00237 
00238    a1 = (*(*fMatSandiaMatrix)[k])[1]; 
00239    a2 = (*(*fMatSandiaMatrix)[k])[2]; 
00240    a3 = (*(*fMatSandiaMatrix)[k])[3]; 
00241    a4 = (*(*fMatSandiaMatrix)[k])[4]; 
00242    // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;   
00243    c1 = (x2 - x1)/x1/x2 ;
00244    c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2 ;
00245    c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2 ;
00246    // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = "<<c2<<"; "<<"c3 = "<<c3<<G4endl;   
00247    
00248    return  a1*log(x2/x1) + a2*c1 + a3*c2/2 + a4*c3/3 ;
00249 
00250 }   // end of RutherfordIntegral


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

G4InitXscPAI Class Reference

Public Member Functions

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation