[u/mrichter/AliRoot.git] / ITS / AliITSstatistics2.cxx

//////////////////////////////////////////////////////////////////////////
//  Alice ITS class to help keep statistical information                //
//                                                                      //
// version: 0.0.0 Draft.                                                //
// Date: April 18 1999                                                  //
// By: Bjorn S. Nilsen                                                  //
//                                                                      //
//////////////////////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#include "TMath.h"
#include "AliITSstatistics2.h"
#include "Riostream.h"

ClassImp(AliITSstatistics2)

//
AliITSstatistics2::AliITSstatistics2() : 
TObject(), // Base Class
fN(-1),    // number of enetries -1 => Uninitilized
fOrder(0), // maximum moment of distributions (^n)
fX(0),    //[fOrder] array of sums of x^n
fYx(0),   //[fOrder] array of sums of (xy)^n
fY(0),    //[fOrder] array of sums of y^n
fW(0){    //[fOrder] array of sums of w^n (weights)
    // default constructor
    // Inputs:
    //   none.
    // Outputs:
    //   none.
    // Return:
    //   A default constructed AliITSstatistics class

    return;
}
//______________________________________________________________________
AliITSstatistics2::AliITSstatistics2(Int_t order) : 
TObject(),     // Base Class
fN(0),         // number of enetries -1 => Uninitilized
fOrder(order), // maximum moment of distributions (^n)
fX(new Double_t[order]),    //[fOrder] array of sums of x^n
fYx(new Double_t[order]),   //[fOrder] array of sums of (xy)^n
fY(new Double_t[order]),    //[fOrder] array of sums of y^n
fW(new Double_t[order]){    //[fOrder] array of sums of w^n (weights)
    // constructor to maximum moment/order order
    // Inputs:
    //   Int_t order   The maximum moment of distributions {for example x^n}
    Int_t i;

    for(i=0;i<order;i++) {
        fX[i]  = 0.0;
        fY[i]  = 0.0;
        fYx[i] = 0.0;
        fW[i]  = 0.0;
    } // end for i
    fN = 0;
    return;
}
//______________________________________________________________________
AliITSstatistics2::~AliITSstatistics2(){
    // destructor
    // Inputs:
    //   none.
    // Outputs:
    //   none.
    // Return:
    //   none.

    if(fX!=0)  delete[] fX;
    if(fY!=0)  delete[] fY;
    if(fYx!=0) delete[] fYx;
    if(fW!=0)  delete[] fW;
    fX  = 0;
    fY  = 0;
    fYx = 0;
    fW  = 0;
    fN  = 0;
    fOrder = 0;
}
//_______________________________________________________________
AliITSstatistics2& AliITSstatistics2::operator=(AliITSstatistics2 &source){
    // operator =
    // Inputs:
    //   AliITSstaticstics2 &source The source of this copy.
    // Outputs:
    //   none.
    // Return:
    //   A copy of the source class

    if(this==&source) return *this;
    TObject::operator=(source);
    Reset(source.GetOrder());
    fN = source.GetN();
    fOrder=source.GetOrder();
    for(Int_t i=0;i<source.fOrder;i++){
        this->fX[i] = source.fX[i];
        this->fYx[i] = source.fYx[i];
        this->fY[i] = source.fY[i];
        this->fW[i] = source.fW[i];
    } // end for i
    return *this;
}
//_______________________________________________________________
AliITSstatistics2::AliITSstatistics2(AliITSstatistics2 &source): 
TObject(source),          // Base Class
fN(source.GetN()),        // number of enetries -1 => Uninitilized
fOrder(source.GetOrder()),// maximum moment of distributions (^n)
fX(new Double_t[source.GetOrder()]),//[fOrder] array of sums of x^n
fYx(new Double_t[source.GetOrder()]),//[fOrder] array of sums of (xy)^n
fY(new Double_t[source.GetOrder()]),//[fOrder] array of sums of y^n
fW(new Double_t[source.GetOrder()]){//[fOrder] array of sums of w^n (weights)
    // Copy constructor
    // Inputs:
    //   AliITSstatistics2 & source the source of this copy
    // Outputs:
    //   none.
    // Return:
    //   A copy of the source.

    if(this==&source) return;
    for(Int_t i=0;i<source.fOrder;i++){
        this->fX[i] = source.fX[i];
        this->fYx[i] = source.fYx[i];
        this->fY[i] = source.fY[i];
        this->fW[i] = source.fW[i];
    } // end for i
}
//______________________________________________________________________
void AliITSstatistics2::Reset(Int_t order){
    // Reset/zero all statistics variables statistics
    // Inputs:
    //   none.
    // Outputs:
    //   none.
    // Return:
    //   none.
    Int_t i;

    for(i=0;i<fOrder;i++) {
        fX[i] = 0.0;
        fY[i] = 0.0;
        fYx[i] = 0.0;
        fW[i] = 0.0;
    } // end for i
    fN = 0;
    if(order<0) return; // just zero
    if(fX!=0)  delete[] fX;
    if(fY!=0)  delete[] fY;
    if(fYx!=0) delete[] fYx;
    if(fW!=0)  delete[] fW;
    fX  = 0;
    fY  = 0;
    fYx = 0;
    fW  = 0;
    fN  = 0;
    fOrder = 0;
    if(order==0) return;
    fOrder = order;
    fX =  new Double_t[fOrder];
    fY =  new Double_t[fOrder];
    fYx = new Double_t[fOrder];
    fW =  new Double_t[fOrder];
    return;
}
//______________________________________________________________________
void AliITSstatistics2::AddValue(Double_t y,Double_t x,Double_t w=1.0){
    // add next x,y pair to statistics
    // Inputs:
    //   Double_t  y    y value of pair
    //   Double_t  x    x value of pair
    //   Double_t  w    weight of pair
    // Outputs:
    //   none.
    // Return:
    //   none.
    Double_t xs=1.0,ys=1.0,yxs=1.0,ws=1.0;
    Int_t i;
    const Double_t kBig=1.0e+38;

    if(y>kBig || x>kBig || w>kBig) return;
    fN++;
    for(i=0;i<fOrder;i++){
	xs  *= x;
	ys  *= y;
	yxs *= x*y;
	ws  *= w;
	fX[i]  += xs*w;
	fY[i]  += ys*w;
	fYx[i] += yxs*w;
	fW[i]  += ws;
    } // end for i
}
//______________________________________________________________________
Double_t AliITSstatistics2::GetXNth(Int_t order)const{
    // This give the unbiased estimator for the RMS.
    // Inputs:
    //   Int_t order   the order of x^n value to be returned
    // Output:
    //   none.
    // Return:
    //   The value sum{x^n}.
    Double_t s;

    if(fW[0]!=0.0&&order<=fOrder) s = fX[order-1]/fW[0];
    else {
	s = 0.0;
	Error("GetXNth","error fOrder=%d fN=%d fW[0]=%f\n",
              fOrder,fN,fW[0]);
    } // end else
    return s;
}
//______________________________________________________________________
Double_t AliITSstatistics2::GetYNth(Int_t order)const{
    // This give the unbiased estimator for the RMS.
    // Inputs:
    //   Int_t order   the order of y^n value to be returned
    // Outputs:
    //   none.
    // Return:
    //  The value sum{y^n}
    Double_t s;

    if(fW[0]!=0.0&&order<=fOrder) s = fY[order-1]/fW[0];
    else {
	s = 0.0;
	Error("GetYNth","fOrder=%d fN=%d fW[0]=%f\n",
              fOrder,fN,fW[0]);
    } // end else
    return s;
}
//______________________________________________________________________
Double_t AliITSstatistics2::GetYXNth(Int_t order)const{
    // This give the unbiased estimator for the RMS.
    // Inputs:
    //   Int_t order   the order of (xy)^n value to be returned
    // Outputs:
    //   none.
    // Return:
    //  The value sum{(xy)^n}
    Double_t s;

    if(fW[0]!=0.0&&order<=fOrder) s = fYx[order-1]/fW[0];
    else {
	s = 0.0;
	Error("GetYXNth","fOrder=%d fN=%d fW[0]=%f\n",
	       fOrder,fN,fW[0]);
    } // end else
    return s;
}
//______________________________________________________________________
Double_t AliITSstatistics2::GetRMSX()const{
    // This give the unbiased estimator for the RMS.
    // Inputs:
    //   none.
    // Outputs:
    //   none.
    // Return:
    //  The rms value
    Double_t x,x2,w,ww,s;

    x  = GetMeanX(); // first order
    x2 = GetXNth(2); // second order
    w  = fW[0];     // first order - 1.
    ww = fW[1];     // second order - 1.

    if(w*w==ww) return (-1.0);
    s = (x2-x*x)*w*w/(w*w-ww);
    return TMath::Sqrt(s);
}
//______________________________________________________________________
Double_t AliITSstatistics2::GetRMSY()const{
    // This give the unbiased estimator for the RMS.
    // Inputs:
    //   none.
    // Outputs:
    //   none.
    // Return:
    //  The rms value
    Double_t x,x2,w,ww,s;

    x  = GetMeanY(); // first order
    x2 = GetYNth(2); // second order
    w  = fW[0];     // first order - 1.
    ww = fW[1];     // second order - 1.

    if(w*w==ww) return (-1.0);
    s = (x2-x*x)*w*w/(w*w-ww);
    return TMath::Sqrt(s);
}
//______________________________________________________________________
Double_t AliITSstatistics2::GetRMSYX()const{
    // This give the unbiased estimator for the RMS.
    // Inputs:
    //   none.
    // Outputs:
    //   none.
    // Return:
    //  The rms value
    Double_t x,x2,w,ww,s;

    x  = GetMeanYX(); // first order
    x2 = GetYXNth(2); // second order
    w  = fW[0];     // first order - 1.
    ww = fW[1];     // second order - 1.

    if(w*w==ww) return (-1.0);
    s = (x2-x*x)*w*w/(w*w-ww);
    return TMath::Sqrt(s);
}
//______________________________________________________________________
Double_t AliITSstatistics2::GetErrorMeanY()const{
    //This is the error in the mean or the square root of the 
    // variance of the mean.
    // Inputs:
    //   none.
    // Outputs:
    //   none.
    // Return:
    //  The error on the mean
    Double_t rms,w,ww,s;

    rms = GetRMSY();
    w   = fW[0];
    ww  = fW[1];
    s   = rms*rms*ww/(w*w);
    return TMath::Sqrt(s);
}
//______________________________________________________________________
Double_t AliITSstatistics2::GetErrorMeanX()const{
    //This is the error in the mean or the square root of the 
    // variance of the mean.
    // Inputs:
    //   none.
    // Outputs:
    //   none.
    // Return:
    //  The error on the mean
    Double_t rms,w,ww,s;

    rms = GetRMSX();
    w   = fW[0];
    ww  = fW[1];
    s   = rms*rms*ww/(w*w);
    return TMath::Sqrt(s);
}
//______________________________________________________________________
Double_t AliITSstatistics2::GetErrorMeanYX()const{
    //This is the error in the mean or the square root of the 
    // variance of the mean.
    // Inputs:
    //   none.
    // Outputs:
    //   none.
    // Return:
    //  The error on the mean
    Double_t rms,w,ww,s;

    rms = GetRMSYX();
    w   = fW[0];
    ww  = fW[1];
    s   = rms*rms*ww/(w*w);
    return TMath::Sqrt(s);
}
//______________________________________________________________________
Double_t AliITSstatistics2::GetErrorRMSY()const{
    // This is the error in the mean or the square root of the variance 
    // of the mean. at this moment this routine is only defined for 
    // weights=1.
    // Inputs:
    //   none.
    // Outputs:
    //   none.
    // Return:
    //  The error on the rms
    Double_t x,x2,x3,x4,w,ww,m2,m4,n,s;

    if(fW[0]!=(Double_t)fN||GetN()<4) return (-1.);
    x  = GetMeanY(); // first order
    x2 = GetYNth(2); // second order
    w  = fW[0];     // first order - 1.
    ww = fW[1];     // second order - 1.
    if(w*w==ww) return (-1.0);
    s = (x2-x*x)*w*w/(w*w-ww);

    m2  = s;
    n   = (Double_t) GetN();
    x3  = GetYNth(3);
    x4  = GetYNth(4);
    // This equation assumes that all of the weights are equal to 1.
    m4  = (n/(n-1.))*(x4-3.*x*x3+6.*x*x*x2-2.*x*x*x*x);
    s   = (m4-(n-3.)*m2*m2/(n-1.))/n;
    return TMath::Sqrt(s);
}
//______________________________________________________________________
Double_t AliITSstatistics2::GetErrorRMSX()const{
    // This is the error in the mean or the square root of the variance 
    // of the mean. at this moment this routine is only defined for 
    // weights=1.
    // Inputs:
    //   none.
    // Outputs:
    //   none.
    // Return:
    //  The error on the rms
    Double_t x,x2,x3,x4,w,ww,m2,m4,n,s;

    if(fW[0]!=(Double_t)fN||GetN()<4) return (-1.);
    x  = GetMeanX(); // first order
    x2 = GetXNth(2); // second order
    w  = fW[0];     // first order - 1.
    ww = fW[1];     // second order - 1.
    if(w*w==ww) return (-1.0);
    s = (x2-x*x)*w*w/(w*w-ww);

    m2  = s;
    n   = (Double_t) GetN();
    x3  = GetXNth(3);
    x4  = GetXNth(4);
// This equation assumes that all of the weights are equal to 1.
    m4  = (n/(n-1.))*(x4-3.*x*x3+6.*x*x*x2-2.*x*x*x*x);
    s   = (m4-(n-3.)*m2*m2/(n-1.))/n;
    return TMath::Sqrt(s);
}
//______________________________________________________________________
Double_t AliITSstatistics2::GetErrorRMSYX()const{
    // This is the error in the mean or the square root of the variance 
    // of the mean. at this moment this routine is only defined for 
    // weights=1.
    // Inputs:
    //   none.
    // Outputs:
    //   none.
    // Return:
    //  The error on the rms
    Double_t x,x2,x3,x4,w,ww,m2,m4,n,s;

    if(fW[0]!=(Double_t)fN||GetN()<4) return (-1.);
    x  = GetMeanYX(); // first order
    x2 = GetYXNth(2); // second order
    w  = fW[0];     // first order - 1.
    ww = fW[1];     // second order - 1.
    if(w*w==ww) return (-1.0);
    s = (x2-x*x)*w*w/(w*w-ww);

    m2  = s;
    n   = (Double_t) GetN();
    x3  = GetYXNth(3);
    x4  = GetYXNth(4);
    // This equation assumes that all of the weights are equal to 1.
    m4  = (n/(n-1.))*(x4-3.*x*x3+6.*x*x*x2-2.*x*x*x*x);
    s   = (m4-(n-3.)*m2*m2/(n-1.))/n;
    return TMath::Sqrt(s);
}
//_______________________________________________________________________
Double_t AliITSstatistics2::FitToLine(Double_t &a,Double_t &b)const{
    // fit to y = ax+b returns Chi squared or -1.0 if an error
    // Inputs:
    //   none.
    // Outputs:
    //   Double_t  a   The slope parameter
    //   Double_t  b   The intercept paramter
    // Return:
    //  The Chi^2 of the fit
    Double_t c,d,e,f,g,h;

    a = b = 0.0;
    if(fOrder<2 || fN<3){
        Error("FitToLine","Order=%d<2 or N=%d<3",fOrder,fN);
        return -1.0;
    } // end if
    c = GetWN(1);
    d = GetYN(1);
    e = GetXN(1);
    f = GetYXN(1);
    g = GetXN(2);
    h = c*g-e*e;
    b = d*g-f*e;
    a = c*f-d*e;
    if(h!=0.0){
	a = a/h;
	b = b/h;
    }else{
	Error("FitToLine","vertical line: fOrder=%d fN=%d "
              "GetWN(1)=%g X GetXN(2)=%g - GetXN(1)=%g^2 = 0",fOrder,fN,c,g,e);
	return -1.0;
    } // end if h
    return GetChiSquared(a,b);
}
//_______________________________________________________________________
Double_t AliITSstatistics2::GetChiSquared(Double_t a,Double_t b)const{
    //  returns Chi^2 value of data to line y=ax+b with given a,b
    // Inputs:
    //   Double_t  a   The slope parameter
    //   Double_t  b   The intercept paramter
    // Outputs::
    //   none.
    // Return:
    //  The Chi^2 of the fit
    Double_t c2;

    c2 = GetYN(2)+b*b*GetWN(1)+
        a*a*GetXN(2)-2.0*b*GetYN(1)-2.0*a*GetYXN(1)+2.0*b*a*GetXN(1);
    c2 /= (Double_t)fN - 2.0;
    return c2;
}
//______________________________________________________________________
void AliITSstatistics2::PrintAscii(ostream *os)const{
    // Print out class data values in Ascii Form to output stream
    // Inputs:
    //   ostream *os   Output stream where Ascii data is to be writen
    // Outputs:
    //   none.
    // Return:
    //   none.
    Int_t i;
#if defined __GNUC__
#if __GNUC__ > 2
    ios::fmtflags fmt;
#else
    Int_t fmt;
#endif
#else
#if defined __ICC || defined __ECC || defined __xlC__
    ios::fmtflags fmt;
#else
    Int_t fmt;
#endif
#endif

    *os << fN <<" "<< fOrder;
    fmt = os->setf(ios::scientific); // set scientific floating point output
    for(i=0;i<fOrder;i++) *os <<" "<< fX[i];
    for(i=0;i<fOrder;i++) *os <<" "<< fYx[i];
    for(i=0;i<fOrder;i++) *os <<" "<< fY[i];
    for(i=0;i<fOrder;i++) *os <<" "<< fW[i];
    os->flags(fmt); // reset back to old Formating.
    return;
}
//______________________________________________________________________
void AliITSstatistics2::ReadAscii(istream *is){
    // Read in class data values in Ascii Form to output stream
    // Inputs:
    //   istream *is   Input stream where Ascii data is to be read in from
    // Outputs:
    //   none.
    // Return:
    //   none.
    Int_t i;

    *is >> i >> fOrder;
    Reset(fOrder);
    fN = i;
    for(i=0;i<fOrder;i++) *is >> fX[i];
    for(i=0;i<fOrder;i++) *is >> fYx[i];
    for(i=0;i<fOrder;i++) *is >> fY[i];
    for(i=0;i<fOrder;i++) *is >> fW[i];
}
//______________________________________________________________________
ostream &operator<<(ostream &os,const AliITSstatistics2 &s){
    // Standard output streaming function
    // Inputs:
    //   ostream            &os  output steam
    //   AliITSstatistics2 &s class to be streamed.
    // Output:
    //   none.
    // Return:
    //   ostream &os  The stream pointer

    s.PrintAscii(&os);
    return os;
}
//______________________________________________________________________
istream &operator>>(istream &is,AliITSstatistics2 &s){
    // Standard inputput streaming function
    // Inputs:
    //   istream            &is  input steam
    //   AliITSstatistics2 &s class to be streamed.
    // Output:
    //   none.
    // Return:
    //   ostream &os  The stream pointer

    s.ReadAscii(&is);
    return is;
}
Commit	Line	Data
b0f5e3fc	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"
f075a801	13	#include "Riostream.h"
b0f5e3fc	14
	15	ClassImp(AliITSstatistics2)
	16
	17	//
