204{
205
206 G4double A[11],ThisPDF[320],ThisCPDF[320];
207 G4double coeff,Ckj,CkjPlus1,CkPlus1j,CkPlus1jPlus1,a,b,mmm,F;
208 G4int Ind0,KIndex=0,JIndex=0,IIndex=0;
209
210
211
212
213 for(
G4int k=0;k<75;k++){
if((lambda>=LAMBDAN[k])&&(lambda<LAMBDAN[k+1])){KIndex=k;
break;}}
214 for(
G4int j=0;j<11;j++){A[j]=scrA[11*KIndex+j];}
215 for(
G4int j=0;j<10;j++){
if((Chia2>=A[j])&&(Chia2<A[j+1])){JIndex=j;
break;}}
216
217
218
219 coeff=(std::log(LAMBDAN[KIndex+1]/LAMBDAN[KIndex]))*(std::log(A[JIndex+1]/A[JIndex]));
220 Ckj=(std::log(LAMBDAN[KIndex+1]/lambda))*(std::log(A[JIndex+1]/Chia2))/coeff;
221 CkjPlus1=(std::log(LAMBDAN[KIndex+1]/lambda))*(std::log(Chia2/A[JIndex]))/coeff;
222 CkPlus1j=(std::log(lambda/LAMBDAN[KIndex]))*(std::log(A[JIndex+1]/Chia2))/coeff;
223 CkPlus1jPlus1=(std::log(lambda/LAMBDAN[KIndex]))*(std::log(Chia2/A[JIndex]))/coeff;
224
225
226 Ind0=320*(11*KIndex+JIndex);
227 for(
G4int i=0 ; i<320 ;i++){
228 ThisPDF[i]=Ckj*PDF[Ind0]+CkjPlus1*PDF[Ind0+320]+CkPlus1j*PDF[Ind0+3520]+CkPlus1jPlus1*PDF[Ind0+3840];
229 ThisCPDF[i]=Ckj*CPDF[Ind0]+CkjPlus1*CPDF[Ind0+320]+CkPlus1j*CPDF[Ind0+3520]+CkPlus1jPlus1*CPDF[Ind0+3840];
230
231 if((i!=0)&&((rndm>=ThisCPDF[i-1])&&(rndm<ThisCPDF[i]))) {IIndex=i-1;break;}
232 Ind0++;
233 }
234
235
236
237 a=uvalues[IIndex];
238 b=uvalues[IIndex+1];
239
240 do{
241 mmm=0.5*(a+b);
242 F=(ThisCPDF[IIndex]+(mmm-uvalues[IIndex])*ThisPDF[IIndex]
243 +((mmm-uvalues[IIndex])*(mmm-uvalues[IIndex])
244 *(ThisPDF[IIndex+1]-ThisPDF[IIndex]))
245 /(2.*(uvalues[IIndex+1]-uvalues[IIndex])))-rndm;
246 if(F>0.)b=mmm;
247 else a=mmm;
248 } while(std::fabs(b-a)>1.0e-9);
249
250 return mmm;
251}