]>
Commit | Line | Data |
---|---|---|
4c039060 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | /* | |
17 | $Log$ | |
959fbac5 | 18 | Revision 1.2 1999/09/29 09:24:28 fca |
19 | Introduction of the Copyright and cvs Log | |
20 | ||
4c039060 | 21 | */ |
22 | ||
959fbac5 | 23 | /////////////////////////////////////////////////////////////////////////// |
24 | // Class AliMath | |
25 | // Various mathematical tools which may be very convenient while | |
26 | // performing physics analysis. | |
27 | // | |
28 | // Example : Probability of a Chi-squared value | |
29 | // ========= | |
30 | // | |
31 | // AliMath M; | |
32 | // Float_t chi2=20; // The chi-squared value | |
33 | // Int_t ndf=12; // The number of degrees of freedom | |
34 | // Float_t p=M.Prob(chi2,ndf); // The probability that at least a Chi-squared | |
35 | // // value of chi2 will be observed, even for a | |
36 | // // correct model | |
37 | // | |
38 | //--- Author: Nick van Eijndhoven 14-nov-1998 UU-SAP Utrecht | |
39 | /////////////////////////////////////////////////////////////////////////// | |
40 | ||
d88f97cc | 41 | #include "AliMath.h" |
42 | ||
43 | ClassImp(AliMath) // Class implementation to enable ROOT I/O | |
44 | ||
45 | AliMath::AliMath() | |
46 | { | |
47 | // Default constructor | |
48 | } | |
49 | /////////////////////////////////////////////////////////////////////////// | |
50 | AliMath::~AliMath() | |
51 | { | |
52 | // Destructor | |
53 | } | |
54 | /////////////////////////////////////////////////////////////////////////// | |
55 | Float_t AliMath::Gamma(Float_t z) | |
56 | { | |
57 | // Computation of gamma(z) for all z>0. | |
58 | // | |
59 | // The algorithm is based on the article by C.Lanczos [1] as denoted in | |
60 | // Numerical Recipes 2nd ed. on p. 207 (W.H.Press et al.). | |
61 | // | |
62 | // [1] C.Lanczos, SIAM Journal of Numerical Analysis B1 (1964), 86. | |
63 | // | |
64 | //--- Nve 14-nov-1998 UU-SAP Utrecht | |
65 | ||
66 | if (z<=0.) | |
67 | { | |
68 | cout << "*Gamma(z)* Wrong argument z = " << z << endl; | |
69 | return 0; | |
70 | } | |
71 | ||
72 | Float_t v=LnGamma(z); | |
73 | return exp(v); | |
74 | } | |
75 | /////////////////////////////////////////////////////////////////////////// | |
76 | Float_t AliMath::Gamma(Float_t a,Float_t x) | |
77 | { | |
78 | // Computation of the incomplete gamma function P(a,x) | |
79 | // | |
80 | // The algorithm is based on the formulas and code as denoted in | |
81 | // Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.). | |
82 | // | |
83 | //--- Nve 14-nov-1998 UU-SAP Utrecht | |
84 | ||
85 | if (a<=0.) | |
86 | { | |
87 | cout << "*Gamma(a,x)* Invalid argument a = " << a << endl; | |
88 | return 0; | |
89 | } | |
90 | ||
91 | if (x<=0.) | |
92 | { | |
93 | if (x<0) cout << "*Gamma(a,x)* Invalid argument x = " << x << endl; | |
94 | return 0; | |
95 | } | |
96 | ||
97 | if (x<(a+1.)) | |
98 | { | |
99 | return GamSer(a,x); | |
100 | } | |
101 | else | |
102 | { | |
103 | return GamCf(a,x); | |
104 | } | |
105 | } | |
106 | /////////////////////////////////////////////////////////////////////////// | |
107 | Float_t AliMath::LnGamma(Float_t z) | |
108 | { | |
109 | // Computation of ln[gamma(z)] for all z>0. | |
110 | // | |
111 | // The algorithm is based on the article by C.Lanczos [1] as denoted in | |
112 | // Numerical Recipes 2nd ed. on p. 207 (W.H.Press et al.). | |
113 | // | |
114 | // [1] C.Lanczos, SIAM Journal of Numerical Analysis B1 (1964), 86. | |
115 | // | |
116 | // The accuracy of the result is better than 2e-10. | |
117 | // | |
118 | //--- Nve 14-nov-1998 UU-SAP Utrecht | |
119 | ||
120 | if (z<=0.) | |
121 | { | |
122 | cout << "*LnGamma(z)* Wrong argument z = " << z << endl; | |
123 | return 0; | |
124 | } | |
125 | ||
126 | // Coefficients for the series expansion | |
127 | Double_t c[7]; | |
128 | c[0]= 2.5066282746310005; | |
129 | c[1]= 76.18009172947146; | |
130 | c[2]=-86.50532032941677; | |
131 | c[3]= 24.01409824083091; | |
132 | c[4]= -1.231739572450155; | |
133 | c[5]= 0.1208650973866179e-2; | |
134 | c[6]= -0.5395239384953e-5; | |
135 | ||
136 | Double_t x=z; | |
137 | Double_t y=x; | |
138 | Double_t tmp=x+5.5; | |
139 | tmp=(x+0.5)*log(tmp)-tmp; | |
140 | Double_t ser=1.000000000190015; | |
141 | for (Int_t i=1; i<7; i++) | |
142 | { | |
143 | y+=1.; | |
144 | ser+=c[i]/y; | |
145 | } | |
146 | Float_t v=tmp+log(c[0]*ser/x); | |
147 | return v; | |
148 | } | |
149 | /////////////////////////////////////////////////////////////////////////// | |
150 | Float_t AliMath::GamSer(Float_t a,Float_t x) | |
151 | { | |
152 | // Computation of the incomplete gamma function P(a,x) | |
153 | // via its series representation. | |
154 | // | |
155 | // The algorithm is based on the formulas and code as denoted in | |
156 | // Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.). | |
157 | // | |
158 | //--- Nve 14-nov-1998 UU-SAP Utrecht | |
159 | ||
160 | Int_t itmax=100; // Maximum number of iterations | |
161 | Float_t eps=3.e-7; // Relative accuracy | |
162 | ||
163 | if (a<=0.) | |
164 | { | |
165 | cout << "*GamSer(a,x)* Invalid argument a = " << a << endl; | |
166 | return 0; | |
167 | } | |
168 | ||
169 | if (x<=0.) | |
170 | { | |
171 | if (x<0) cout << "*GamSer(a,x)* Invalid argument x = " << x << endl; | |
172 | return 0; | |
173 | } | |
174 | ||
175 | Float_t gln=LnGamma(a); | |
176 | Float_t ap=a; | |
177 | Float_t sum=1./a; | |
178 | Float_t del=sum; | |
179 | for (Int_t n=1; n<=itmax; n++) | |
180 | { | |
181 | ap+=1.; | |
182 | del=del*x/ap; | |
183 | sum+=del; | |
184 | if (fabs(del)<fabs(sum*eps)) break; | |
185 | if (n==itmax) cout << "*GamSer(a,x)* a too large or itmax too small" << endl; | |
186 | } | |
187 | Float_t v=sum*exp(-x+a*log(x)-gln); | |
188 | return v; | |
189 | } | |
190 | /////////////////////////////////////////////////////////////////////////// | |
191 | Float_t AliMath::GamCf(Float_t a,Float_t x) | |
192 | { | |
193 | // Computation of the incomplete gamma function P(a,x) | |
194 | // via its continued fraction representation. | |
195 | // | |
196 | // The algorithm is based on the formulas and code as denoted in | |
197 | // Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.). | |
198 | // | |
199 | //--- Nve 14-nov-1998 UU-SAP Utrecht | |
200 | ||
201 | Int_t itmax=100; // Maximum number of iterations | |
202 | Float_t eps=3.e-7; // Relative accuracy | |
203 | Float_t fpmin=1.e-30; // Smallest Float_t value allowed here | |
204 | ||
205 | if (a<=0.) | |
206 | { | |
207 | cout << "*GamCf(a,x)* Invalid argument a = " << a << endl; | |
208 | return 0; | |
209 | } | |
210 | ||
211 | if (x<=0.) | |
212 | { | |
213 | if (x<0) cout << "*GamCf(a,x)* Invalid argument x = " << x << endl; | |
214 | return 0; | |
215 | } | |
216 | ||
217 | Float_t gln=LnGamma(a); | |
218 | Float_t b=x+1.-a; | |
219 | Float_t c=1./fpmin; | |
220 | Float_t d=1./b; | |
221 | Float_t h=d; | |
222 | Float_t an,del; | |
223 | for (Int_t i=1; i<=itmax; i++) | |
224 | { | |
225 | an=float(-i)*(float(i)-a); | |
226 | b+=2.; | |
227 | d=an*d+b; | |
228 | if (fabs(d)<fpmin) d=fpmin; | |
229 | c=b+an/c; | |
230 | if (fabs(c)<fpmin) c=fpmin; | |
231 | d=1./d; | |
232 | del=d*c; | |
233 | h=h*del; | |
234 | if (fabs(del-1.)<eps) break; | |
235 | if (i==itmax) cout << "*GamCf(a,x)* a too large or itmax too small" << endl; | |
236 | } | |
237 | Float_t v=exp(-x+a*log(x)-gln)*h; | |
238 | return (1.-v); | |
239 | } | |
240 | /////////////////////////////////////////////////////////////////////////// | |
241 | Float_t AliMath::Erf(Float_t x) | |
242 | { | |
243 | // Computation of the error function erf(x). | |
244 | // | |
245 | //--- NvE 14-nov-1998 UU-SAP Utrecht | |
246 | ||
247 | return (1.-Erfc(x)); | |
248 | } | |
249 | /////////////////////////////////////////////////////////////////////////// | |
250 | Float_t AliMath::Erfc(Float_t x) | |
251 | { | |
252 | // Computation of the complementary error function erfc(x). | |
253 | // | |
254 | // The algorithm is based on a Chebyshev fit as denoted in | |
255 | // Numerical Recipes 2nd ed. on p. 214 (W.H.Press et al.). | |
256 | // | |
257 | // The fractional error is always less than 1.2e-7. | |
258 | // | |
259 | //--- Nve 14-nov-1998 UU-SAP Utrecht | |
260 | ||
261 | // The parameters of the Chebyshev fit | |
262 | const Float_t a1=-1.26551223, a2=1.00002368, | |
263 | a3= 0.37409196, a4=0.09678418, | |
264 | a5=-0.18628806, a6=0.27886807, | |
265 | a7=-1.13520398, a8=1.48851587, | |
266 | a9=-0.82215223, a10=0.17087277; | |
267 | ||
268 | Float_t v=1.; // The return value | |
269 | ||
270 | Float_t z=fabs(x); | |
271 | ||
272 | if (z <= 0.) return v; // erfc(0)=1 | |
273 | ||
274 | Float_t t=1./(1.+0.5*z); | |
275 | ||
276 | v=t*exp((-z*z) | |
277 | +a1+t*(a2+t*(a3+t*(a4+t*(a5+t*(a6+t*(a7+t*(a8+t*(a9+t*a10))))))))); | |
278 | ||
279 | if (x < 0.) v=2.-v; // erfc(-x)=2-erfc(x) | |
280 | ||
281 | return v; | |
282 | } | |
283 | /////////////////////////////////////////////////////////////////////////// | |
284 | Float_t AliMath::Prob(Float_t chi2,Int_t ndf) | |
285 | { | |
286 | // Computation of the probability for a certain Chi-squared (chi2) | |
287 | // and number of degrees of freedom (ndf). | |
288 | // | |
289 | // Calculations are based on the incomplete gamma function P(a,x), | |
290 | // where a=ndf/2 and x=chi2/2. | |
291 | // | |
292 | // P(a,x) represents the probability that the observed Chi-squared | |
293 | // for a correct model should be less than the value chi2. | |
294 | // | |
295 | // The returned probability corresponds to 1-P(a,x), | |
296 | // which denotes the probability that an observed Chi-squared exceeds | |
297 | // the value chi2 by chance, even for a correct model. | |
298 | // | |
299 | //--- NvE 14-nov-1998 UU-SAP Utrecht | |
300 | ||
301 | if (ndf <= 0) return 0; // Set CL to zero in case ndf<=0 | |
302 | ||
303 | if (chi2 <= 0.) | |
304 | { | |
305 | if (chi2 < 0.) | |
306 | { | |
307 | return 0; | |
308 | } | |
309 | else | |
310 | { | |
311 | return 1; | |
312 | } | |
313 | } | |
314 | ||
315 | // Alternative which is exact | |
316 | // This code may be activated in case the gamma function gives problems | |
317 | // if (ndf==1) | |
318 | // { | |
319 | // Float_t v=1.-Erf(sqrt(chi2)/sqrt(2.)); | |
320 | // return v; | |
321 | // } | |
322 | ||
323 | // Gaussian approximation for large ndf | |
324 | // This code may be activated in case the gamma function shows a problem | |
325 | // Float_t q=sqrt(2.*chi2)-sqrt(float(2*ndf-1)); | |
326 | // if (n>30 && q>0.) | |
327 | // { | |
328 | // Float_t v=0.5*(1.-Erf(q/sqrt(2.))); | |
329 | // return v; | |
330 | // } | |
331 | ||
332 | // Evaluate the incomplete gamma function | |
333 | Float_t a=float(ndf)/2.; | |
334 | Float_t x=chi2/2.; | |
335 | return (1.-Gamma(a,x)); | |
336 | } | |
337 | /////////////////////////////////////////////////////////////////////////// |