492b8715	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
b0f5e3fc	34	return;
b0f5e3fc	35	}
492b8715	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;
b0f5e3fc	49
492b8715	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
b0f5e3fc	56	fN = 0;
	57	return;
	58	}
492b8715	59	//______________________________________________________________________
b0f5e3fc	60	AliITSstatistics2::~AliITSstatistics2(){
492b8715	61	// destructor
	62	// Inputs:
	63	// none.
	64	// Outputs:
	65	// none.
	66	// Return:
	67	// none.
	68
07025d6d	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;
b0f5e3fc	77	fN = 0;
	78	fOrder = 0;
	79	}
b0f5e3fc	80	//_______________________________________________________________
b0f5e3fc	81	AliITSstatistics2& AliITSstatistics2::operator=(AliITSstatistics2 &source){
492b8715	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;
b0f5e3fc	102	}
b0f5e3fc	103	//_______________________________________________________________
492b8715	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
b0f5e3fc	127	}
492b8715	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;
b0f5e3fc	138
492b8715	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
b0f5e3fc	145	fN = 0;
492b8715	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];
b0f5e3fc	163	return;
b0f5e3fc	164	}
492b8715	165	//______________________________________________________________________
b0f5e3fc	166	void AliITSstatistics2::AddValue(Double_t y,Double_t x,Double_t w=1.0){
492b8715	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.
b0f5e3fc	176	Double_t xs=1.0,ys=1.0,yxs=1.0,ws=1.0;
b0f5e3fc	177	Int_t i;
e8189707	178	const Double_t kBig=1.0e+38;
	179
	180	if(y>kBig \|\| x>kBig \|\| w>kBig) return;
b0f5e3fc	181	fN++;
	182	for(i=0;i<fOrder;i++){
	183	xs *= x;
	184	ys *= y;
	185	yxs = xy;
	186	ws *= w;
07025d6d	187	fX[i] += xs*w;
	188	fY[i] += ys*w;
	189	fYx[i] += yxs*w;
	190	fW[i] += ws;
b0f5e3fc	191	} // end for i
b0f5e3fc	192	}
492b8715	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}.
b0f5e3fc	202	Double_t s;
b0f5e3fc	203
07025d6d	204	if(fW[0]!=0.0&&order<=fOrder) s = fX[order-1]/fW[0];
b0f5e3fc	205	else {
b0f5e3fc	206	s = 0.0;
492b8715	207	Error("GetXNth","error fOrder=%d fN=%d fW[0]=%f\n",
492b8715	208	fOrder,fN,fW[0]);
b0f5e3fc	209	} // end else
	210	return s;
	211	}
