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 // Copy constructor
  87
  88      if(this==&source) return;
  89           if(source.fOrder!=0){
  90                this->fOrder = source.fOrder;
  91                          this->fN = source.fN;
  92                          this->fX = new Double_t[this->fOrder];
  93                          this->fW = new Double_t[this->fOrder];
  94                          for(Int_t i=0;i<source.fOrder;i++){
  95                               this->fX[i] = source.fX[i];
  96                                         this->fW[i] = source.fW[i];
  97                          } // end for i
  98           }else{
  99                this->fX = 0;
 100                          this->fW = 0;
 101                          this->fN = 0;
 102                          this->fOrder = 0;
 103           }// end if source.fOrder!=0
 104 }
 105
 106 void AliITSstatistics2::Reset(){
 107 //
 108 // Reset/zero statistics
 109 //
 110     for(Int_t i=0;i<fOrder;i++) {fX[i] = 0.0;fY[i] = 0.0;
 111                                 fYx[i] = 0.0; fW[i] = 0.0;}
 112     fN = 0;
 113     return;
 114 }
 115
 116 void AliITSstatistics2::AddValue(Double_t y,Double_t x,Double_t w=1.0){
 117 //
 118 // add next x,y pair to statistics
 119 //
 120     Double_t xs=1.0,ys=1.0,yxs=1.0,ws=1.0;
 121     Int_t i;
 122
 123     const Double_t kBig=1.0e+38;
 124
 125     if(y>kBig || x>kBig || w>kBig) return;
 126
 127
 128     fN++;
 129     for(i=0;i<fOrder;i++){
 130         xs  *= x;
 131         ys  *= y;
 132         yxs *= x*y;
 133         ws  *= w;
 134         fX[i]  += xs*w;
 135         fY[i]  += ys*w;
 136         fYx[i] += yxs*w;
 137         fW[i]  += ws;
 138     } // end for i
 139 }
 140
 141 Double_t AliITSstatistics2::GetXNth(Int_t order){
 142 //
 143 // This give the unbiased estimator for the RMS.
 144 //
 145
 146     Double_t s;
 147
 148     if(fW[0]!=0.0&&order<=fOrder) s = fX[order-1]/fW[0];
 149     else {
 150         s = 0.0;
 151         printf("AliITSstatistics2: error in GetNth: fOrder=%d fN=%d fW[0]=%f\n",
 152                fOrder,fN,fW[0]);
 153     } // end else
 154     return s;
 155 }
 156 Double_t AliITSstatistics2::GetYNth(Int_t order){
 157 //
 158 // This give the unbiased estimator for the RMS.
 159 //
 160     Double_t s;
 161
 162     if(fW[0]!=0.0&&order<=fOrder) s = fY[order-1]/fW[0];
 163     else {
 164         s = 0.0;
 165         printf("AliITSstatistics2: error in GetNth: fOrder=%d fN=%d fW[0]=%f\n",
 166                fOrder,fN,fW[0]);
 167     } // end else
 168     return s;
 169 }
 170 Double_t AliITSstatistics2::GetYXNth(Int_t order){
 171 // This give the unbiased estimator for the RMS.
 172     Double_t s;
 173
 174     if(fW[0]!=0.0&&order<=fOrder) s = fYx[order-1]/fW[0];
 175     else {
 176         s = 0.0;
 177         printf("AliITSstatistics2: error in GetNth: fOrder=%d fN=%d fW[0]=%f\n",
 178                fOrder,fN,fW[0]);
 179     } // end else
 180     return s;
 181 }
 182 Double_t AliITSstatistics2::GetRMSX(){
 183 // This give the unbiased estimator for the RMS.
 184     Double_t x,x2,w,ww,s;
 185
 186     x  = GetMeanX(); // first order
 187     x2 = GetXNth(2); // second order
 188     w  = fW[0];     // first order - 1.
 189     ww = fW[1];     // second order - 1.
 190
 191     if(w*w==ww) return (-1.0);
 192     s = (x2-x*x)*w*w/(w*w-ww);
 193     return TMath::Sqrt(s);
 194 }
 195 Double_t AliITSstatistics2::GetRMSY(){
 196 // This give the unbiased estimator for the RMS.
 197     Double_t x,x2,w,ww,s;
 198
 199     x  = GetMeanY(); // first order
 200     x2 = GetYNth(2); // second order
 201     w  = fW[0];     // first order - 1.
 202     ww = fW[1];     // second order - 1.
 203
 204     if(w*w==ww) return (-1.0);
 205     s = (x2-x*x)*w*w/(w*w-ww);
 206     return TMath::Sqrt(s);
 207 }
 208 Double_t AliITSstatistics2::GetRMSYX(){
 209 // This give the unbiased estimator for the RMS.
 210     Double_t x,x2,w,ww,s;
 211
 212     x  = GetMeanYX(); // first order
 213     x2 = GetYXNth(2); // second order
 214     w  = fW[0];     // first order - 1.
 215     ww = fW[1];     // second order - 1.
 216
 217     if(w*w==ww) return (-1.0);
 218     s = (x2-x*x)*w*w/(w*w-ww);
 219     return TMath::Sqrt(s);
 220 }
 221 Double_t AliITSstatistics2::GetErrorMeanY(){
 222 //This is the error in the mean or the square root of the variance of the mean.
 223     Double_t rms,w,ww,s;
 224
 225     rms = GetRMSY();
 226     w   = fW[0];
 227     ww  = fW[1];
 228     s   = rms*rms*ww/(w*w);
 229     return TMath::Sqrt(s);
 230 }
 231 Double_t AliITSstatistics2::GetErrorMeanX(){
 232 //This is the error in the mean or the square root of the variance of the mean.
 233     Double_t rms,w,ww,s;
 234
 235     rms = GetRMSX();
 236     w   = fW[0];
 237     ww  = fW[1];
 238     s   = rms*rms*ww/(w*w);
 239     return TMath::Sqrt(s);
 240 }
 241 Double_t AliITSstatistics2::GetErrorMeanYX(){
 242 //This is the error in the mean or the square root of the variance of the mean.
 243     Double_t rms,w,ww,s;
 244
 245     rms = GetRMSYX();
 246     w   = fW[0];
 247     ww  = fW[1];
 248     s   = rms*rms*ww/(w*w);
 249     return TMath::Sqrt(s);
 250 }
 251
 252
 253 Double_t AliITSstatistics2::GetErrorRMSY(){
 254 //This is the error in the mean or the square root of the variance of the mean.
 255 // at this moment this routine is only defined for weights=1.
 256     Double_t x,x2,x3,x4,w,ww,m2,m4,n,s;
 257
 258     if(fW[0]!=(Double_t)fN||GetN()<4) return (-1.);
 259     x  = GetMeanY(); // first order
 260     x2 = GetYNth(2); // second order
 261     w  = fW[0];     // first order - 1.
 262     ww = fW[1];     // second order - 1.
 263     if(w*w==ww) return (-1.0);
 264     s = (x2-x*x)*w*w/(w*w-ww);
 265
 266     m2  = s;
 267     n   = (Double_t) GetN();
 268     x3  = GetYNth(3);
 269     x4  = GetYNth(4);
 270 // This equation assumes that all of the weights are equal to 1.
 271     m4  = (n/(n-1.))*(x4-3.*x*x3+6.*x*x*x2-2.*x*x*x*x);
 272     s   = (m4-(n-3.)*m2*m2/(n-1.))/n;
 273     return TMath::Sqrt(s);
 274 }
 275 Double_t AliITSstatistics2::GetErrorRMSX(){
 276 //This is the error in the mean or the square root of the variance of the mean.
 277 // at this moment this routine is only defined for weights=1.
 278     Double_t x,x2,x3,x4,w,ww,m2,m4,n,s;
 279
 280     if(fW[0]!=(Double_t)fN||GetN()<4) return (-1.);
 281     x  = GetMeanX(); // first order
 282     x2 = GetXNth(2); // second order
 283     w  = fW[0];     // first order - 1.
 284     ww = fW[1];     // second order - 1.
 285     if(w*w==ww) return (-1.0);
 286     s = (x2-x*x)*w*w/(w*w-ww);
 287
 288     m2  = s;
 289     n   = (Double_t) GetN();
 290     x3  = GetXNth(3);
 291     x4  = GetXNth(4);
 292 // This equation assumes that all of the weights are equal to 1.
 293     m4  = (n/(n-1.))*(x4-3.*x*x3+6.*x*x*x2-2.*x*x*x*x);
 294     s   = (m4-(n-3.)*m2*m2/(n-1.))/n;
 295     return TMath::Sqrt(s);
 296 }
 297 Double_t AliITSstatistics2::GetErrorRMSYX(){
 298 //This is the error in the mean or the square root of the variance of the mean.
 299 // at this moment this routine is only defined for weights=1.
 300     Double_t x,x2,x3,x4,w,ww,m2,m4,n,s;
 301
 302     if(fW[0]!=(Double_t)fN||GetN()<4) return (-1.);
 303     x  = GetMeanYX(); // first order
 304     x2 = GetYXNth(2); // second order
 305     w  = fW[0];     // first order - 1.
 306     ww = fW[1];     // second order - 1.
 307     if(w*w==ww) return (-1.0);
 308     s = (x2-x*x)*w*w/(w*w-ww);
 309
 310     m2  = s;
 311     n   = (Double_t) GetN();
 312     x3  = GetYXNth(3);
 313     x4  = GetYXNth(4);
 314 // This equation assumes that all of the weights are equal to 1.
 315     m4  = (n/(n-1.))*(x4-3.*x*x3+6.*x*x*x2-2.*x*x*x*x);
 316     s   = (m4-(n-3.)*m2*m2/(n-1.))/n;
 317     return TMath::Sqrt(s);
 318 }
 319 //_______________________________________________________________________
 320 Double_t AliITSstatistics2::FitToLine(Double_t &a,Double_t &b){
 321 // fit to y = a+bx returns Chi squared or -1.0 if an error
 322     Double_t c,d,e,f,g,h;
 323
 324     a = b = 0.0;
 325     if(fOrder<2 || fN<3) return -1.0;
 326     c = GetWN(1);
 327     d = GetYN(1);
 328     e = GetXN(1);
 329     f = GetYXN(1);
 330     g = GetXN(2);
 331     h = c*g-e*e;
 332     a = d*g-f*e;
 333     b = c*f-d*e;
 334     if(h!=0.0){
 335         a = a/h;
 336         b = b/h;
 337     }else{
 338         printf("AliITSstatistics2: Error in FitToLine vertical line\n");
 339         return -1.0;
 340     } // end if h
 341     h = GetYN(2)+a*a*c+b*b*g-2.0*a*d-2.0*b*f+2.0*a*b*e;
 342     h /= (Double_t)fN - 2.0;
 343     return h;
 344 }