ITS/AliITSstatistics.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
  13 #include "AliITSstatistics.h"
  14
  15 ClassImp(AliITSstatistics)
  16
  17 //
  18 AliITSstatistics::AliITSstatistics() : TObject(){
  19 //
  20 // default contructor
  21 //
  22     fX = 0;
  23     fW = 0;
  24     fN = 0;
  25     fOrder = 0;
  26     return;
  27 }
  28
  29
  30 AliITSstatistics::AliITSstatistics(Int_t order) : TObject(){
  31 //
  32 // contructor to a specific order in the moments
  33 //
  34     fOrder = order;
  35     fX = new Double_t[order];
  36     fW = new Double_t[order];
  37     for(Int_t i=0;i<order;i++) {fX[i] = 0.0; fW[i] = 0.0;}
  38     fN = 0;
  39     return;
  40 }
  41
  42 AliITSstatistics::~AliITSstatistics(){
  43 //
  44 // default destructor
  45 //
  46     if(fX!=0) delete[] fX;
  47     if(fW!=0) delete[] fW;
  48     fX = 0;
  49     fW = 0;
  50     fN = 0;
  51     fOrder = 0;
  52 }
  53 //_______________________________________________________________
  54 AliITSstatistics& AliITSstatistics::operator=(AliITSstatistics &source){
  55 // operator =
  56
  57      if(this==&source) return *this;
  58           if(source.fOrder!=0){
  59                this->fOrder = source.fOrder;
  60                          this->fN = source.fN;
  61                          this->fX = new Double_t[this->fOrder];
  62                          this->fW = new Double_t[this->fOrder];
  63                          for(Int_t i=0;i<source.fOrder;i++){
  64                               this->fX[i] = source.fX[i];
  65                                         this->fW[i] = source.fW[i];
  66                          } // end for i
  67           }else{
  68                this->fX = 0;
  69                          this->fW = 0;
  70                          this->fN = 0;
  71                          this->fOrder = 0;
  72           }// end if source.fOrder!=0
  73           return *this;
  74 }
  75 //_______________________________________________________________
  76 AliITSstatistics::AliITSstatistics(AliITSstatistics &source){
  77 // Copy constructor
  78
  79      if(this==&source) return;
  80           if(source.fOrder!=0){
  81                this->fOrder = source.fOrder;
  82                          this->fN = source.fN;
  83                          this->fX = new Double_t[this->fOrder];
  84                          this->fW = new Double_t[this->fOrder];
  85                          for(Int_t i=0;i<source.fOrder;i++){
  86                               this->fX[i] = source.fX[i];
  87                                         this->fW[i] = source.fW[i];
  88                          } // end for i
  89           }else{
  90                this->fX = 0;
  91                          this->fW = 0;
  92                          this->fN = 0;
  93                          this->fOrder = 0;
  94           }// end if source.fOrder!=0
  95 }
  96 //_______________________________________________________________
  97 void AliITSstatistics::Reset(){
  98 //
  99 // reset all values to zero
 100 //
 101     for(Int_t i=0;i<fOrder;i++) {fX[i] = 0.0; fW[i] = 0.0;}
 102     fN = 0;
 103     return;
 104 }
 105
 106 //_______________________________________________________________
 107 void AliITSstatistics::AddValue(Double_t x,Double_t w){
 108 //
 109 // accumulate element x with weight w.
 110 //
 111
 112   //it was AddValue(Double_t x,Double_t w=1.0);
 113
 114     Double_t y=1.0,z=1.0;
 115
 116     Int_t i;
 117     const Double_t kBig=1.0e+38;
 118
 119     if(y>kBig || x>kBig || w>kBig) return;
 120
 121     fN++;
 122     for(i=0;i<fOrder;i++){
 123         y *= x;
 124         z *= w;
 125         fX[i] += y*w;
 126         fW[i] += z;
 127     } // end for i
 128 }
 129
 130 Double_t AliITSstatistics::GetNth(Int_t order){
 131 // This give the unbiased estimator for the RMS.
 132     Double_t s;
 133
 134     if(fW[0]!=0.0&&order<=fOrder) s = fX[order-1]/fW[0];
 135     else {
 136         s = 0.0;
 137         printf("AliITSstatistics: error in GetNth: fOrder=%d fN=%d fW[0]=%f\n",
 138                fOrder,fN,fW[0]);
 139     } // end else
 140     return s;
 141 }
 142
 143 Double_t AliITSstatistics::GetRMS(){
 144 // This give the unbiased estimator for the RMS.
 145     Double_t x,x2,w,ww,s;
 146
 147     x  = GetMean(); // first order
 148     x2 = GetNth(2); // second order
 149     w  = fW[0];     // first order - 1.
 150     ww = fW[1];     // second order - 1.
 151
 152     if(w*w==ww) return (-1.0);
 153     s = (x2-x*x)*w*w/(w*w-ww);
 154     return TMath::Sqrt(s);
 155 }
 156
 157 Double_t AliITSstatistics::GetErrorMean(){
 158 //This is the error in the mean or the square root of the variance of the mean.
 159     Double_t rms,w,ww,s;
 160
 161     rms = GetRMS();
 162     w   = fW[0];
 163     ww  = fW[1];
 164     s   = rms*rms*ww/(w*w);
 165     return TMath::Sqrt(s);
 166 }
 167
 168
 169 Double_t AliITSstatistics::GetErrorRMS(){
 170 //This is the error in the mean or the square root of the variance of the mean.
 171 // at this moment this routine is only defined for weights=1.
 172     Double_t x,x2,x3,x4,w,ww,m2,m4,n,s;
 173
 174     if(fW[0]!=(Double_t)fN||GetN()<4) return (-1.);
 175     x  = GetMean(); // first order
 176     x2 = GetNth(2); // second order
 177     w  = fW[0];     // first order - 1.
 178     ww = fW[1];     // second order - 1.
 179     if(w*w==ww) return (-1.0);
 180     s = (x2-x*x)*w*w/(w*w-ww);
 181
 182     m2  = s;
 183     n   = (Double_t) GetN();
 184     x3  = GetNth(3);
 185     x4  = GetNth(4);
 186 // This equation assumes that all of the weights are equal to 1.
 187     m4  = (n/(n-1.))*(x4-3.*x*x3+6.*x*x*x2-2.*x*x*x*x);
 188     s   = (m4-(n-3.)*m2*m2/(n-1.))/n;
 189     return TMath::Sqrt(s);
 190 }