492b8715	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}
b0f5e3fc	221	Double_t s;
b0f5e3fc	222
07025d6d	223	if(fW[0]!=0.0&&order<=fOrder) s = fY[order-1]/fW[0];
b0f5e3fc	224	else {
b0f5e3fc	225	s = 0.0;
492b8715	226	Error("GetYNth","fOrder=%d fN=%d fW[0]=%f\n",
492b8715	227	fOrder,fN,fW[0]);
b0f5e3fc	228	} // end else
	229	return s;
	230	}
492b8715	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}
b0f5e3fc	240	Double_t s;
b0f5e3fc	241
07025d6d	242	if(fW[0]!=0.0&&order<=fOrder) s = fYx[order-1]/fW[0];
b0f5e3fc	243	else {
b0f5e3fc	244	s = 0.0;
492b8715	245	Error("GetYXNth","fOrder=%d fN=%d fW[0]=%f\n",
07025d6d	246	fOrder,fN,fW[0]);
b0f5e3fc	247	} // end else
	248	return s;
	249	}
492b8715	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
b0f5e3fc	259	Double_t x,x2,w,ww,s;
	260
	261	x = GetMeanX(); // first order
	262	x2 = GetXNth(2); // second order
07025d6d	263	w = fW[0]; // first order - 1.
07025d6d	264	ww = fW[1]; // second order - 1.
b0f5e3fc	265
	266	if(w*w==ww) return (-1.0);
	267	s = (x2-xx)ww/(ww-ww);
	268	return TMath::Sqrt(s);
	269	}
