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     fOrder = order;
  36     fX     = new Double_t[order];
  37     fY     = new Double_t[order];
  38     fYx    = new Double_t[order];
  39     fW     = new Double_t[order];
  40     for(Int_t i=0;i<order;i++) {fX[i] = 0.0;fY[i] = 0.0;
  41                                 fYx[i] = 0.0; fW[i] = 0.0;}
  42     fN = 0;
  43     return;
  44 }
  45
  46 AliITSstatistics2::~AliITSstatistics2(){
  47 //
  48 // destructor
  49 //
  50     if(fX!=0)  delete[] fX;
  51     if(fY!=0)  delete[] fY;
  52     if(fYx!=0) delete[] fYx;
  53     if(fW!=0)  delete[] fW;
  54     fX  = 0;
  55     fY  = 0;
  56     fYx = 0;
  57     fW  = 0;
  58     fN  = 0;
  59     fOrder = 0;
  60 }
  61
  62 //_______________________________________________________________
  63 AliITSstatistics2& AliITSstatistics2::operator=(AliITSstatistics2 &source){
  64 // operator =
  65
  66      if(this==&source) return *this;
  67           if(source.fOrder!=0){
  68                this->fOrder = source.fOrder;
  69                          this->fN = source.fN;
  70                          this->fX = new Double_t[this->fOrder];
  71                          this->fW = new Double_t[this->fOrder];
  72                          for(Int_t i=0;i<source.fOrder;i++){
  73                               this->fX[i] = source.fX[i];
  74                                         this->fW[i] = source.fW[i];
  75                          } // end for i
  76           }else{
  77                this->fX = 0;
  78                          this->fW = 0;
  79                          this->fN = 0;
  80                          this->fOrder = 0;
  81           }// end if source.fOrder!=0
  82           return *this;
  83 }
  84 //_______________________________________________________________
  85 AliITSstatistics2::AliITSstatistics2(AliITSstatistics2 &source):
  86     TObject(source){
  87 // Copy constructor
  88
  89      if(this==&source) return;
  90           if(source.fOrder!=0){
  91                this->fOrder = source.fOrder;
  92                          this->fN = source.fN;
  93                          this->fX = new Double_t[this->fOrder];
  94                          this->fW = new Double_t[this->fOrder];
  95                          for(Int_t i=0;i<source.fOrder;i++){
  96                               this->fX[i] = source.fX[i];
  97                                         this->fW[i] = source.fW[i];
  98                          } // end for i
  99           }else{
 100                this->fX = 0;
 101                          this->fW = 0;
 102                          this->fN = 0;
 103                          this->fOrder = 0;
 104           }// end if source.fOrder!=0
 105 }
 106
 107 void AliITSstatistics2::Reset(){
 108 //
 109 // Reset/zero statistics
 110 //
 111     for(Int_t i=0;i<fOrder;i++) {fX[i] = 0.0;fY[i] = 0.0;
 112                                 fYx[i] = 0.0; fW[i] = 0.0;}
 113     fN = 0;
 114     return;
 115 }
 116
 117 void AliITSstatistics2::AddValue(Double_t y,Double_t x,Double_t w=1.0){
 118 //
 119 // add next x,y pair to statistics
 120 //
 121     Double_t xs=1.0,ys=1.0,yxs=1.0,ws=1.0;
 122     Int_t i;
 123
 124     const Double_t kBig=1.0e+38;
 125
 126     if(y>kBig || x>kBig || w>kBig) return;
 127
 128
 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 }