[u/mrichter/AliRoot.git] / RALICE / AliMath.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/*
$Log$
*/

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