492b8715	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
b0f5e3fc	279	Double_t x,x2,w,ww,s;
	280
	281	x = GetMeanY(); // first order
	282	x2 = GetYNth(2); // second order
07025d6d	283	w = fW[0]; // first order - 1.
07025d6d	284	ww = fW[1]; // second order - 1.
b0f5e3fc	285
	286	if(w*w==ww) return (-1.0);
	287	s = (x2-xx)ww/(ww-ww);
	288	return TMath::Sqrt(s);
	289	}
492b8715	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
b0f5e3fc	299	Double_t x,x2,w,ww,s;
	300
	301	x = GetMeanYX(); // first order
	302	x2 = GetYXNth(2); // second order
07025d6d	303	w = fW[0]; // first order - 1.
07025d6d	304	ww = fW[1]; // second order - 1.
b0f5e3fc	305
	306	if(w*w==ww) return (-1.0);
	307	s = (x2-xx)ww/(ww-ww);
	308	return TMath::Sqrt(s);
	309	}
492b8715	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
b0f5e3fc	320	Double_t rms,w,ww,s;
	321
	322	rms = GetRMSY();
07025d6d	323	w = fW[0];
07025d6d	324	ww = fW[1];
b0f5e3fc	325	s = rmsrmsww/(w*w);
	326	return TMath::Sqrt(s);
	327	}
