#include <G4GoudsmitSaundersonTable.hh>
Public Member Functions | |
G4GoudsmitSaundersonTable () | |
~G4GoudsmitSaundersonTable () | |
G4double | SampleTheta (G4double, G4double, G4double) |
Definition at line 48 of file G4GoudsmitSaundersonTable.hh.
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 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 }