ITS/AliITSstatistics2.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2007-2009, 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 /* $Id$ */
  17
  18 //////////////////////////////////////////////////////////////////////////
  19 //  Alice ITS class to help keep statistical information. Can also be   //
  20 // used to fit data to lines and other 2 dimentional sytistical         //
  21 // operations.                                                          //
  22 //                                                                      //
  23 // version: 0.0.0 Draft.                                                //
  24 // Date: April 18 1999                                                  //
  25 // By: Bjorn S. Nilsen                                                  //
  26 // Updated: 1.0.0, Date: September 6 2007, By: Bjorn S. Nilsen          //
  27 //                                                                      //
  28 //////////////////////////////////////////////////////////////////////////
  29 #include <stdio.h>     //  ios::fmtflags fmt used in PrintAscii
  30 #include "Riostream.h" // IO functions.
  31 #include "TMath.h"     // TMath::Sqrt() function used.
  32 #include "AliITSstatistics2.h" // Also defined TObject {base class}
  33
  34 ClassImp(AliITSstatistics2)
  35
  36 //
  37 AliITSstatistics2::AliITSstatistics2() :
  38 TObject(), // Base Class
  39 fN(-1),    // number of enetries -1 => Uninitilized
  40 fOrder(0), // maximum moment of distributions (^n)
  41 fX(0),     //[fOrder] array of sums of x^n
  42 fYx(0),    //[fOrder] array of sums of (xy)^n
  43 fY(0),     //[fOrder] array of sums of y^n
  44 fW(0)      //[fOrder] array of sums of w^n (weights)
  45 //,fDig(5)   // The number of significant digits to keep
  46 //,fOver(0)  //! In case of numerical precistion problems
  47 {
  48     // default constructor
  49     // Inputs:
  50     //   none.
  51     // Outputs:
  52     //   none.
  53     // Return:
  54     //   A default constructed AliITSstatistics class
  55
  56     return;
  57 }
  58 //______________________________________________________________________
  59 AliITSstatistics2::AliITSstatistics2(Int_t order) :
  60 TObject(),     // Base Class
  61 fN(0),         // number of enetries -1 => Uninitilized
  62 fOrder(order), // maximum moment of distributions (^n)
  63 fX(new Double_t[order]),    //[fOrder] array of sums of x^n
  64 fYx(new Double_t[order]),   //[fOrder] array of sums of (xy)^n
  65 fY(new Double_t[order]),    //[fOrder] array of sums of y^n
  66 fW(new Double_t[order])     //[fOrder] array of sums of w^n (weights)
  67 //,fDig(5)   // The number of significant digits to keep
  68 //,fOver(0)                   //! In case of numeerical precistion problems
  69 {   // constructor to maximum moment/order order
  70     // Inputs:
  71     //   Int_t order   The maximum moment of distributions {for example x^n}
  72     Int_t i;
  73
  74     for(i=0;i<order;i++) {
  75         fX[i]  = 0.0;
  76         fY[i]  = 0.0;
  77         fYx[i] = 0.0;
  78         fW[i]  = 0.0;
  79     } // end for i
  80     fN = 0;
  81     return;
  82 }
  83 //______________________________________________________________________
  84 AliITSstatistics2::~AliITSstatistics2(){
  85     // destructor
  86     // Inputs:
  87     //   none.
  88     // Outputs:
  89     //   none.
  90     // Return:
  91     //   none.
  92
  93     if(fX!=0)  delete[] fX;
  94     if(fY!=0)  delete[] fY;
  95     if(fYx!=0) delete[] fYx;
  96     if(fW!=0)  delete[] fW;
  97     // if(fOver!=0) delete fOver; fOver=0;
  98     // fDig=0;
  99     fX  = 0;
 100     fY  = 0;
 101     fYx = 0;
 102     fW  = 0;
 103     fN  = 0;
 104     fOrder = 0;
 105 }
 106 //_______________________________________________________________
 107 AliITSstatistics2& AliITSstatistics2::operator=(AliITSstatistics2 &source){
 108     // operator =
 109     // Inputs:
 110     //   AliITSstaticstics2 &source The source of this copy.
 111     // Outputs:
 112     //   none.
 113     // Return:
 114     //   A copy of the source class
 115
 116     if(this==&source) return *this;
 117     TObject::operator=(source);
 118     Reset(source.GetOrder());
 119     fN = source.GetN();
 120     fOrder=source.GetOrder();
 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     // this->fDig = source.fDig;
 128     // if(fOver!=0) this->fOver = new AliITSstatistics2(*(source.fOver));
 129     // else fOver=0;
 130     return *this;
 131 }
 132 //_______________________________________________________________
 133 AliITSstatistics2::AliITSstatistics2(AliITSstatistics2 &source):
 134 TObject(source),          // Base Class
 135 fN(source.GetN()),        // number of enetries -1 => Uninitilized
 136 fOrder(source.GetOrder()),// maximum moment of distributions (^n)
 137 fX(new Double_t[source.GetOrder()]),//[fOrder] array of sums of x^n
 138 fYx(new Double_t[source.GetOrder()]),//[fOrder] array of sums of (xy)^n
 139 fY(new Double_t[source.GetOrder()]),//[fOrder] array of sums of y^n
 140 fW(new Double_t[source.GetOrder()]) //[fOrder] array of sums of w^n (weights)
 141 //,fDig(source.fDig)                  // The number of significant digits to keep
 142 //,fOver(0)             //! In case of numerical precistion problems
 143 {
 144     // Copy constructor
 145     // Inputs:
 146     //   AliITSstatistics2 & source the source of this copy
 147     // Outputs:
 148     //   none.
 149     // Return:
 150     //   A copy of the source.
 151
 152     for(Int_t i=0;i<source.fOrder;i++){
 153         this->fX[i] = source.fX[i];
 154         this->fYx[i] = source.fYx[i];
 155         this->fY[i] = source.fY[i];
 156         this->fW[i] = source.fW[i];
 157     } // end for i
 158     //if(fOver!=0) this->fOver = new AliITSstatistics2(*(source.fOver));
 159     return;
 160 }
 161 //______________________________________________________________________
 162 void AliITSstatistics2::Reset(Int_t order){
 163     // Reset/zero all statistics variables statistics
 164     // Inputs:
 165     //   none.
 166     // Outputs:
 167     //   none.
 168     // Return:
 169     //   none.
 170     Int_t i;
 171
 172     for(i=0;i<fOrder;i++) {
 173         fX[i] = 0.0;
 174         fY[i] = 0.0;
 175         fYx[i] = 0.0;
 176         fW[i] = 0.0;
 177     } // end for i
 178     fN = 0;
 179     if(order<0) return; // just zero
 180     if(fX!=0)  delete[] fX;
 181     if(fY!=0)  delete[] fY;
 182     if(fYx!=0) delete[] fYx;
 183     if(fW!=0)  delete[] fW;
 184     fX  = 0;
 185     fY  = 0;
 186     fYx = 0;
 187     fW  = 0;
 188     fN  = 0;
 189     fOrder = 0;
 190     if(order==0) return;
 191     fOrder = order;
 192     fX =  new Double_t[fOrder];
 193     fY =  new Double_t[fOrder];
 194     fYx = new Double_t[fOrder];
 195     fW =  new Double_t[fOrder];
 196     //if(fOver!=0) delete fOver; fOver = 0;
 197     return;
 198 }
 199 //----------------------------------------------------------------------
 200 /*
 201 void SetSignificantDigits(Int_t d){
 202     // Sets the number of significant digits. If adding a value to
 203     // one of this class' arrays looses significance at the fDig
 204     // level, a new instance of this class is created to keep
 205     // signigicance at or better than fDig level. if fDig<0, then
 206     // this feature it disabled and significance can be lost.
 207     // Inputs:
 208     //    Int_t  d   The new significance level
 209     // Outputs:
 210     //    none.
 211     // Return:
 212     //    none.
 213
 214     fDig = d;
 215 }
 216  */
 217 //______________________________________________________________________
 218 void AliITSstatistics2::AddValue(Double_t y,Double_t x,Double_t w=1.0){
 219     // add next x,y pair to statistics
 220     // Inputs:
 221     //   Double_t  y    y value of pair
 222     //   Double_t  x    x value of pair
 223     //   Double_t  w    weight of pair
 224     // Outputs:
 225     //   none.
 226     // Return:
 227     //   none.
 228     Double_t xs=1.0,ys=1.0,yxs=1.0,ws=1.0;
 229     Int_t i;
 230     const Double_t kBig=1.0e+38;
 231
 232     if(y>kBig || x>kBig || w>kBig) return;
 233     /*  If problem with precision, then creat/fill fOver
 234         as a partical sum to be added to "this" later.
 235         if(????fDig){
 236         if(fOver==0){
 237             fOver = new AliITSstatistics2(fOrder);
 238         } // end if fOver==0
 239         fOver->AddValue(y,x,w);
 240         return;
 241         } // end if(???)
 242      */
 243     fN++;
 244     for(i=0;i<GetOrder();i++){
 245         xs  *= x;
 246         ys  *= y;
 247         yxs *= x*y;
 248         ws  *= w;
 249         fX[i]  += xs*w;
 250         fY[i]  += ys*w;
 251         fYx[i] += yxs*w;
 252         fW[i]  += ws;
 253     } // end for i
 254 }
 255 //______________________________________________________________________
 256 Double_t AliITSstatistics2::GetXNth(Int_t order)const{
 257     // This give the unbiased estimator for the RMS.
 258     // Inputs:
 259     //   Int_t order   the order of x^n value to be returned
 260     // Output:
 261     //   none.
 262     // Return:
 263     //   The value sum{x^n}.
 264     Double_t s;
 265
 266     if(GetWN(1)!=0.0 && order<=GetOrder()) s = GetXN(order)/GetWN(1);
 267     else {
 268         s = 0.0;
 269         Error("GetXNth","error fOrder=%d fN=%d fW[0]=%f\n",
 270               GetOrder(),GetN(),GetWN(1));
 271     } // end else
 272     return s;
 273 }
 274 //______________________________________________________________________
 275 Double_t AliITSstatistics2::GetYNth(Int_t order)const{
 276     // This give the unbiased estimator for the RMS.
 277     // Inputs:
 278     //   Int_t order   the order of y^n value to be returned
 279     // Outputs:
 280     //   none.
 281     // Return:
 282     //  The value sum{y^n}
 283     Double_t s;
 284
 285     if(GetWN(1)!=0.0&&order<=GetOrder()) s = GetYN(order)/GetWN(1);
 286     else {
 287         s = 0.0;
 288         Error("GetYNth","fOrder=%d fN=%d fW[0]=%f\n",
 289               GetOrder(),GetN(),GetWN(1));
 290     } // end else
 291     return s;
 292 }
 293 //______________________________________________________________________
 294 Double_t AliITSstatistics2::GetYXNth(Int_t order)const{
 295     // This give the unbiased estimator for the RMS.
 296     // Inputs:
 297     //   Int_t order   the order of (xy)^n value to be returned
 298     // Outputs:
 299     //   none.
 300     // Return:
 301     //  The value sum{(xy)^n}
 302     Double_t s;
 303
 304     if(GetWN(1)!=0.0&&order<=GetOrder()) s = GetYXN(order)/GetWN(1);
 305     else {
 306         s = 0.0;
 307         Error("GetYXNth","fOrder=%d fN=%d fW[0]=%f\n",
 308               GetOrder(),GetN(),GetWN(1));
 309     } // end else
 310     return s;
 311 }
 312 //______________________________________________________________________
 313 Double_t AliITSstatistics2::GetRMSX()const{
 314     // This give the unbiased estimator for the RMS.
 315     // Inputs:
 316     //   none.
 317     // Outputs:
 318     //   none.
 319     // Return:
 320     //  The rms value
 321     Double_t x,x2,w,ww,s;
 322
 323     x  = GetMeanX(); // first order
 324     x2 = GetXNth(2); // second order
 325     w  = GetWN(1);   // first order
 326     ww = GetWN(2);   // second order
 327
 328     if(w*w==ww) return (-1.0);
 329     s = (x2-x*x)*w*w/(w*w-ww);
 330     return TMath::Sqrt(s);
 331 }
 332 //______________________________________________________________________
 333 Double_t AliITSstatistics2::GetRMSY()const{
 334     // This give the unbiased estimator for the RMS.
 335     // Inputs:
 336     //   none.
 337     // Outputs:
 338     //   none.
 339     // Return:
 340     //  The rms value
 341     Double_t x,x2,w,ww,s;
 342
 343     x  = GetMeanY(); // first order
 344     x2 = GetYNth(2); // second order
 345     w  = GetWN(1);   // first order
 346     ww = GetWN(2);   // second order
 347
 348     if(w*w==ww) return (-1.0);
 349     s = (x2-x*x)*w*w/(w*w-ww);
 350     return TMath::Sqrt(s);
 351 }
 352 //______________________________________________________________________
 353 Double_t AliITSstatistics2::GetRMSYX()const{
 354     // This give the unbiased estimator for the RMS.
 355     // Inputs:
 356     //   none.
 357     // Outputs:
 358     //   none.
 359     // Return:
 360     //  The rms value
 361     Double_t x,x2,w,ww,s;
 362
 363     x  = GetMeanYX(); // first order
 364     x2 = GetYXNth(2); // second order
 365     w  = GetWN(1);   // first order
 366     ww = GetWN(2);     // second order
 367
 368     if(w*w==ww) return (-1.0);
 369     s = (x2-x*x)*w*w/(w*w-ww);
 370     return TMath::Sqrt(s);
 371 }
 372 //______________________________________________________________________
 373 Double_t AliITSstatistics2::GetErrorMeanY()const{
 374     //This is the error in the mean or the square root of the
 375     // variance of the mean.
 376     // Inputs:
 377     //   none.
 378     // Outputs:
 379     //   none.
 380     // Return:
 381     //  The error on the mean
 382     Double_t rms,w,ww,s;
 383
 384     rms = GetRMSY();
 385     w   = GetWN(1);
 386     ww  = GetWN(2);
 387     s   = rms*rms*ww/(w*w);
 388     return TMath::Sqrt(s);
 389 }
 390 //______________________________________________________________________
 391 Double_t AliITSstatistics2::GetErrorMeanX()const{
 392     //This is the error in the mean or the square root of the
 393     // variance of the mean.
 394     // Inputs:
 395     //   none.
 396     // Outputs:
 397     //   none.
 398     // Return:
 399     //  The error on the mean
 400     Double_t rms,w,ww,s;
 401
 402     rms = GetRMSX();
 403     w   = GetWN(1);
 404     ww  = GetWN(2);
 405     s   = rms*rms*ww/(w*w);
 406     return TMath::Sqrt(s);
 407 }
 408 //______________________________________________________________________
 409 Double_t AliITSstatistics2::GetErrorMeanYX()const{
 410     //This is the error in the mean or the square root of the
 411     // variance of the mean.
 412     // Inputs:
 413     //   none.
 414     // Outputs:
 415     //   none.
 416     // Return:
 417     //  The error on the mean
 418     Double_t rms,w,ww,s;
 419
 420     rms = GetRMSYX();
 421     w   = GetWN(1);
 422     ww  = GetWN(2);
 423     s   = rms*rms*ww/(w*w);
 424     return TMath::Sqrt(s);
 425 }
 426 //______________________________________________________________________
 427 Double_t AliITSstatistics2::GetErrorRMSY()const{
 428     // This is the error in the mean or the square root of the variance
 429     // of the mean. at this moment this routine is only defined for
 430     // weights=1.
 431     // Inputs:
 432     //   none.
 433     // Outputs:
 434     //   none.
 435     // Return:
 436     //  The error on the rms
 437     Double_t x,x2,x3,x4,w,ww,m2,m4,n,s;
 438
 439     if(GetWN(1)!=(Double_t)GetN()||GetN()<4) return (-1.);
 440     x  = GetMeanY(); // first order
 441     x2 = GetYNth(2); // second order
 442     w  = GetWN(1);   // first order
 443     ww = GetWN(2);   // second order
 444     if(w*w==ww) return (-1.0);
 445     s = (x2-x*x)*w*w/(w*w-ww);
 446
 447     m2  = s;
 448     n   = (Double_t) GetN();
 449     x3  = GetYNth(3);
 450     x4  = GetYNth(4);
 451     // This equation assumes that all of the weights are equal to 1.
 452     m4  = (n/(n-1.))*(x4-3.*x*x3+6.*x*x*x2-2.*x*x*x*x);
 453     s   = (m4-(n-3.)*m2*m2/(n-1.))/n;
 454     return TMath::Sqrt(s);
 455 }
 456 //______________________________________________________________________
 457 Double_t AliITSstatistics2::GetErrorRMSX()const{
 458     // This is the error in the mean or the square root of the variance
 459     // of the mean. at this moment this routine is only defined for
 460     // weights=1.
 461     // Inputs:
 462     //   none.
 463     // Outputs:
 464     //   none.
 465     // Return:
 466     //  The error on the rms
 467     Double_t x,x2,x3,x4,w,ww,m2,m4,n,s;
 468
 469     if(GetWN(1)!=(Double_t)GetN()||GetN()<4) return (-1.);
 470     x  = GetMeanX(); // first order
 471     x2 = GetXNth(2); // second order
 472     w  = GetWN(1);   // first order
 473     ww = GetWN(2);   // second order
 474     if(w*w==ww) return (-1.0);
 475     s = (x2-x*x)*w*w/(w*w-ww);
 476
 477     m2  = s;
 478     n   = (Double_t) GetN();
 479     x3  = GetXNth(3);
 480     x4  = GetXNth(4);
 481 // This equation assumes that all of the weights are equal to 1.
 482     m4  = (n/(n-1.))*(x4-3.*x*x3+6.*x*x*x2-2.*x*x*x*x);
 483     s   = (m4-(n-3.)*m2*m2/(n-1.))/n;
 484     return TMath::Sqrt(s);
 485 }
 486 //______________________________________________________________________
 487 Double_t AliITSstatistics2::GetErrorRMSYX()const{
 488     // This is the error in the mean or the square root of the variance
 489     // of the mean. at this moment this routine is only defined for
 490     // weights=1.
 491     // Inputs:
 492     //   none.
 493     // Outputs:
 494     //   none.
 495     // Return:
 496     //  The error on the rms
 497     Double_t x,x2,x3,x4,w,ww,m2,m4,n,s;
 498
 499     if(GetWN(1)!=(Double_t)GetN()||GetN()<4) return (-1.);
 500     x  = GetMeanYX(); // first order
 501     x2 = GetYXNth(2); // second order
 502     w  = GetWN(1);    // first order
 503     ww = GetWN(2);    // second order
 504     if(w*w==ww) return (-1.0);
 505     s = (x2-x*x)*w*w/(w*w-ww);
 506
 507     m2  = s;
 508     n   = (Double_t) GetN();
 509     x3  = GetYXNth(3);
 510     x4  = GetYXNth(4);
 511     // This equation assumes that all of the weights are equal to 1.
 512     m4  = (n/(n-1.))*(x4-3.*x*x3+6.*x*x*x2-2.*x*x*x*x);
 513     s   = (m4-(n-3.)*m2*m2/(n-1.))/n;
 514     return TMath::Sqrt(s);
 515 }
 516 //_______________________________________________________________________
 517 Double_t AliITSstatistics2::FitToLine(Double_t &a,Double_t &ea,
 518                                       Double_t &b,Double_t &eb)const{
 519     // fit to y = ax+b returns Chi squared or -1.0 if an error.
 520     // The fitting is done by analitically minimizing
 521     /*
 522       Begin_Latex
 523       \begin{equation*}
 524       \Chi^{2}=\sum_{i} (y_{i}-a x_{i} -b)^{2} w_{i}
 525       \end{equation*}
 526       Where if the weight used in
 527       AliITSstatistics2::AddValue(Double_t y,Double_t x,Double_t w=1.0)
 528       is of the form
 529       \begin{equation*}
 530       w_{i}=\frac{1}{\delta y^{2}}.
 531       \end{equation*}
 532       Then we get the typicall chi square minimization.
 533       End_Latex
 534      */
 535     // Inputs:
 536     //   none.
 537     // Outputs:
 538     //   Double_t  a   The slope parameter
 539     //   Double_t  ea  Error on fitted slope parameter
 540     //   Double_t  b   The intercept paramter
 541     //   Double_t  eb  Error on fitted intercept parameter
 542     // Return:
 543     //  The Chi^2 of the fit
 544     Double_t c,d,e,f,g,h;
 545
 546     a = ea = b = eb = 0.0;
 547     if(GetOrder()<2 || GetN()<3){
 548         Error("FitToLine","Order=%d<2 or N=%d<3",GetOrder(),GetN());
 549         return -1.0;
 550     } // end if
 551     c = GetWN(1);
 552     d = GetYN(1);
 553     e = GetXN(1);
 554     f = GetYXN(1);
 555     g = GetXN(2);
 556     h = c*g-e*e;
 557     b = d*g-f*e;
 558     a = c*f-d*e;
 559     if(h==0.0){
 560         Error("FitToLine","vertical line: fOrder=%d fN=%d "
 561               "GetWN(1)=%g X GetXN(2)=%g - GetXN(1)=%g^2 = 0",
 562               GetOrder(),GetN(),c,g,e);
 563         return -1.0;
 564     } // end if h
 565     a = a/h;
 566     b = b/h;
 567     // Now for the errors.
 568     ea = c*c*g+(a*a-1.0)*c*e*e;
 569     ea = ea/(h*h);
 570     if(ea<0.0){
 571       Error("FitToLine","ea=%g is less than zero",ea);
 572       return -2.0;
 573     } // end if ea<0
 574     ea = TMath::Sqrt(ea);
 575     eb = c*g*g-2.0*d*e*g-2.0*(1.0-b)*c*e*e*g+2.0*(1.0-b)*d*e*e*e+
 576           GetYN(2)*e*e+(1.0-b)*(1.0-b)*c*e*e*e*e;
 577     eb = eb/(h*h);
 578     if(eb<0.0){
 579       Error("FitToLine","eb=%g is less than zero",eb);
 580       return -2.0;
 581     } // end if ea<0
 582     eb = TMath::Sqrt(eb);
 583     c = GetChiSquared(a,b);
 584     if(c<=0.0){ // must be a numerical precision problem.
 585     } // end if
 586     return c;
 587 }
 588 //_______________________________________________________________________
 589 Double_t AliITSstatistics2::GetChiSquared(Double_t a,Double_t b)const{
 590     //  Returns Chi^2 value of data to line y=ax+b with given a,b.
 591     /*
 592       Begin_Latex
 593       Note: The Chi^2 value is computed from the expression
 594       \begin{equation*}
 595       \chi^{2}=\sum_{i}{w_{i}y_{i}^{2}} + b^{2}\sum_{i}{w_{i}}
 596                 -2b\sum_{i}{w_{i}y_{i}}-2a\sum_{i}{w_{i}y_{i}x_{i}}
 597                 +2ab\sum_{i}{w_{i}x_{i}}
 598                 +a^{2}\sum_{i}w_{i}x_{i}^{2}
 599       \end{equation*}
 600       and not form the expression
 601       \begin{equation*}
 602       \chi^{2}= \sum_{i}{(y_{i}-ax_{i}-b)^{2}w_{i}.
 603       \end{equation*}
 604       Consiquently, there are occations when numerically these
 605       two expressions will not agree. In fact the form code here
 606       can give negitive values. This happens when the numerical
 607       significance is larger than the $\chi^{2}$ value. This should
 608       not be confused the the error values which can be returned.
 609       At present there is no check on the numberical significance
 610       of any results.
 611       End_Latex
 612      */
 613     // Inputs:
 614     //   Double_t  a   The slope parameter
 615     //   Double_t  b   The intercept paramter
 616     // Outputs::
 617     //   none.
 618     // Return:
 619     //  The Chi^2 of the fit
 620     Double_t c2;
 621
 622     c2 = GetYN(2)+b*b*GetWN(1)+
 623         a*a*GetXN(2)-2.0*b*GetYN(1)-2.0*a*GetYXN(1)+2.0*b*a*GetXN(1);
 624     c2 /= (Double_t)GetN() - 2.0;
 625     return c2;
 626 }
 627 //______________________________________________________________________
 628 void AliITSstatistics2::PrintAscii(ostream *os)const{
 629     // Print out class data values in Ascii Form to output stream
 630     // Inputs:
 631     //   ostream *os   Output stream where Ascii data is to be writen
 632     // Outputs:
 633     //   none.
 634     // Return:
 635     //   none.
 636     Int_t i;
 637 #if defined __GNUC__
 638 #if __GNUC__ > 2
 639     ios::fmtflags fmt;  // Standard IO format object, required for output.
 640 #else
 641     Int_t fmt;
 642 #endif
 643 #else
 644 #if defined __ICC || defined __ECC || defined __xlC__
 645     ios::fmtflags fmt;
 646 #else
 647     Int_t fmt;
 648 #endif
 649 #endif
 650
 651     *os << fN <<" "<< GetOrder();
 652     fmt = os->setf(ios::scientific); // set scientific floating point output
 653     for(i=0;i<GetOrder();i++) *os <<" "<< GetXN(i+1);
 654     for(i=0;i<GetOrder();i++) *os <<" "<< GetYXN(i+1);
 655     for(i=0;i<GetOrder();i++) *os <<" "<< GetYN(i+1);
 656     for(i=0;i<GetOrder();i++) *os <<" "<< GetWN(i+1);
 657     //if(fOver!=0) { *os << " " << fDig;
 658     //*os << " " << fOver;
 659     // } else *os << " " << fDig;
 660     os->flags(fmt); // reset back to old Formating.
 661     return;
 662 }
 663 //______________________________________________________________________
 664 void AliITSstatistics2::ReadAscii(istream *is){
 665     // Read in class data values in Ascii Form to output stream
 666     // Inputs:
 667     //   istream *is   Input stream where Ascii data is to be read in from
 668     // Outputs:
 669     //   none.
 670     // Return:
 671     //   none.
 672     Int_t i;
 673
 674     *is >> i >> fOrder;
 675     Reset(fOrder);
 676     fN = i;
 677     for(i=0;i<fOrder;i++) *is >> fX[i];
 678     for(i=0;i<fOrder;i++) *is >> fYx[i];
 679     for(i=0;i<fOrder;i++) *is >> fY[i];
 680     for(i=0;i<fOrder;i++) *is >> fW[i];
 681     //*is >> fDig;
 682     // if(fDig>0) *is >> fOver;
 683     // else fDig *= -1;
 684 }
 685 //______________________________________________________________________
 686 ostream &operator<<(ostream &os,const AliITSstatistics2 &s){
 687     // Standard output streaming function
 688     // Inputs:
 689     //   ostream            &os  output steam
 690     //   AliITSstatistics2 &s class to be streamed.
 691     // Output:
 692     //   none.
 693     // Return:
 694     //   ostream &os  The stream pointer
 695
 696     s.PrintAscii(&os);
 697     return os;
 698 }
 699 //______________________________________________________________________
 700 istream &operator>>(istream &is,AliITSstatistics2 &s){
 701     // Standard inputput streaming function
 702     // Inputs:
 703     //   istream            &is  input steam
 704     //   AliITSstatistics2 &s class to be streamed.
 705     // Output:
 706     //   none.
 707     // Return:
 708     //   ostream &os  The stream pointer
 709
 710     s.ReadAscii(&is);
 711     return is;
 712 }
 713