492b8715	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
b0f5e3fc	338	Double_t rms,w,ww,s;
	339
	340	rms = GetRMSX();
07025d6d	341	w = fW[0];
07025d6d	342	ww = fW[1];
b0f5e3fc	343	s = rmsrmsww/(w*w);
	344	return TMath::Sqrt(s);
	345	}
492b8715	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
b0f5e3fc	356	Double_t rms,w,ww,s;
	357
	358	rms = GetRMSYX();
07025d6d	359	w = fW[0];
07025d6d	360	ww = fW[1];
b0f5e3fc	361	s = rmsrmsww/(w*w);
	362	return TMath::Sqrt(s);
	363	}
492b8715	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
b0f5e3fc	375	Double_t x,x2,x3,x4,w,ww,m2,m4,n,s;
b0f5e3fc	376
07025d6d	377	if(fW[0]!=(Double_t)fN\|\|GetN()<4) return (-1.);
b0f5e3fc	378	x = GetMeanY(); // first order
b0f5e3fc	379	x2 = GetYNth(2); // second order
07025d6d	380	w = fW[0]; // first order - 1.
07025d6d	381	ww = fW[1]; // second order - 1.
b0f5e3fc	382	if(w*w==ww) return (-1.0);
	383	s = (x2-xx)ww/(ww-ww);
	384
	385	m2 = s;
	386	n = (Double_t) GetN();
	387	x3 = GetYNth(3);
	388	x4 = GetYNth(4);
