[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$
Revision 1.2  1999/09/29 09:24:28  fca
Introduction of the Copyright and cvs Log

*/

///////////////////////////////////////////////////////////////////////////
// Class AliMath
// Various mathematical tools which may be very convenient while
// performing physics analysis.
//
// Example : Probability of a Chi-squared value
// =========
//
// AliMath M;
// Float_t chi2=20;            // The chi-squared value
// Int_t ndf=12;               // The number of degrees of freedom
// Float_t p=M.Prob(chi2,ndf); // The probability that at least a Chi-squared
//                             // value of chi2 will be observed, even for a
//                             // correct model
//
//--- Author: Nick van Eijndhoven 14-nov-1998 UU-SAP Utrecht
///////////////////////////////////////////////////////////////////////////

#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$
959fbac5	18	Revision 1.2 1999/09/29 09:24:28 fca
	19	Introduction of the Copyright and cvs Log
	20
4c039060	21	*/
4c039060	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=sumexp(-x+alog(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+alog(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=texp((-zz)
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(2ndf-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	///////////////////////////////////////////////////////////////////////////