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