ITS/AliITSstatistics2.cxx

   1 //////////////////////////////////////////////////////////////////////////
   2 //  Alice ITS class to help keep statistical information                //
   3 //                                                                      //
   4 // version: 0.0.0 Draft.                                                //
   5 // Date: April 18 1999                                                  //
   6 // By: Bjorn S. Nilsen                                                  //
   7 //                                                                      //
   8 //////////////////////////////////////////////////////////////////////////
   9 #include <stdio.h>
  10 #include <math.h>
  11 #include "TMath.h"
  12 #include "AliITSstatistics2.h"
  13 #include "Riostream.h"
  14
  15 ClassImp(AliITSstatistics2)
  16
  17 //
  18 AliITSstatistics2::AliITSstatistics2() :
  19 TObject(), // Base Class
  20 fN(-1),    // number of enetries -1 => Uninitilized
  21 fOrder(0), // maximum moment of distributions (^n)
  22 fX(0),    //[fOrder] array of sums of x^n
  23 fYx(0),   //[fOrder] array of sums of (xy)^n
  24 fY(0),    //[fOrder] array of sums of y^n
  25 fW(0){    //[fOrder] array of sums of w^n (weights)
  26     // default constructor
  27     // Inputs:
  28     //   none.
  29     // Outputs:
  30     //   none.
  31     // Return:
  32     //   A default constructed AliITSstatistics class
  33
  34     return;
  35 }
  36 //______________________________________________________________________
  37 AliITSstatistics2::AliITSstatistics2(Int_t order) :
  38 TObject(),     // Base Class
  39 fN(0),         // number of enetries -1 => Uninitilized
  40 fOrder(order), // maximum moment of distributions (^n)
  41 fX(new Double_t[order]),    //[fOrder] array of sums of x^n
  42 fYx(new Double_t[order]),   //[fOrder] array of sums of (xy)^n
  43 fY(new Double_t[order]),    //[fOrder] array of sums of y^n
  44 fW(new Double_t[order]){    //[fOrder] array of sums of w^n (weights)
  45     // constructor to maximum moment/order order
  46     // Inputs:
  47     //   Int_t order   The maximum moment of distributions {for example x^n}
  48     Int_t i;
  49
  50     for(i=0;i<order;i++) {
  51         fX[i]  = 0.0;
  52         fY[i]  = 0.0;
  53         fYx[i] = 0.0;
  54         fW[i]  = 0.0;
  55     } // end for i
  56     fN = 0;
  57     return;
  58 }
  59 //______________________________________________________________________
  60 AliITSstatistics2::~AliITSstatistics2(){
  61     // destructor
  62     // Inputs:
  63     //   none.
  64     // Outputs:
  65     //   none.
  66     // Return:
  67     //   none.
  68
  69     if(fX!=0)  delete[] fX;
  70     if(fY!=0)  delete[] fY;
  71     if(fYx!=0) delete[] fYx;
  72     if(fW!=0)  delete[] fW;
  73     fX  = 0;
  74     fY  = 0;
  75     fYx = 0;
  76     fW  = 0;
  77     fN  = 0;
  78     fOrder = 0;
  79 }
  80 //_______________________________________________________________
  81 AliITSstatistics2& AliITSstatistics2::operator=(AliITSstatistics2 &source){
  82     // operator =
  83     // Inputs:
  84     //   AliITSstaticstics2 &source The source of this copy.
  85     // Outputs:
  86     //   none.
  87     // Return:
  88     //   A copy of the source class
  89
  90     if(this==&source) return *this;
  91     TObject::operator=(source);
  92     Reset(source.GetOrder());
  93     fN = source.GetN();
  94     fOrder=source.GetOrder();
  95     for(Int_t i=0;i<source.fOrder;i++){
  96         this->fX[i] = source.fX[i];
  97         this->fYx[i] = source.fYx[i];
  98         this->fY[i] = source.fY[i];
  99         this->fW[i] = source.fW[i];
 100     } // end for i
 101     return *this;
 102 }
 103 //_______________________________________________________________
 104 AliITSstatistics2::AliITSstatistics2(AliITSstatistics2 &source):
 105 TObject(source),          // Base Class
 106 fN(source.GetN()),        // number of enetries -1 => Uninitilized
 107 fOrder(source.GetOrder()),// maximum moment of distributions (^n)
 108 fX(new Double_t[source.GetOrder()]),//[fOrder] array of sums of x^n
 109 fYx(new Double_t[source.GetOrder()]),//[fOrder] array of sums of (xy)^n
 110 fY(new Double_t[source.GetOrder()]),//[fOrder] array of sums of y^n
 111 fW(new Double_t[source.GetOrder()]){//[fOrder] array of sums of w^n (weights)
 112     // Copy constructor
 113     // Inputs:
 114     //   AliITSstatistics2 & source the source of this copy
 115     // Outputs:
 116     //   none.
 117     // Return:
 118     //   A copy of the source.
 119
 120     if(this==&source) return;
 121     for(Int_t i=0;i<source.fOrder;i++){
 122         this->fX[i] = source.fX[i];
 123         this->fYx[i] = source.fYx[i];
 124         this->fY[i] = source.fY[i];
 125         this->fW[i] = source.fW[i];
 126     } // end for i
 127 }
 128 //______________________________________________________________________
 129 void AliITSstatistics2::Reset(Int_t order){
 130     // Reset/zero all statistics variables statistics
 131     // Inputs:
 132     //   none.
 133     // Outputs:
 134     //   none.
 135     // Return:
 136     //   none.
 137     Int_t i;
 138
 139     for(i=0;i<fOrder;i++) {
 140         fX[i] = 0.0;
 141         fY[i] = 0.0;
 142         fYx[i] = 0.0;
 143         fW[i] = 0.0;
 144     } // end for i
 145     fN = 0;
 146     if(order<0) return; // just zero
 147     if(fX!=0)  delete[] fX;
 148     if(fY!=0)  delete[] fY;
 149     if(fYx!=0) delete[] fYx;
 150     if(fW!=0)  delete[] fW;
 151     fX  = 0;
 152     fY  = 0;
 153     fYx = 0;
 154     fW  = 0;
 155     fN  = 0;
 156     fOrder = 0;
 157     if(order==0) return;
 158     fOrder = order;
 159     fX =  new Double_t[fOrder];
 160     fY =  new Double_t[fOrder];
 161     fYx = new Double_t[fOrder];
 162     fW =  new Double_t[fOrder];
 163     return;
 164 }
 165 //______________________________________________________________________
 166 void AliITSstatistics2::AddValue(Double_t y,Double_t x,Double_t w=1.0){
 167     // add next x,y pair to statistics
 168     // Inputs:
 169     //   Double_t  y    y value of pair
 170     //   Double_t  x    x value of pair
 171     //   Double_t  w    weight of pair
 172     // Outputs:
 173     //   none.
 174     // Return:
 175     //   none.
 176     Double_t xs=1.0,ys=1.0,yxs=1.0,ws=1.0;
 177     Int_t i;
 178     const Double_t kBig=1.0e+38;
 179
 180     if(y>kBig || x>kBig || w>kBig) return;
 181     fN++;
 182     for(i=0;i<fOrder;i++){
 183         xs  *= x;
 184         ys  *= y;
 185         yxs *= x*y;
 186         ws  *= w;
 187         fX[i]  += xs*w;
 188         fY[i]  += ys*w;
 189         fYx[i] += yxs*w;
 190         fW[i]  += ws;
 191     } // end for i
 192 }
 193 //______________________________________________________________________
 194 Double_t AliITSstatistics2::GetXNth(Int_t order)const{
 195     // This give the unbiased estimator for the RMS.
 196     // Inputs:
 197     //   Int_t order   the order of x^n value to be returned
 198     // Output:
 199     //   none.
 200     // Return:
 201     //   The value sum{x^n}.
 202     Double_t s;
 203
 204     if(fW[0]!=0.0&&order<=fOrder) s = fX[order-1]/fW[0];
 205     else {
 206         s = 0.0;
 207         Error("GetXNth","error fOrder=%d fN=%d fW[0]=%f\n",
 208               fOrder,fN,fW[0]);
 209     } // end else
 210     return s;
 211 }
 212 //______________________________________________________________________
 213 Double_t AliITSstatistics2::GetYNth(Int_t order)const{
 214     // This give the unbiased estimator for the RMS.
 215     // Inputs:
 216     //   Int_t order   the order of y^n value to be returned
 217     // Outputs:
 218     //   none.
 219     // Return:
 220     //  The value sum{y^n}
 221     Double_t s;
 222
 223     if(fW[0]!=0.0&&order<=fOrder) s = fY[order-1]/fW[0];
 224     else {
 225         s = 0.0;
 226         Error("GetYNth","fOrder=%d fN=%d fW[0]=%f\n",
 227               fOrder,fN,fW[0]);
 228     } // end else
 229     return s;
 230 }
 231 //______________________________________________________________________
 232 Double_t AliITSstatistics2::GetYXNth(Int_t order)const{
 233     // This give the unbiased estimator for the RMS.
 234     // Inputs:
 235     //   Int_t order   the order of (xy)^n value to be returned
 236     // Outputs:
 237     //   none.
 238     // Return:
 239     //  The value sum{(xy)^n}
 240     Double_t s;
 241
 242     if(fW[0]!=0.0&&order<=fOrder) s = fYx[order-1]/fW[0];
 243     else {
 244         s = 0.0;
 245         Error("GetYXNth","fOrder=%d fN=%d fW[0]=%f\n",
 246                fOrder,fN,fW[0]);
 247     } // end else
 248     return s;
 249 }
 250 //______________________________________________________________________
 251 Double_t AliITSstatistics2::GetRMSX()const{
 252     // This give the unbiased estimator for the RMS.
 253     // Inputs:
 254     //   none.
 255     // Outputs:
 256     //   none.
 257     // Return:
 258     //  The rms value
 259     Double_t x,x2,w,ww,s;
 260
 261     x  = GetMeanX(); // first order
 262     x2 = GetXNth(2); // second order
 263     w  = fW[0];     // first order - 1.
 264     ww = fW[1];     // second order - 1.
 265
 266     if(w*w==ww) return (-1.0);
 267     s = (x2-x*x)*w*w/(w*w-ww);
 268     return TMath::Sqrt(s);
 269 }
 270 //______________________________________________________________________
 271 Double_t AliITSstatistics2::GetRMSY()const{
 272     // This give the unbiased estimator for the RMS.
 273     // Inputs:
 274     //   none.
 275     // Outputs:
 276     //   none.
 277     // Return:
 278     //  The rms value
 279     Double_t x,x2,w,ww,s;
 280
 281     x  = GetMeanY(); // first order
 282     x2 = GetYNth(2); // second order
 283     w  = fW[0];     // first order - 1.
 284     ww = fW[1];     // second order - 1.
 285
 286     if(w*w==ww) return (-1.0);
 287     s = (x2-x*x)*w*w/(w*w-ww);
 288     return TMath::Sqrt(s);
 289 }
 290 //______________________________________________________________________
 291 Double_t AliITSstatistics2::GetRMSYX()const{
 292     // This give the unbiased estimator for the RMS.
 293     // Inputs:
 294     //   none.
 295     // Outputs:
 296     //   none.
 297     // Return:
 298     //  The rms value
 299     Double_t x,x2,w,ww,s;
 300
 301     x  = GetMeanYX(); // first order
 302     x2 = GetYXNth(2); // second order
 303     w  = fW[0];     // first order - 1.
 304     ww = fW[1];     // second order - 1.
 305
 306     if(w*w==ww) return (-1.0);
 307     s = (x2-x*x)*w*w/(w*w-ww);
 308     return TMath::Sqrt(s);
 309 }
 310 //______________________________________________________________________
 311 Double_t AliITSstatistics2::GetErrorMeanY()const{
 312     //This is the error in the mean or the square root of the
 313     // variance of the mean.
 314     // Inputs:
 315     //   none.
 316     // Outputs:
 317     //   none.
 318     // Return:
 319     //  The error on the mean
 320     Double_t rms,w,ww,s;
 321
 322     rms = GetRMSY();
 323     w   = fW[0];
 324     ww  = fW[1];
 325     s   = rms*rms*ww/(w*w);
 326     return TMath::Sqrt(s);
 327 }
 328 //______________________________________________________________________
 329 Double_t AliITSstatistics2::GetErrorMeanX()const{
 330     //This is the error in the mean or the square root of the
 331     // variance of the mean.
 332     // Inputs:
 333     //   none.
 334     // Outputs:
 335     //   none.
 336     // Return:
 337     //  The error on the mean
 338     Double_t rms,w,ww,s;
 339
 340     rms = GetRMSX();
 341     w   = fW[0];
 342     ww  = fW[1];
 343     s   = rms*rms*ww/(w*w);
 344     return TMath::Sqrt(s);
 345 }
 346 //______________________________________________________________________
 347 Double_t AliITSstatistics2::GetErrorMeanYX()const{
 348     //This is the error in the mean or the square root of the
 349     // variance of the mean.
 350     // Inputs:
 351     //   none.
 352     // Outputs:
 353     //   none.
 354     // Return:
 355     //  The error on the mean
 356     Double_t rms,w,ww,s;
 357
 358     rms = GetRMSYX();
 359     w   = fW[0];
 360     ww  = fW[1];
 361     s   = rms*rms*ww/(w*w);
 362     return TMath::Sqrt(s);
 363 }
 364 //______________________________________________________________________
 365 Double_t AliITSstatistics2::GetErrorRMSY()const{
 366     // This is the error in the mean or the square root of the variance
 367     // of the mean. at this moment this routine is only defined for
 368     // weights=1.
 369     // Inputs:
 370     //   none.
 371     // Outputs:
 372     //   none.
 373     // Return:
 374     //  The error on the rms
 375     Double_t x,x2,x3,x4,w,ww,m2,m4,n,s;
 376
 377     if(fW[0]!=(Double_t)fN||GetN()<4) return (-1.);
 378     x  = GetMeanY(); // first order
 379     x2 = GetYNth(2); // second order
 380     w  = fW[0];     // first order - 1.
 381     ww = fW[1];     // second order - 1.
 382     if(w*w==ww) return (-1.0);
 383     s = (x2-x*x)*w*w/(w*w-ww);
 384
 385     m2  = s;
 386     n   = (Double_t) GetN();
 387     x3  = GetYNth(3);
 388     x4  = GetYNth(4);
 389     // This equation assumes that all of the weights are equal to 1.
 390     m4  = (n/(n-1.))*(x4-3.*x*x3+6.*x*x*x2-2.*x*x*x*x);
 391     s   = (m4-(n-3.)*m2*m2/(n-1.))/n;
 392     return TMath::Sqrt(s);
 393 }
 394 //______________________________________________________________________
 395 Double_t AliITSstatistics2::GetErrorRMSX()const{
 396     // This is the error in the mean or the square root of the variance
 397     // of the mean. at this moment this routine is only defined for
 398     // weights=1.
 399     // Inputs:
 400     //   none.
 401     // Outputs:
 402     //   none.
 403     // Return:
 404     //  The error on the rms
 405     Double_t x,x2,x3,x4,w,ww,m2,m4,n,s;
 406
 407     if(fW[0]!=(Double_t)fN||GetN()<4) return (-1.);
 408     x  = GetMeanX(); // first order
 409     x2 = GetXNth(2); // second order
 410     w  = fW[0];     // first order - 1.
 411     ww = fW[1];     // second order - 1.
 412     if(w*w==ww) return (-1.0);
 413     s = (x2-x*x)*w*w/(w*w-ww);
 414
 415     m2  = s;
 416     n   = (Double_t) GetN();
 417     x3  = GetXNth(3);
 418     x4  = GetXNth(4);
 419 // This equation assumes that all of the weights are equal to 1.
 420     m4  = (n/(n-1.))*(x4-3.*x*x3+6.*x*x*x2-2.*x*x*x*x);
 421     s   = (m4-(n-3.)*m2*m2/(n-1.))/n;
 422     return TMath::Sqrt(s);
 423 }
 424 //______________________________________________________________________
 425 Double_t AliITSstatistics2::GetErrorRMSYX()const{
 426     // This is the error in the mean or the square root of the variance
 427     // of the mean. at this moment this routine is only defined for
 428     // weights=1.
 429     // Inputs:
 430     //   none.
 431     // Outputs:
 432     //   none.
 433     // Return:
 434     //  The error on the rms
 435     Double_t x,x2,x3,x4,w,ww,m2,m4,n,s;
 436
 437     if(fW[0]!=(Double_t)fN||GetN()<4) return (-1.);
 438     x  = GetMeanYX(); // first order
 439     x2 = GetYXNth(2); // second order
 440     w  = fW[0];     // first order - 1.
 441     ww = fW[1];     // second order - 1.
 442     if(w*w==ww) return (-1.0);
 443     s = (x2-x*x)*w*w/(w*w-ww);
 444
 445     m2  = s;
 446     n   = (Double_t) GetN();
 447     x3  = GetYXNth(3);
 448     x4  = GetYXNth(4);
 449     // This equation assumes that all of the weights are equal to 1.
 450     m4  = (n/(n-1.))*(x4-3.*x*x3+6.*x*x*x2-2.*x*x*x*x);
 451     s   = (m4-(n-3.)*m2*m2/(n-1.))/n;
 452     return TMath::Sqrt(s);
 453 }
 454 //_______________________________________________________________________
 455 Double_t AliITSstatistics2::FitToLine(Double_t &a,Double_t &b)const{
 456     // fit to y = ax+b returns Chi squared or -1.0 if an error
 457     // Inputs:
 458     //   none.
 459     // Outputs:
 460     //   Double_t  a   The slope parameter
 461     //   Double_t  b   The intercept paramter
 462     // Return:
 463     //  The Chi^2 of the fit
 464     Double_t c,d,e,f,g,h;
 465
 466     a = b = 0.0;
 467     if(fOrder<2 || fN<3){
 468         Error("FitToLine","Order=%d<2 or N=%d<3",fOrder,fN);
 469         return -1.0;
 470     } // end if
 471     c = GetWN(1);
 472     d = GetYN(1);
 473     e = GetXN(1);
 474     f = GetYXN(1);
 475     g = GetXN(2);
 476     h = c*g-e*e;
 477     b = d*g-f*e;
 478     a = c*f-d*e;
 479     if(h!=0.0){
 480         a = a/h;
 481         b = b/h;
 482     }else{
 483         Error("FitToLine","vertical line: fOrder=%d fN=%d "
 484               "GetWN(1)=%g X GetXN(2)=%g - GetXN(1)=%g^2 = 0",fOrder,fN,c,g,e);
 485         return -1.0;
 486     } // end if h
 487     return GetChiSquared(a,b);
 488 }
 489 //_______________________________________________________________________
 490 Double_t AliITSstatistics2::GetChiSquared(Double_t a,Double_t b)const{
 491     //  returns Chi^2 value of data to line y=ax+b with given a,b
 492     // Inputs:
 493     //   Double_t  a   The slope parameter
 494     //   Double_t  b   The intercept paramter
 495     // Outputs::
 496     //   none.
 497     // Return:
 498     //  The Chi^2 of the fit
 499     Double_t c2;
 500
 501     c2 = GetYN(2)+b*b*GetWN(1)+
 502         a*a*GetXN(2)-2.0*b*GetYN(1)-2.0*a*GetYXN(1)+2.0*b*a*GetXN(1);
 503     c2 /= (Double_t)fN - 2.0;
 504     return c2;
 505 }
 506 //______________________________________________________________________
 507 void AliITSstatistics2::PrintAscii(ostream *os)const{
 508     // Print out class data values in Ascii Form to output stream
 509     // Inputs:
 510     //   ostream *os   Output stream where Ascii data is to be writen
 511     // Outputs:
 512     //   none.
 513     // Return:
 514     //   none.
 515     Int_t i;
 516 #if defined __GNUC__
 517 #if __GNUC__ > 2
 518     ios::fmtflags fmt;
 519 #else
 520     Int_t fmt;
 521 #endif
 522 #else
 523 #if defined __ICC || defined __ECC || defined __xlC__
 524     ios::fmtflags fmt;
 525 #else
 526     Int_t fmt;
 527 #endif
 528 #endif
 529
 530     *os << fN <<" "<< fOrder;
 531     fmt = os->setf(ios::scientific); // set scientific floating point output
 532     for(i=0;i<fOrder;i++) *os <<" "<< fX[i];
 533     for(i=0;i<fOrder;i++) *os <<" "<< fYx[i];
 534     for(i=0;i<fOrder;i++) *os <<" "<< fY[i];
 535     for(i=0;i<fOrder;i++) *os <<" "<< fW[i];
 536     os->flags(fmt); // reset back to old Formating.
 537     return;
 538 }
 539 //______________________________________________________________________
 540 void AliITSstatistics2::ReadAscii(istream *is){
 541     // Read in class data values in Ascii Form to output stream
 542     // Inputs:
 543     //   istream *is   Input stream where Ascii data is to be read in from
 544     // Outputs:
 545     //   none.
 546     // Return:
 547     //   none.
 548     Int_t i;
 549
 550     *is >> i >> fOrder;
 551     Reset(fOrder);
 552     fN = i;
 553     for(i=0;i<fOrder;i++) *is >> fX[i];
 554     for(i=0;i<fOrder;i++) *is >> fYx[i];
 555     for(i=0;i<fOrder;i++) *is >> fY[i];
 556     for(i=0;i<fOrder;i++) *is >> fW[i];
 557 }
 558 //______________________________________________________________________
 559 ostream &operator<<(ostream &os,const AliITSstatistics2 &s){
 560     // Standard output streaming function
 561     // Inputs:
 562     //   ostream            &os  output steam
 563     //   AliITSstatistics2 &s class to be streamed.
 564     // Output:
 565     //   none.
 566     // Return:
 567     //   ostream &os  The stream pointer
 568
 569     s.PrintAscii(&os);
 570     return os;
 571 }
 572 //______________________________________________________________________
 573 istream &operator>>(istream &is,AliITSstatistics2 &s){
 574     // Standard inputput streaming function
 575     // Inputs:
 576     //   istream            &is  input steam
 577     //   AliITSstatistics2 &s class to be streamed.
 578     // Output:
 579     //   none.
 580     // Return:
 581     //   ostream &os  The stream pointer
 582
 583     s.ReadAscii(&is);
 584     return is;
 585 }