Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions
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 182 of file G4GoudsmitSaundersonTable.cc.

References G4cout, and G4endl.

183 {
184  if(PDF == 0){
185  G4cout << "### G4GoudsmitSaundersonTable loading PDF data" << G4endl;
186 
187  PDF = new G4double [76*11*320];
188  CPDF= new G4double [76*11*320];
189 
190  LoadPDFandCPDFdata();
191  }
192 }
G4GLOB_DLL std::ostream G4cout
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4GoudsmitSaundersonTable::~G4GoudsmitSaundersonTable ( )

Definition at line 195 of file G4GoudsmitSaundersonTable.cc.

196 {
197  if(PDF) {
198  delete [] PDF;
199  delete [] CPDF;
200  PDF = 0;
201  }
202 }

Member Function Documentation

G4double G4GoudsmitSaundersonTable::SampleTheta ( G4double  lambda,
G4double  Chia2,
G4double  rndm 
)

Definition at line 205 of file G4GoudsmitSaundersonTable.cc.

References test::a, test::b, G4Log(), and rndm().

206 {
207  //Benedito's procedure
208  G4double A[11],ThisPDF[320],ThisCPDF[320];
209  G4double coeff,Ckj,CkjPlus1,CkPlus1j,CkPlus1jPlus1,a,b,mmm,F;
210  G4int Ind0,KIndex=0,JIndex=0,IIndex=0;
211 
212 
213  ///////////////////////////////////////////////////////////////////////////
214  // Find Lambda and Chia2 Index
215  for(G4int k=0;k<75;k++){if((lambda>=LAMBDAN[k])&&(lambda<LAMBDAN[k+1])){KIndex=k;break;}}
216  for(G4int j=0;j<11;j++){A[j]=scrA[11*KIndex+j];}
217  for(G4int j=0;j<10;j++){if((Chia2>=A[j])&&(Chia2<A[j+1])){JIndex=j;break;}}
218 
219  ///////////////////////////////////////////////////////////////////////////
220  // Calculate some necessary coefficients for PDF and CPDF interpolation
221  coeff=(G4Log(LAMBDAN[KIndex+1]/LAMBDAN[KIndex]))*(G4Log(A[JIndex+1]/A[JIndex]));
222  Ckj=(G4Log(LAMBDAN[KIndex+1]/lambda))*(G4Log(A[JIndex+1]/Chia2))/coeff;
223  CkjPlus1=(G4Log(LAMBDAN[KIndex+1]/lambda))*(G4Log(Chia2/A[JIndex]))/coeff;
224  CkPlus1j=(G4Log(lambda/LAMBDAN[KIndex]))*(G4Log(A[JIndex+1]/Chia2))/coeff;
225  CkPlus1jPlus1=(G4Log(lambda/LAMBDAN[KIndex]))*(G4Log(Chia2/A[JIndex]))/coeff;
226  ///////////////////////////////////////////////////////////////////////////
227  // Calculate Interpolated PDF and CPDF arrays
228  Ind0=320*(11*KIndex+JIndex);
229  for(G4int i=0 ; i<320 ;i++){
230  ThisPDF[i]=Ckj*PDF[Ind0]+CkjPlus1*PDF[Ind0+320]+CkPlus1j*PDF[Ind0+3520]+CkPlus1jPlus1*PDF[Ind0+3840];
231  ThisCPDF[i]=Ckj*CPDF[Ind0]+CkjPlus1*CPDF[Ind0+320]+CkPlus1j*CPDF[Ind0+3520]+CkPlus1jPlus1*CPDF[Ind0+3840];
232  // Find u Index using secant method
233  if((i!=0)&&((rndm>=ThisCPDF[i-1])&&(rndm<ThisCPDF[i]))) {IIndex=i-1;break;}
234  Ind0++;
235  }
236 
237  ///////////////////////////////////////////////////////////////////////////
238  //CPDF^-1(rndm)=x ==> CPDF(x)=rndm;
239  a=uvalues[IIndex];
240  b=uvalues[IIndex+1];
241 
242  do{
243  mmm=0.5*(a+b);
244  F=(ThisCPDF[IIndex]+(mmm-uvalues[IIndex])*ThisPDF[IIndex]
245  +((mmm-uvalues[IIndex])*(mmm-uvalues[IIndex])
246  *(ThisPDF[IIndex+1]-ThisPDF[IIndex]))
247  /(2.*(uvalues[IIndex+1]-uvalues[IIndex])))-rndm;
248  if(F>0.)b=mmm;
249  else a=mmm;
250  } while(std::fabs(b-a)>1.0e-9);
251 
252  return mmm;
253 }
int G4int
Definition: G4Types.hh:78
double precision function rndm(RDUMMY)
Definition: dpm25nulib.f:1460
G4double G4Log(G4double x)
Definition: G4Log.hh:227
double G4double
Definition: G4Types.hh:76

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