1 /**************************************************************************
2 * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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 //
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 //
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}
34 ClassImp(AliITSstatistics2)
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
48 // default constructor
54 // A default constructed AliITSstatistics class
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
71 // Int_t order The maximum moment of distributions {for example x^n}
74 for(i=0;i<order;i++) {
83 //______________________________________________________________________
84 AliITSstatistics2::~AliITSstatistics2(){
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;
106 //_______________________________________________________________
107 AliITSstatistics2& AliITSstatistics2::operator=(AliITSstatistics2 &source){
110 // AliITSstaticstics2 &source The source of this copy.
114 // A copy of the source class
116 if(this==&source) return *this;
117 TObject::operator=(source);
118 Reset(source.GetOrder());
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];
127 // this->fDig = source.fDig;
128 // if(fOver!=0) this->fOver = new AliITSstatistics2(*(source.fOver));
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
146 // AliITSstatistics2 & source the source of this copy
150 // A copy of the source.
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];
158 //if(fOver!=0) this->fOver = new AliITSstatistics2(*(source.fOver));
161 //______________________________________________________________________
162 void AliITSstatistics2::Reset(Int_t order){
163 // Reset/zero all statistics variables statistics
172 for(i=0;i<fOrder;i++) {
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;
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;
199 //----------------------------------------------------------------------
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.
208 // Int_t d The new significance level
217 //______________________________________________________________________
218 void AliITSstatistics2::AddValue(Double_t y,Double_t x,Double_t w=1.0){
219 // add next x,y pair to statistics
221 // Double_t y y value of pair
222 // Double_t x x value of pair
223 // Double_t w weight of pair
228 Double_t xs=1.0,ys=1.0,yxs=1.0,ws=1.0;
230 const Double_t kBig=1.0e+38;
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.
237 fOver = new AliITSstatistics2(fOrder);
239 fOver->AddValue(y,x,w);
244 for(i=0;i<GetOrder();i++){
255 //______________________________________________________________________
256 Double_t AliITSstatistics2::GetXNth(Int_t order)const{
257 // This give the unbiased estimator for the RMS.
259 // Int_t order the order of x^n value to be returned
263 // The value sum{x^n}.
266 if(GetWN(1)!=0.0 && order<=GetOrder()) s = GetXN(order)/GetWN(1);
269 Error("GetXNth","error fOrder=%d fN=%d fW[0]=%f\n",
270 GetOrder(),GetN(),GetWN(1));
274 //______________________________________________________________________
275 Double_t AliITSstatistics2::GetYNth(Int_t order)const{
276 // This give the unbiased estimator for the RMS.
278 // Int_t order the order of y^n value to be returned
282 // The value sum{y^n}
285 if(GetWN(1)!=0.0&&order<=GetOrder()) s = GetYN(order)/GetWN(1);
288 Error("GetYNth","fOrder=%d fN=%d fW[0]=%f\n",
289 GetOrder(),GetN(),GetWN(1));
293 //______________________________________________________________________
294 Double_t AliITSstatistics2::GetYXNth(Int_t order)const{
295 // This give the unbiased estimator for the RMS.
297 // Int_t order the order of (xy)^n value to be returned
301 // The value sum{(xy)^n}
304 if(GetWN(1)!=0.0&&order<=GetOrder()) s = GetYXN(order)/GetWN(1);
307 Error("GetYXNth","fOrder=%d fN=%d fW[0]=%f\n",
308 GetOrder(),GetN(),GetWN(1));
312 //______________________________________________________________________
313 Double_t AliITSstatistics2::GetRMSX()const{
314 // This give the unbiased estimator for the RMS.
321 Double_t x,x2,w,ww,s;
323 x = GetMeanX(); // first order
324 x2 = GetXNth(2); // second order
325 w = GetWN(1); // first order
326 ww = GetWN(2); // second order
328 if(w*w==ww) return (-1.0);
329 s = (x2-x*x)*w*w/(w*w-ww);
330 return TMath::Sqrt(s);
332 //______________________________________________________________________
333 Double_t AliITSstatistics2::GetRMSY()const{
334 // This give the unbiased estimator for the RMS.
341 Double_t x,x2,w,ww,s;
343 x = GetMeanY(); // first order
344 x2 = GetYNth(2); // second order
345 w = GetWN(1); // first order
346 ww = GetWN(2); // second order
348 if(w*w==ww) return (-1.0);
349 s = (x2-x*x)*w*w/(w*w-ww);
350 return TMath::Sqrt(s);
352 //______________________________________________________________________
353 Double_t AliITSstatistics2::GetRMSYX()const{
354 // This give the unbiased estimator for the RMS.
361 Double_t x,x2,w,ww,s;
363 x = GetMeanYX(); // first order
364 x2 = GetYXNth(2); // second order
365 w = GetWN(1); // first order
366 ww = GetWN(2); // second order
368 if(w*w==ww) return (-1.0);
369 s = (x2-x*x)*w*w/(w*w-ww);
370 return TMath::Sqrt(s);
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.
381 // The error on the mean
387 s = rms*rms*ww/(w*w);
388 return TMath::Sqrt(s);
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.
399 // The error on the mean
405 s = rms*rms*ww/(w*w);
406 return TMath::Sqrt(s);
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.
417 // The error on the mean
423 s = rms*rms*ww/(w*w);
424 return TMath::Sqrt(s);
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
436 // The error on the rms
437 Double_t x,x2,x3,x4,w,ww,m2,m4,n,s;
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);
448 n = (Double_t) GetN();
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);
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
466 // The error on the rms
467 Double_t x,x2,x3,x4,w,ww,m2,m4,n,s;
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);
478 n = (Double_t) GetN();
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);
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
496 // The error on the rms
497 Double_t x,x2,x3,x4,w,ww,m2,m4,n,s;
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);
508 n = (Double_t) GetN();
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);
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
524 \Chi^{2}=\sum_{i} (y_{i}-a x_{i} -b)^{2} w_{i}
526 Where if the weight used in
527 AliITSstatistics2::AddValue(Double_t y,Double_t x,Double_t w=1.0)
530 w_{i}=\frac{1}{\delta y^{2}}.
532 Then we get the typicall chi square minimization.
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
543 // The Chi^2 of the fit
544 Double_t c,d,e,f,g,h;
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());
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);
567 // Now for the errors.
568 ea = c*c*g+(a*a-1.0)*c*e*e;
571 Error("FitToLine","ea=%g is less than zero",ea);
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;
579 Error("FitToLine","eb=%g is less than zero",eb);
582 eb = TMath::Sqrt(eb);
583 c = GetChiSquared(a,b);
584 if(c<=0.0){ // must be a numerical precision problem.
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.
593 Note: The Chi^2 value is computed from the expression
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}
600 and not form the expression
602 \chi^{2}= \sum_{i}{(y_{i}-ax_{i}-b)^{2}w_{i}.
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
614 // Double_t a The slope parameter
615 // Double_t b The intercept paramter
619 // The Chi^2 of the fit
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;
627 //______________________________________________________________________
628 void AliITSstatistics2::PrintAscii(ostream *os)const{
629 // Print out class data values in Ascii Form to output stream
631 // ostream *os Output stream where Ascii data is to be writen
639 ios::fmtflags fmt; // Standard IO format object, required for output.
644 #if defined __ICC || defined __ECC || defined __xlC__
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.
663 //______________________________________________________________________
664 void AliITSstatistics2::ReadAscii(istream *is){
665 // Read in class data values in Ascii Form to output stream
667 // istream *is Input stream where Ascii data is to be read in from
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];
682 // if(fDig>0) *is >> fOver;
685 //______________________________________________________________________
686 ostream &operator<<(ostream &os,const AliITSstatistics2 &s){
687 // Standard output streaming function
689 // ostream &os output steam
690 // AliITSstatistics2 &s class to be streamed.
694 // ostream &os The stream pointer
699 //______________________________________________________________________
700 istream &operator>>(istream &is,AliITSstatistics2 &s){
701 // Standard inputput streaming function
703 // istream &is input steam
704 // AliITSstatistics2 &s class to be streamed.
708 // ostream &os The stream pointer