1 //////////////////////////////////////////////////////////////////////////
2 // Alice ITS class to help keep statistical information //
4 // version: 0.0.0 Draft. //
5 // Date: April 18 1999 //
6 // By: Bjorn S. Nilsen //
8 //////////////////////////////////////////////////////////////////////////
12 #include "AliITSstatistics2.h"
13 #include "Riostream.h"
15 ClassImp(AliITSstatistics2)
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
32 // A default constructed AliITSstatistics class
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
47 // Int_t order The maximum moment of distributions {for example x^n}
50 for(i=0;i<order;i++) {
59 //______________________________________________________________________
60 AliITSstatistics2::~AliITSstatistics2(){
69 if(fX!=0) delete[] fX;
70 if(fY!=0) delete[] fY;
71 if(fYx!=0) delete[] fYx;
72 if(fW!=0) delete[] fW;
80 //_______________________________________________________________
81 AliITSstatistics2& AliITSstatistics2::operator=(AliITSstatistics2 &source){
84 // AliITSstaticstics2 &source The source of this copy.
88 // A copy of the source class
90 if(this==&source) return *this;
91 TObject::operator=(source);
92 Reset(source.GetOrder());
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];
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)
114 // AliITSstatistics2 & source the source of this copy
118 // A copy of the source.
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];
128 //______________________________________________________________________
129 void AliITSstatistics2::Reset(Int_t order){
130 // Reset/zero all statistics variables statistics
139 for(i=0;i<fOrder;i++) {
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;
159 fX = new Double_t[fOrder];
160 fY = new Double_t[fOrder];
161 fYx = new Double_t[fOrder];
162 fW = new Double_t[fOrder];
165 //______________________________________________________________________
166 void AliITSstatistics2::AddValue(Double_t y,Double_t x,Double_t w=1.0){
167 // add next x,y pair to statistics
169 // Double_t y y value of pair
170 // Double_t x x value of pair
171 // Double_t w weight of pair
176 Double_t xs=1.0,ys=1.0,yxs=1.0,ws=1.0;
178 const Double_t kBig=1.0e+38;
180 if(y>kBig || x>kBig || w>kBig) return;
182 for(i=0;i<fOrder;i++){
193 //______________________________________________________________________
194 Double_t AliITSstatistics2::GetXNth(Int_t order)const{
195 // This give the unbiased estimator for the RMS.
197 // Int_t order the order of x^n value to be returned
201 // The value sum{x^n}.
204 if(fW[0]!=0.0&&order<=fOrder) s = fX[order-1]/fW[0];
207 Error("GetXNth","error fOrder=%d fN=%d fW[0]=%f\n",
212 //______________________________________________________________________
213 Double_t AliITSstatistics2::GetYNth(Int_t order)const{
214 // This give the unbiased estimator for the RMS.
216 // Int_t order the order of y^n value to be returned
220 // The value sum{y^n}
223 if(fW[0]!=0.0&&order<=fOrder) s = fY[order-1]/fW[0];
226 Error("GetYNth","fOrder=%d fN=%d fW[0]=%f\n",
231 //______________________________________________________________________
232 Double_t AliITSstatistics2::GetYXNth(Int_t order)const{
233 // This give the unbiased estimator for the RMS.
235 // Int_t order the order of (xy)^n value to be returned
239 // The value sum{(xy)^n}
242 if(fW[0]!=0.0&&order<=fOrder) s = fYx[order-1]/fW[0];
245 Error("GetYXNth","fOrder=%d fN=%d fW[0]=%f\n",
250 //______________________________________________________________________
251 Double_t AliITSstatistics2::GetRMSX()const{
252 // This give the unbiased estimator for the RMS.
259 Double_t x,x2,w,ww,s;
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.
266 if(w*w==ww) return (-1.0);
267 s = (x2-x*x)*w*w/(w*w-ww);
268 return TMath::Sqrt(s);
270 //______________________________________________________________________
271 Double_t AliITSstatistics2::GetRMSY()const{
272 // This give the unbiased estimator for the RMS.
279 Double_t x,x2,w,ww,s;
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.
286 if(w*w==ww) return (-1.0);
287 s = (x2-x*x)*w*w/(w*w-ww);
288 return TMath::Sqrt(s);
290 //______________________________________________________________________
291 Double_t AliITSstatistics2::GetRMSYX()const{
292 // This give the unbiased estimator for the RMS.
299 Double_t x,x2,w,ww,s;
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.
306 if(w*w==ww) return (-1.0);
307 s = (x2-x*x)*w*w/(w*w-ww);
308 return TMath::Sqrt(s);
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.
319 // The error on the mean
325 s = rms*rms*ww/(w*w);
326 return TMath::Sqrt(s);
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.
337 // The error on the mean
343 s = rms*rms*ww/(w*w);
344 return TMath::Sqrt(s);
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.
355 // The error on the mean
361 s = rms*rms*ww/(w*w);
362 return TMath::Sqrt(s);
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
374 // The error on the rms
375 Double_t x,x2,x3,x4,w,ww,m2,m4,n,s;
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);
386 n = (Double_t) GetN();
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);
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
404 // The error on the rms
405 Double_t x,x2,x3,x4,w,ww,m2,m4,n,s;
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);
416 n = (Double_t) GetN();
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);
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
434 // The error on the rms
435 Double_t x,x2,x3,x4,w,ww,m2,m4,n,s;
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);
446 n = (Double_t) GetN();
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);
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
460 // Double_t a The slope parameter
461 // Double_t b The intercept paramter
463 // The Chi^2 of the fit
464 Double_t c,d,e,f,g,h;
467 if(fOrder<2 || fN<3){
468 Error("FitToLine","Order=%d<2 or N=%d<3",fOrder,fN);
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);
487 return GetChiSquared(a,b);
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
493 // Double_t a The slope parameter
494 // Double_t b The intercept paramter
498 // The Chi^2 of the fit
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;
506 //______________________________________________________________________
507 void AliITSstatistics2::PrintAscii(ostream *os)const{
508 // Print out class data values in Ascii Form to output stream
510 // ostream *os Output stream where Ascii data is to be writen
523 #if defined __ICC || defined __ECC || defined __xlC__
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.
539 //______________________________________________________________________
540 void AliITSstatistics2::ReadAscii(istream *is){
541 // Read in class data values in Ascii Form to output stream
543 // istream *is Input stream where Ascii data is to be read in from
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];
558 //______________________________________________________________________
559 ostream &operator<<(ostream &os,const AliITSstatistics2 &s){
560 // Standard output streaming function
562 // ostream &os output steam
563 // AliITSstatistics2 &s class to be streamed.
567 // ostream &os The stream pointer
572 //______________________________________________________________________
573 istream &operator>>(istream &is,AliITSstatistics2 &s){
574 // Standard inputput streaming function
576 // istream &is input steam
577 // AliITSstatistics2 &s class to be streamed.
581 // ostream &os The stream pointer