3 ClassImp(AliMath) // Class implementation to enable ROOT I/O
9 ///////////////////////////////////////////////////////////////////////////
14 ///////////////////////////////////////////////////////////////////////////
15 Float_t AliMath::Gamma(Float_t z)
17 // Computation of gamma(z) for all z>0.
19 // The algorithm is based on the article by C.Lanczos [1] as denoted in
20 // Numerical Recipes 2nd ed. on p. 207 (W.H.Press et al.).
22 // [1] C.Lanczos, SIAM Journal of Numerical Analysis B1 (1964), 86.
24 //--- Nve 14-nov-1998 UU-SAP Utrecht
28 cout << "*Gamma(z)* Wrong argument z = " << z << endl;
35 ///////////////////////////////////////////////////////////////////////////
36 Float_t AliMath::Gamma(Float_t a,Float_t x)
38 // Computation of the incomplete gamma function P(a,x)
40 // The algorithm is based on the formulas and code as denoted in
41 // Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.).
43 //--- Nve 14-nov-1998 UU-SAP Utrecht
47 cout << "*Gamma(a,x)* Invalid argument a = " << a << endl;
53 if (x<0) cout << "*Gamma(a,x)* Invalid argument x = " << x << endl;
66 ///////////////////////////////////////////////////////////////////////////
67 Float_t AliMath::LnGamma(Float_t z)
69 // Computation of ln[gamma(z)] for all z>0.
71 // The algorithm is based on the article by C.Lanczos [1] as denoted in
72 // Numerical Recipes 2nd ed. on p. 207 (W.H.Press et al.).
74 // [1] C.Lanczos, SIAM Journal of Numerical Analysis B1 (1964), 86.
76 // The accuracy of the result is better than 2e-10.
78 //--- Nve 14-nov-1998 UU-SAP Utrecht
82 cout << "*LnGamma(z)* Wrong argument z = " << z << endl;
86 // Coefficients for the series expansion
88 c[0]= 2.5066282746310005;
89 c[1]= 76.18009172947146;
90 c[2]=-86.50532032941677;
91 c[3]= 24.01409824083091;
92 c[4]= -1.231739572450155;
93 c[5]= 0.1208650973866179e-2;
94 c[6]= -0.5395239384953e-5;
99 tmp=(x+0.5)*log(tmp)-tmp;
100 Double_t ser=1.000000000190015;
101 for (Int_t i=1; i<7; i++)
106 Float_t v=tmp+log(c[0]*ser/x);
109 ///////////////////////////////////////////////////////////////////////////
110 Float_t AliMath::GamSer(Float_t a,Float_t x)
112 // Computation of the incomplete gamma function P(a,x)
113 // via its series representation.
115 // The algorithm is based on the formulas and code as denoted in
116 // Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.).
118 //--- Nve 14-nov-1998 UU-SAP Utrecht
120 Int_t itmax=100; // Maximum number of iterations
121 Float_t eps=3.e-7; // Relative accuracy
125 cout << "*GamSer(a,x)* Invalid argument a = " << a << endl;
131 if (x<0) cout << "*GamSer(a,x)* Invalid argument x = " << x << endl;
135 Float_t gln=LnGamma(a);
139 for (Int_t n=1; n<=itmax; n++)
144 if (fabs(del)<fabs(sum*eps)) break;
145 if (n==itmax) cout << "*GamSer(a,x)* a too large or itmax too small" << endl;
147 Float_t v=sum*exp(-x+a*log(x)-gln);
150 ///////////////////////////////////////////////////////////////////////////
151 Float_t AliMath::GamCf(Float_t a,Float_t x)
153 // Computation of the incomplete gamma function P(a,x)
154 // via its continued fraction representation.
156 // The algorithm is based on the formulas and code as denoted in
157 // Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.).
159 //--- Nve 14-nov-1998 UU-SAP Utrecht
161 Int_t itmax=100; // Maximum number of iterations
162 Float_t eps=3.e-7; // Relative accuracy
163 Float_t fpmin=1.e-30; // Smallest Float_t value allowed here
167 cout << "*GamCf(a,x)* Invalid argument a = " << a << endl;
173 if (x<0) cout << "*GamCf(a,x)* Invalid argument x = " << x << endl;
177 Float_t gln=LnGamma(a);
183 for (Int_t i=1; i<=itmax; i++)
185 an=float(-i)*(float(i)-a);
188 if (fabs(d)<fpmin) d=fpmin;
190 if (fabs(c)<fpmin) c=fpmin;
194 if (fabs(del-1.)<eps) break;
195 if (i==itmax) cout << "*GamCf(a,x)* a too large or itmax too small" << endl;
197 Float_t v=exp(-x+a*log(x)-gln)*h;
200 ///////////////////////////////////////////////////////////////////////////
201 Float_t AliMath::Erf(Float_t x)
203 // Computation of the error function erf(x).
205 //--- NvE 14-nov-1998 UU-SAP Utrecht
209 ///////////////////////////////////////////////////////////////////////////
210 Float_t AliMath::Erfc(Float_t x)
212 // Computation of the complementary error function erfc(x).
214 // The algorithm is based on a Chebyshev fit as denoted in
215 // Numerical Recipes 2nd ed. on p. 214 (W.H.Press et al.).
217 // The fractional error is always less than 1.2e-7.
219 //--- Nve 14-nov-1998 UU-SAP Utrecht
221 // The parameters of the Chebyshev fit
222 const Float_t a1=-1.26551223, a2=1.00002368,
223 a3= 0.37409196, a4=0.09678418,
224 a5=-0.18628806, a6=0.27886807,
225 a7=-1.13520398, a8=1.48851587,
226 a9=-0.82215223, a10=0.17087277;
228 Float_t v=1.; // The return value
232 if (z <= 0.) return v; // erfc(0)=1
234 Float_t t=1./(1.+0.5*z);
237 +a1+t*(a2+t*(a3+t*(a4+t*(a5+t*(a6+t*(a7+t*(a8+t*(a9+t*a10)))))))));
239 if (x < 0.) v=2.-v; // erfc(-x)=2-erfc(x)
243 ///////////////////////////////////////////////////////////////////////////
244 Float_t AliMath::Prob(Float_t chi2,Int_t ndf)
246 // Computation of the probability for a certain Chi-squared (chi2)
247 // and number of degrees of freedom (ndf).
249 // Calculations are based on the incomplete gamma function P(a,x),
250 // where a=ndf/2 and x=chi2/2.
252 // P(a,x) represents the probability that the observed Chi-squared
253 // for a correct model should be less than the value chi2.
255 // The returned probability corresponds to 1-P(a,x),
256 // which denotes the probability that an observed Chi-squared exceeds
257 // the value chi2 by chance, even for a correct model.
259 //--- NvE 14-nov-1998 UU-SAP Utrecht
261 if (ndf <= 0) return 0; // Set CL to zero in case ndf<=0
275 // Alternative which is exact
276 // This code may be activated in case the gamma function gives problems
279 // Float_t v=1.-Erf(sqrt(chi2)/sqrt(2.));
283 // Gaussian approximation for large ndf
284 // This code may be activated in case the gamma function shows a problem
285 // Float_t q=sqrt(2.*chi2)-sqrt(float(2*ndf-1));
288 // Float_t v=0.5*(1.-Erf(q/sqrt(2.)));
292 // Evaluate the incomplete gamma function
293 Float_t a=float(ndf)/2.;
295 return (1.-Gamma(a,x));
297 ///////////////////////////////////////////////////////////////////////////