G4GoudsmitSaundersonTable Class Reference

#include <G4GoudsmitSaundersonTable.hh>


Public Member Functions

 G4GoudsmitSaundersonTable ()
 ~G4GoudsmitSaundersonTable ()
G4double SampleTheta (G4double, G4double, G4double)


Detailed Description

Definition at line 48 of file G4GoudsmitSaundersonTable.hh.


Constructor & Destructor Documentation

G4GoudsmitSaundersonTable::G4GoudsmitSaundersonTable (  ) 

Definition at line 180 of file G4GoudsmitSaundersonTable.cc.

References G4cout, and G4endl.

00181 {
00182   if(PDF == 0){ 
00183     G4cout << "### G4GoudsmitSaundersonTable loading PDF data" << G4endl;
00184 
00185     PDF = new G4double [76*11*320];
00186     CPDF= new G4double [76*11*320];
00187 
00188     LoadPDFandCPDFdata();
00189   }
00190 }

G4GoudsmitSaundersonTable::~G4GoudsmitSaundersonTable (  ) 

Definition at line 193 of file G4GoudsmitSaundersonTable.cc.

00194 {
00195   if(PDF) {
00196     delete [] PDF;
00197     delete [] CPDF;
00198     PDF = 0;
00199   }
00200 }


Member Function Documentation

G4double G4GoudsmitSaundersonTable::SampleTheta ( G4double  ,
G4double  ,
G4double   
)

Definition at line 203 of file G4GoudsmitSaundersonTable.cc.

00204 { 
00205   //Benedito's procedure
00206   G4double A[11],ThisPDF[320],ThisCPDF[320];
00207   G4double coeff,Ckj,CkjPlus1,CkPlus1j,CkPlus1jPlus1,a,b,mmm,F;
00208   G4int Ind0,KIndex=0,JIndex=0,IIndex=0;
00209 
00210 
00212   // Find Lambda and Chia2 Index
00213   for(G4int k=0;k<75;k++){if((lambda>=LAMBDAN[k])&&(lambda<LAMBDAN[k+1])){KIndex=k;break;}}  
00214   for(G4int j=0;j<11;j++){A[j]=scrA[11*KIndex+j];}  
00215   for(G4int j=0;j<10;j++){if((Chia2>=A[j])&&(Chia2<A[j+1])){JIndex=j;break;}}
00216 
00218   // Calculate some necessary coefficients for PDF and CPDF interpolation 
00219   coeff=(std::log(LAMBDAN[KIndex+1]/LAMBDAN[KIndex]))*(std::log(A[JIndex+1]/A[JIndex]));
00220   Ckj=(std::log(LAMBDAN[KIndex+1]/lambda))*(std::log(A[JIndex+1]/Chia2))/coeff;
00221   CkjPlus1=(std::log(LAMBDAN[KIndex+1]/lambda))*(std::log(Chia2/A[JIndex]))/coeff;
00222   CkPlus1j=(std::log(lambda/LAMBDAN[KIndex]))*(std::log(A[JIndex+1]/Chia2))/coeff;
00223   CkPlus1jPlus1=(std::log(lambda/LAMBDAN[KIndex]))*(std::log(Chia2/A[JIndex]))/coeff;
00225   // Calculate Interpolated PDF and CPDF arrays
00226   Ind0=320*(11*KIndex+JIndex);
00227   for(G4int i=0 ; i<320 ;i++){
00228     ThisPDF[i]=Ckj*PDF[Ind0]+CkjPlus1*PDF[Ind0+320]+CkPlus1j*PDF[Ind0+3520]+CkPlus1jPlus1*PDF[Ind0+3840];
00229     ThisCPDF[i]=Ckj*CPDF[Ind0]+CkjPlus1*CPDF[Ind0+320]+CkPlus1j*CPDF[Ind0+3520]+CkPlus1jPlus1*CPDF[Ind0+3840];  
00230   // Find u Index using secant method
00231     if((i!=0)&&((rndm>=ThisCPDF[i-1])&&(rndm<ThisCPDF[i])))  {IIndex=i-1;break;}
00232     Ind0++;
00233   }
00234 
00236   //CPDF^-1(rndm)=x ==> CPDF(x)=rndm;
00237   a=uvalues[IIndex];
00238   b=uvalues[IIndex+1];
00239   
00240   do{
00241     mmm=0.5*(a+b);
00242     F=(ThisCPDF[IIndex]+(mmm-uvalues[IIndex])*ThisPDF[IIndex]
00243        +((mmm-uvalues[IIndex])*(mmm-uvalues[IIndex])
00244          *(ThisPDF[IIndex+1]-ThisPDF[IIndex]))
00245        /(2.*(uvalues[IIndex+1]-uvalues[IIndex])))-rndm;
00246     if(F>0.)b=mmm;
00247     else a=mmm;
00248   } while(std::fabs(b-a)>1.0e-9);
00249 
00250   return mmm;
00251 }


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:52:07 2013 for Geant4 by  doxygen 1.4.7