0826a0a3233b97befdbde623e89b532d921f1fff
[u/mrichter/AliRoot.git] / RALICE / AliMath.cxx
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$
18 Revision 1.2  1999/09/29 09:24:28  fca
19 Introduction of the Copyright and cvs Log
20
21 */
22
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
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 ///////////////////////////////////////////////////////////////////////////