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 ///////////////////////////////////////////////////////////////////////////