492b8715	389	// This equation assumes that all of the weights are equal to 1.
b0f5e3fc	390	m4 = (n/(n-1.))(x4-3.xx3+6.xxx2-2.xxxx);
	391	s = (m4-(n-3.)m2m2/(n-1.))/n;
	392	return TMath::Sqrt(s);
	393	}
492b8715	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
b0f5e3fc	405	Double_t x,x2,x3,x4,w,ww,m2,m4,n,s;
b0f5e3fc	406
07025d6d	407	if(fW[0]!=(Double_t)fN\|\|GetN()<4) return (-1.);
b0f5e3fc	408	x = GetMeanX(); // first order
b0f5e3fc	409	x2 = GetXNth(2); // second order
07025d6d	410	w = fW[0]; // first order - 1.
07025d6d	411	ww = fW[1]; // second order - 1.
b0f5e3fc	412	if(w*w==ww) return (-1.0);
	413	s = (x2-xx)ww/(ww-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.xx3+6.xxx2-2.xxxx);
	421	s = (m4-(n-3.)m2m2/(n-1.))/n;
	422	return TMath::Sqrt(s);
	423	}
492b8715	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
b0f5e3fc	435	Double_t x,x2,x3,x4,w,ww,m2,m4,n,s;
b0f5e3fc	436
07025d6d	437	if(fW[0]!=(Double_t)fN\|\|GetN()<4) return (-1.);
b0f5e3fc	438	x = GetMeanYX(); // first order
b0f5e3fc	439	x2 = GetYXNth(2); // second order
07025d6d	440	w = fW[0]; // first order - 1.
07025d6d	441	ww = fW[1]; // second order - 1.
b0f5e3fc	442	if(w*w==ww) return (-1.0);
	443	s = (x2-xx)ww/(ww-ww);
	444
	445	m2 = s;
	446	n = (Double_t) GetN();
	447	x3 = GetYXNth(3);
	448	x4 = GetYXNth(4);
