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
  14 ClassImp(AliITSstatistics2)
  15
  16 //
  17 AliITSstatistics2::AliITSstatistics2() : TObject(){
  18 //
  19 // default constructor
  20 //
  21     fx  = 0;
  22     fy  = 0;
  23     fyx = 0;
  24     fw  = 0;
  25     fN  = 0;
  26     fOrder = 0;
  27     return;
  28 }
  29
  30
  31 AliITSstatistics2::AliITSstatistics2(Int_t order) : TObject(){
  32 //
  33 // constructor to maximum moment/order order
  34 //
  35     Int_t i;
  36     fOrder = order;
  37     fx     = new Double_t[order];
  38     fy     = new Double_t[order];
  39     fyx    = new Double_t[order];
  40     fw     = new Double_t[order];
  41     for(i=0;i<order;i++) {fx[i] = 0.0;fy[i] = 0.0;
  42                           fyx[i] = 0.0; fw[i] = 0.0;}
  43     fN = 0;
  44     return;
  45 }
  46
  47 AliITSstatistics2::~AliITSstatistics2(){
  48 //
  49 // destructor
  50 //
  51     if(fx!=0)  delete[] fx;
  52     if(fy!=0)  delete[] fy;
  53     if(fyx!=0) delete[] fyx;
  54     if(fw!=0)  delete[] fw;
  55     fx  = 0;
  56     fy  = 0;
  57     fyx = 0;
  58     fw  = 0;
  59     fN  = 0;
  60     fOrder = 0;
  61 }
  62
  63 //_______________________________________________________________
  64 AliITSstatistics2& AliITSstatistics2::operator=(AliITSstatistics2 &source){
  65 // operator =
  66
  67   Int_t i;
  68   if(this==&source) return *this;
  69   if(source.fOrder!=0){
  70     this->fOrder = source.fOrder;
  71     this->fN = source.fN;
  72     this->fx = new Double_t[this->fOrder];
  73     this->fw = new Double_t[this->fOrder];
  74     for(i=0;i<source.fOrder;i++){
  75       this->fx[i] = source.fx[i];
  76       this->fw[i] = source.fw[i];
  77     } // end for i
  78   }else{
  79     this->fx = 0;
  80     this->fw = 0;
  81     this->fN = 0;
  82     this->fOrder = 0;
  83   }// end if source.fOrder!=0
  84   return *this;
  85 }
  86 //_______________________________________________________________
  87 AliITSstatistics2::AliITSstatistics2(AliITSstatistics2 &source){
  88 // Copy constructor
  89
  90   Int_t i;
  91   if(this==&source) return;
  92   if(source.fOrder!=0){
  93     this->fOrder = source.fOrder;
  94     this->fN = source.fN;
  95     this->fx = new Double_t[this->fOrder];
  96     this->fw = new Double_t[this->fOrder];
  97     for(i=0;i<source.fOrder;i++){
  98       this->fx[i] = source.fx[i];
  99       this->fw[i] = source.fw[i];
 100     } // end for i
 101   }else{
 102     this->fx = 0;
 103     this->fw = 0;
 104     this->fN = 0;
 105     this->fOrder = 0;
 106   }// end if source.fOrder!=0
 107 }
 108
 109 void AliITSstatistics2::Reset(){
 110 //
 111 // Reset/zero statistics
 112 //
 113     Int_t i;
 114     for(i=0;i<fOrder;i++) {fx[i] = 0.0;fy[i] = 0.0;
 115                            fyx[i] = 0.0; fw[i] = 0.0;}
 116     fN = 0;
 117     return;
 118 }
 119
 120 void AliITSstatistics2::AddValue(Double_t y,Double_t x,Double_t w=1.0){
 121 //
 122 // add next x,y pair to statistics
 123 //
 124     Double_t xs=1.0,ys=1.0,yxs=1.0,ws=1.0;
 125     Int_t i;
 126
 127     if(isinf(y)!=0||isinf(x)!=0||isinf(w)!=0) return;
 128     if(isnan(y)!=0||isnan(x)!=0||isnan(w)!=0) return;
 129     fN++;
 130     for(i=0;i<fOrder;i++){
 131         xs  *= x;
 132         ys  *= y;
 133         yxs *= x*y;
 134         ws  *= w;
 135         fx[i]  += xs*w;
 136         fy[i]  += ys*w;
 137         fyx[i] += yxs*w;
 138         fw[i]  += ws;
 139     } // end for i
 140 }
 141
 142 Double_t AliITSstatistics2::GetXNth(Int_t order){
 143 //
 144 // This give the unbiased estimator for the RMS.
 145 //
 146
 147     Double_t s;
 148
 149     if(fw[0]!=0.0&&order<=fOrder) s = fx[order-1]/fw[0];
 150     else {
 151         s = 0.0;
 152         printf("AliITSstatistics2: error in GetNth: fOrder=%d fN=%d fw[0]=%f\n",
 153                fOrder,fN,fw[0]);
 154     } // end else
 155     return s;
 156 }
 157 Double_t AliITSstatistics2::GetYNth(Int_t order){
 158 //
 159 // This give the unbiased estimator for the RMS.
 160 //
 161     Double_t s;
 162
 163     if(fw[0]!=0.0&&order<=fOrder) s = fy[order-1]/fw[0];
 164     else {
 165         s = 0.0;
 166         printf("AliITSstatistics2: error in GetNth: fOrder=%d fN=%d fw[0]=%f\n",
 167                fOrder,fN,fw[0]);
 168     } // end else
 169     return s;
 170 }
 171 Double_t AliITSstatistics2::GetYXNth(Int_t order){
 172 // This give the unbiased estimator for the RMS.
 173     Double_t s;
 174
 175     if(fw[0]!=0.0&&order<=fOrder) s = fyx[order-1]/fw[0];
 176     else {
 177         s = 0.0;
 178         printf("AliITSstatistics2: error in GetNth: fOrder=%d fN=%d fw[0]=%f\n",
 179                fOrder,fN,fw[0]);
 180     } // end else
 181     return s;
 182 }
 183 Double_t AliITSstatistics2::GetRMSX(){
 184 // This give the unbiased estimator for the RMS.
 185     Double_t x,x2,w,ww,s;
 186
 187     x  = GetMeanX(); // first order
 188     x2 = GetXNth(2); // second order
 189     w  = fw[0];     // first order - 1.
 190     ww = fw[1];     // second order - 1.
 191
 192     if(w*w==ww) return (-1.0);
 193     s = (x2-x*x)*w*w/(w*w-ww);
 194     return TMath::Sqrt(s);
 195 }
 196 Double_t AliITSstatistics2::GetRMSY(){
 197 // This give the unbiased estimator for the RMS.
 198     Double_t x,x2,w,ww,s;
 199
 200     x  = GetMeanY(); // first order
 201     x2 = GetYNth(2); // second order
 202     w  = fw[0];     // first order - 1.
 203     ww = fw[1];     // second order - 1.
 204
 205     if(w*w==ww) return (-1.0);
 206     s = (x2-x*x)*w*w/(w*w-ww);
 207     return TMath::Sqrt(s);
 208 }
 209 Double_t AliITSstatistics2::GetRMSYX(){
 210 // This give the unbiased estimator for the RMS.
 211     Double_t x,x2,w,ww,s;
 212
 213     x  = GetMeanYX(); // first order
 214     x2 = GetYXNth(2); // second order
 215     w  = fw[0];     // first order - 1.
 216     ww = fw[1];     // second order - 1.
 217
 218     if(w*w==ww) return (-1.0);
 219     s = (x2-x*x)*w*w/(w*w-ww);
 220     return TMath::Sqrt(s);
 221 }
 222 Double_t AliITSstatistics2::GetErrorMeanY(){
 223 //This is the error in the mean or the square root of the variance of the mean.
 224     Double_t rms,w,ww,s;
 225
 226     rms = GetRMSY();
 227     w   = fw[0];
 228     ww  = fw[1];
 229     s   = rms*rms*ww/(w*w);
 230     return TMath::Sqrt(s);
 231 }
 232 Double_t AliITSstatistics2::GetErrorMeanX(){
 233 //This is the error in the mean or the square root of the variance of the mean.
 234     Double_t rms,w,ww,s;
 235
 236     rms = GetRMSX();
 237     w   = fw[0];
 238     ww  = fw[1];
 239     s   = rms*rms*ww/(w*w);
 240     return TMath::Sqrt(s);
 241 }
 242 Double_t AliITSstatistics2::GetErrorMeanYX(){
 243 //This is the error in the mean or the square root of the variance of the mean.
 244     Double_t rms,w,ww,s;
 245
 246     rms = GetRMSYX();
 247     w   = fw[0];
 248     ww  = fw[1];
 249     s   = rms*rms*ww/(w*w);
 250     return TMath::Sqrt(s);
 251 }
 252
 253
 254 Double_t AliITSstatistics2::GetErrorRMSY(){
 255 //This is the error in the mean or the square root of the variance of the mean.
 256 // at this moment this routine is only defined for weights=1.
 257     Double_t x,x2,x3,x4,w,ww,m2,m4,n,s;
 258
 259     if(fw[0]!=(Double_t)fN||GetN()<4) return (-1.);
 260     x  = GetMeanY(); // first order
 261     x2 = GetYNth(2); // second order
 262     w  = fw[0];     // first order - 1.
 263     ww = fw[1];     // second order - 1.
 264     if(w*w==ww) return (-1.0);
 265     s = (x2-x*x)*w*w/(w*w-ww);
 266
 267     m2  = s;
 268     n   = (Double_t) GetN();
 269     x3  = GetYNth(3);
 270     x4  = GetYNth(4);
 271 // This equation assumes that all of the weights are equal to 1.
 272     m4  = (n/(n-1.))*(x4-3.*x*x3+6.*x*x*x2-2.*x*x*x*x);
 273     s   = (m4-(n-3.)*m2*m2/(n-1.))/n;
 274     return TMath::Sqrt(s);
 275 }
 276 Double_t AliITSstatistics2::GetErrorRMSX(){
 277 //This is the error in the mean or the square root of the variance of the mean.
 278 // at this moment this routine is only defined for weights=1.
 279     Double_t x,x2,x3,x4,w,ww,m2,m4,n,s;
 280
 281     if(fw[0]!=(Double_t)fN||GetN()<4) return (-1.);
 282     x  = GetMeanX(); // first order
 283     x2 = GetXNth(2); // second order
 284     w  = fw[0];     // first order - 1.
 285     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
 289     m2  = s;
 290     n   = (Double_t) GetN();
 291     x3  = GetXNth(3);
 292     x4  = GetXNth(4);
 293 // This equation assumes that all of the weights are equal to 1.
 294     m4  = (n/(n-1.))*(x4-3.*x*x3+6.*x*x*x2-2.*x*x*x*x);
 295     s   = (m4-(n-3.)*m2*m2/(n-1.))/n;
 296     return TMath::Sqrt(s);
 297 }
 298 Double_t AliITSstatistics2::GetErrorRMSYX(){
 299 //This is the error in the mean or the square root of the variance of the mean.
 300 // at this moment this routine is only defined for weights=1.
 301     Double_t x,x2,x3,x4,w,ww,m2,m4,n,s;
 302
 303     if(fw[0]!=(Double_t)fN||GetN()<4) return (-1.);
 304     x  = GetMeanYX(); // first order
 305     x2 = GetYXNth(2); // second order
 306     w  = fw[0];     // first order - 1.
 307     ww = fw[1];     // second order - 1.
 308     if(w*w==ww) return (-1.0);
 309     s = (x2-x*x)*w*w/(w*w-ww);
 310
 311     m2  = s;
 312     n   = (Double_t) GetN();
 313     x3  = GetYXNth(3);
 314     x4  = GetYXNth(4);
 315 // This equation assumes that all of the weights are equal to 1.
 316     m4  = (n/(n-1.))*(x4-3.*x*x3+6.*x*x*x2-2.*x*x*x*x);
 317     s   = (m4-(n-3.)*m2*m2/(n-1.))/n;
 318     return TMath::Sqrt(s);
 319 }
 320 //_______________________________________________________________________
 321 Double_t AliITSstatistics2::FitToLine(Double_t &a,Double_t &b){
 322 // fit to y = a+bx returns Chi squared or -1.0 if an error
 323     Double_t c,d,e,f,g,h;
 324
 325     a = b = 0.0;
 326     if(fOrder<2 || fN<3) return -1.0;
 327     c = GetWN(1);
 328     d = GetYN(1);
 329     e = GetXN(1);
 330     f = GetYXN(1);
 331     g = GetXN(2);
 332     h = c*g-e*e;
 333     a = d*g-f*e;
 334     b = c*f-d*e;
 335     if(h!=0.0){
 336         a = a/h;
 337         b = b/h;
 338     }else{
 339         printf("AliITSstatistics2: Error in FitToLine vertical line\n");
 340         return -1.0;
 341     } // end if h
 342     h = GetYN(2)+a*a*c+b*b*g-2.0*a*d-2.0*b*f+2.0*a*b*e;
 343     h /= (Double_t)fN - 2.0;
 344     return h;
 345 }
 346 //_______________________________________________________________________
 347 void AliITSstatistics2::Streamer(TBuffer &R__b){
 348    // Stream an object of class AliITSstatistics2.
 349
 350    if (R__b.IsReading()) {
 351       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
 352       TObject::Streamer(R__b);
 353       R__b >> fN;
 354       R__b >> fOrder;
 355       R__b.ReadArray(fy);
 356       R__b.ReadArray(fx);
 357       R__b.ReadArray(fyx);
 358       R__b.ReadArray(fw);
 359    } else {
 360       R__b.WriteVersion(AliITSstatistics2::IsA());
 361       TObject::Streamer(R__b);
 362       R__b << fN;
 363       R__b << fOrder;
 364       R__b.WriteArray(fy,fOrder);
 365       R__b.WriteArray(fx,fOrder);
 366       R__b.WriteArray(fyx,fOrder);
 367       R__b.WriteArray(fw,fOrder);
 368    }
 369 }