RALICE/AliMath.cxx

   1 #include "AliMath.h"
   2
   3 ClassImp(AliMath) // Class implementation to enable ROOT I/O
   4
   5 AliMath::AliMath()
   6 {
   7 // Default constructor
   8 }
   9 ///////////////////////////////////////////////////////////////////////////
  10 AliMath::~AliMath()
  11 {
  12 // Destructor
  13 }
  14 ///////////////////////////////////////////////////////////////////////////
  15 Float_t AliMath::Gamma(Float_t z)
  16 {
  17 // Computation of gamma(z) for all z>0.
  18 //
  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.).
  21 //
  22 // [1] C.Lanczos, SIAM Journal of Numerical Analysis B1 (1964), 86.
  23 //
  24 //--- Nve 14-nov-1998 UU-SAP Utrecht
  25
  26  if (z<=0.)
  27  {
  28   cout << "*Gamma(z)* Wrong argument z = " << z << endl;
  29   return 0;
  30  }
  31
  32  Float_t v=LnGamma(z);
  33  return exp(v);
  34 }
  35 ///////////////////////////////////////////////////////////////////////////
  36 Float_t AliMath::Gamma(Float_t a,Float_t x)
  37 {
  38 // Computation of the incomplete gamma function P(a,x)
  39 //
  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.).
  42 //
  43 //--- Nve 14-nov-1998 UU-SAP Utrecht
  44
  45  if (a<=0.)
  46  {
  47   cout << "*Gamma(a,x)* Invalid argument a = " << a << endl;
  48   return 0;
  49  }
  50
  51  if (x<=0.)
  52  {
  53   if (x<0) cout << "*Gamma(a,x)* Invalid argument x = " << x << endl;
  54   return 0;
  55  }
  56
  57  if (x<(a+1.))
  58  {
  59   return GamSer(a,x);
  60  }
  61  else
  62  {
  63   return GamCf(a,x);
  64  }
  65 }
  66 ///////////////////////////////////////////////////////////////////////////
  67 Float_t AliMath::LnGamma(Float_t z)
  68 {
  69 // Computation of ln[gamma(z)] for all z>0.
  70 //
  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.).
  73 //
  74 // [1] C.Lanczos, SIAM Journal of Numerical Analysis B1 (1964), 86.
  75 //
  76 // The accuracy of the result is better than 2e-10.
  77 //
  78 //--- Nve 14-nov-1998 UU-SAP Utrecht
  79
  80  if (z<=0.)
  81  {
  82   cout << "*LnGamma(z)* Wrong argument z = " << z << endl;
  83   return 0;
  84  }
  85
  86  // Coefficients for the series expansion
  87  Double_t c[7];
  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;
  95
  96  Double_t x=z;
  97  Double_t y=x;
  98  Double_t tmp=x+5.5;
  99  tmp=(x+0.5)*log(tmp)-tmp;
 100  Double_t ser=1.000000000190015;
 101  for (Int_t i=1; i<7; i++)
 102  {
 103   y+=1.;
 104   ser+=c[i]/y;
 105  }
 106  Float_t v=tmp+log(c[0]*ser/x);
 107  return v;
 108 }
 109 ///////////////////////////////////////////////////////////////////////////
 110 Float_t AliMath::GamSer(Float_t a,Float_t x)
 111 {
 112 // Computation of the incomplete gamma function P(a,x)
 113 // via its series representation.
 114 //
 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.).
 117 //
 118 //--- Nve 14-nov-1998 UU-SAP Utrecht
 119
 120  Int_t itmax=100;   // Maximum number of iterations
 121  Float_t eps=3.e-7; // Relative accuracy
 122
 123  if (a<=0.)
 124  {
 125   cout << "*GamSer(a,x)* Invalid argument a = " << a << endl;
 126   return 0;
 127  }
 128
 129  if (x<=0.)
 130  {
 131   if (x<0) cout << "*GamSer(a,x)* Invalid argument x = " << x << endl;
 132   return 0;
 133  }
 134
 135  Float_t gln=LnGamma(a);
 136  Float_t ap=a;
 137  Float_t sum=1./a;
 138  Float_t del=sum;
 139  for (Int_t n=1; n<=itmax; n++)
 140  {
 141   ap+=1.;
 142   del=del*x/ap;
 143   sum+=del;
 144   if (fabs(del)<fabs(sum*eps)) break;
 145   if (n==itmax) cout << "*GamSer(a,x)* a too large or itmax too small" << endl;
 146  }
 147  Float_t v=sum*exp(-x+a*log(x)-gln);
 148  return v;
 149 }
 150 ///////////////////////////////////////////////////////////////////////////
 151 Float_t AliMath::GamCf(Float_t a,Float_t x)
 152 {
 153 // Computation of the incomplete gamma function P(a,x)
 154 // via its continued fraction representation.
 155 //
 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.).
 158 //
 159 //--- Nve 14-nov-1998 UU-SAP Utrecht
 160
 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
 164
 165  if (a<=0.)
 166  {
 167   cout << "*GamCf(a,x)* Invalid argument a = " << a << endl;
 168   return 0;
 169  }
 170
 171  if (x<=0.)
 172  {
 173   if (x<0) cout << "*GamCf(a,x)* Invalid argument x = " << x << endl;
 174   return 0;
 175  }
 176
 177  Float_t gln=LnGamma(a);
 178  Float_t b=x+1.-a;
 179  Float_t c=1./fpmin;
 180  Float_t d=1./b;
 181  Float_t h=d;
 182  Float_t an,del;
 183  for (Int_t i=1; i<=itmax; i++)
 184  {
 185   an=float(-i)*(float(i)-a);
 186   b+=2.;
 187   d=an*d+b;
 188   if (fabs(d)<fpmin) d=fpmin;
 189   c=b+an/c;
 190   if (fabs(c)<fpmin) c=fpmin;
 191   d=1./d;
 192   del=d*c;
 193   h=h*del;
 194   if (fabs(del-1.)<eps) break;
 195   if (i==itmax) cout << "*GamCf(a,x)* a too large or itmax too small" << endl;
 196  }
 197  Float_t v=exp(-x+a*log(x)-gln)*h;
 198  return (1.-v);
 199 }
 200 ///////////////////////////////////////////////////////////////////////////
 201 Float_t AliMath::Erf(Float_t x)
 202 {
 203 // Computation of the error function erf(x).
 204 //
 205 //--- NvE 14-nov-1998 UU-SAP Utrecht
 206
 207  return (1.-Erfc(x));
 208 }
 209 ///////////////////////////////////////////////////////////////////////////
 210 Float_t AliMath::Erfc(Float_t x)
 211 {
 212 // Computation of the complementary error function erfc(x).
 213 //
 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.).
 216 //
 217 // The fractional error is always less than 1.2e-7.
 218 //
 219 //--- Nve 14-nov-1998 UU-SAP Utrecht
 220
 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;
 227
 228  Float_t v=1.; // The return value
 229
 230  Float_t z=fabs(x);
 231
 232  if (z <= 0.) return v; // erfc(0)=1
 233
 234  Float_t t=1./(1.+0.5*z);
 235
 236  v=t*exp((-z*z)
 237          +a1+t*(a2+t*(a3+t*(a4+t*(a5+t*(a6+t*(a7+t*(a8+t*(a9+t*a10)))))))));
 238
 239  if (x < 0.) v=2.-v; // erfc(-x)=2-erfc(x)
 240
 241  return v;
 242 }
 243 ///////////////////////////////////////////////////////////////////////////
 244 Float_t AliMath::Prob(Float_t chi2,Int_t ndf)
 245 {
 246 // Computation of the probability for a certain Chi-squared (chi2)
 247 // and number of degrees of freedom (ndf).
 248 //
 249 // Calculations are based on the incomplete gamma function P(a,x),
 250 // where a=ndf/2 and x=chi2/2.
 251 //
 252 // P(a,x) represents the probability that the observed Chi-squared
 253 // for a correct model should be less than the value chi2.
 254 //
 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.
 258 //
 259 //--- NvE 14-nov-1998 UU-SAP Utrecht
 260
 261  if (ndf <= 0) return 0; // Set CL to zero in case ndf<=0
 262
 263  if (chi2 <= 0.)
 264  {
 265   if (chi2 < 0.)
 266   {
 267     return 0;
 268   }
 269   else
 270   {
 271    return 1;
 272   }
 273  }
 274
 275 // Alternative which is exact
 276 // This code may be activated in case the gamma function gives problems
 277 // if (ndf==1)
 278 // {
 279 //  Float_t v=1.-Erf(sqrt(chi2)/sqrt(2.));
 280 //  return v;
 281 // }
 282
 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));
 286 // if (n>30 && q>0.)
 287 // {
 288 //  Float_t v=0.5*(1.-Erf(q/sqrt(2.)));
 289 //  return v;
 290 // }
 291
 292  // Evaluate the incomplete gamma function
 293  Float_t a=float(ndf)/2.;
 294  Float_t x=chi2/2.;
 295  return (1.-Gamma(a,x));
 296 }
 297 ///////////////////////////////////////////////////////////////////////////