492b8715	449	// This equation assumes that all of the weights are equal to 1.
b0f5e3fc	450	m4 = (n/(n-1.))(x4-3.xx3+6.xxx2-2.xxxx);
	451	s = (m4-(n-3.)m2m2/(n-1.))/n;
	452	return TMath::Sqrt(s);
	453	}
	454	//_______________________________________________________________________
492b8715	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
b0f5e3fc	464	Double_t c,d,e,f,g,h;
	465
	466	a = b = 0.0;
492b8715	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
b0f5e3fc	471	c = GetWN(1);
	472	d = GetYN(1);
	473	e = GetXN(1);
	474	f = GetYXN(1);
	475	g = GetXN(2);
	476	h = cg-ee;
492b8715	477	b = dg-fe;
492b8715	478	a = cf-de;
b0f5e3fc	479	if(h!=0.0){
	480	a = a/h;
	481	b = b/h;
	482	}else{
492b8715	483	Error("FitToLine","vertical line: fOrder=%d fN=%d "
492b8715	484	"GetWN(1)=%g X GetXN(2)=%g - GetXN(1)=%g^2 = 0",fOrder,fN,c,g,e);
b0f5e3fc	485	return -1.0;
b0f5e3fc	486	} // end if h
492b8715	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)+bbGetWN(1)+
	502	aaGetXN(2)-2.0bGetYN(1)-2.0aGetYXN(1)+2.0ba*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;
b0f5e3fc	585	}