STEER/AliMathBase.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 ///////////////////////////////////////////////////////////////////////////
  18 // Class AliMathBase
  19 //
  20 // Subset of  matheamtical functions  not included in the TMath
  21 //
  22
  23 ///////////////////////////////////////////////////////////////////////////
  24 #include "TMath.h"
  25 #include "AliMathBase.h"
  26 #include "Riostream.h"
  27 #include "TH1F.h"
  28 #include "TF1.h"
  29 #include "TLinearFitter.h"
  30
  31 ClassImp(AliMathBase) // Class implementation to enable ROOT I/O
  32
  33 AliMathBase::AliMathBase() : TObject()
  34 {
  35 // Default constructor
  36 }
  37 ///////////////////////////////////////////////////////////////////////////
  38 AliMathBase::~AliMathBase()
  39 {
  40 // Destructor
  41 }
  42
  43
  44 //_____________________________________________________________________________
  45 void AliMathBase::EvaluateUni(Int_t nvectors, Double_t *data, Double_t &mean
  46                            , Double_t &sigma, Int_t hh)
  47 {
  48   //
  49   // Robust estimator in 1D case MI version - (faster than ROOT version)
  50   //
  51   // For the univariate case
  52   // estimates of location and scatter are returned in mean and sigma parameters
  53   // the algorithm works on the same principle as in multivariate case -
  54   // it finds a subset of size hh with smallest sigma, and then returns mean and
  55   // sigma of this subset
  56   //
  57
  58   if (hh==0)
  59     hh=(nvectors+2)/2;
  60   Double_t faclts[]={2.6477,2.5092,2.3826,2.2662,2.1587,2.0589,1.9660,1.879,1.7973,1.7203,1.6473};
  61   Int_t *index=new Int_t[nvectors];
  62   TMath::Sort(nvectors, data, index, kFALSE);
  63
  64   Int_t    nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11);
  65   Double_t factor = faclts[TMath::Max(0,nquant-1)];
  66
  67   Double_t sumx  =0;
  68   Double_t sumx2 =0;
  69   Int_t    bestindex = -1;
  70   Double_t bestmean  = 0;
  71   Double_t bestsigma = data[index[nvectors-1]]-data[index[0]];   // maximal possible sigma
  72   for (Int_t i=0; i<hh; i++){
  73     sumx  += data[index[i]];
  74     sumx2 += data[index[i]]*data[index[i]];
  75   }
  76
  77   Double_t norm = 1./Double_t(hh);
  78   Double_t norm2 = 1./Double_t(hh-1);
  79   for (Int_t i=hh; i<nvectors; i++){
  80     Double_t cmean  = sumx*norm;
  81     Double_t csigma = (sumx2 - hh*cmean*cmean)*norm2;
  82     if (csigma<bestsigma){
  83       bestmean  = cmean;
  84       bestsigma = csigma;
  85       bestindex = i-hh;
  86     }
  87
  88     sumx  += data[index[i]]-data[index[i-hh]];
  89     sumx2 += data[index[i]]*data[index[i]]-data[index[i-hh]]*data[index[i-hh]];
  90   }
  91
  92   Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma));
  93   mean  = bestmean;
  94   sigma = bstd;
  95   delete [] index;
  96
  97 }
  98
  99
 100
 101 void AliMathBase::EvaluateUniExternal(Int_t nvectors, Double_t *data, Double_t &mean, Double_t &sigma, Int_t hh,  Float_t externalfactor)
 102 {
 103   // Modified version of ROOT robust EvaluateUni
 104   // robust estimator in 1D case MI version
 105   // added external factor to include precision of external measurement
 106   //
 107
 108   if (hh==0)
 109     hh=(nvectors+2)/2;
 110   Double_t faclts[]={2.6477,2.5092,2.3826,2.2662,2.1587,2.0589,1.9660,1.879,1.7973,1.7203,1.6473};
 111   Int_t *index=new Int_t[nvectors];
 112   TMath::Sort(nvectors, data, index, kFALSE);
 113   //
 114   Int_t    nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11);
 115   Double_t factor = faclts[0];
 116   if (nquant>0){
 117     // fix proper normalization - Anja
 118     factor = faclts[nquant-1];
 119   }
 120
 121   //
 122   //
 123   Double_t sumx  =0;
 124   Double_t sumx2 =0;
 125   Int_t    bestindex = -1;
 126   Double_t bestmean  = 0;
 127   Double_t bestsigma = -1;
 128   for (Int_t i=0; i<hh; i++){
 129     sumx  += data[index[i]];
 130     sumx2 += data[index[i]]*data[index[i]];
 131   }
 132   //
 133   Double_t kfactor = 2.*externalfactor - externalfactor*externalfactor;
 134   Double_t norm = 1./Double_t(hh);
 135   for (Int_t i=hh; i<nvectors; i++){
 136     Double_t cmean  = sumx*norm;
 137     Double_t csigma = (sumx2*norm - cmean*cmean*kfactor);
 138     if (csigma<bestsigma ||  bestsigma<0){
 139       bestmean  = cmean;
 140       bestsigma = csigma;
 141       bestindex = i-hh;
 142     }
 143     //
 144     //
 145     sumx  += data[index[i]]-data[index[i-hh]];
 146     sumx2 += data[index[i]]*data[index[i]]-data[index[i-hh]]*data[index[i-hh]];
 147   }
 148
 149   Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma));
 150   mean  = bestmean;
 151   sigma = bstd;
 152   delete [] index;
 153 }
 154
 155
 156 //_____________________________________________________________________________
 157 Int_t AliMathBase::Freq(Int_t n, const Int_t *inlist
 158                         , Int_t *outlist, Bool_t down)
 159 {
 160   //
 161   //  Sort eleements according occurancy
 162   //  The size of output array has is 2*n
 163   //
 164
 165   Int_t * sindexS = new Int_t[n];     // temp array for sorting
 166   Int_t * sindexF = new Int_t[2*n];
 167   for (Int_t i=0;i<n;i++) sindexF[i]=0;
 168   //
 169   TMath::Sort(n,inlist, sindexS, down);
 170   Int_t last      = inlist[sindexS[0]];
 171   Int_t val       = last;
 172   sindexF[0]      = 1;
 173   sindexF[0+n]    = last;
 174   Int_t countPos  = 0;
 175   //
 176   //  find frequency
 177   for(Int_t i=1;i<n; i++){
 178     val = inlist[sindexS[i]];
 179     if (last == val)   sindexF[countPos]++;
 180     else{
 181       countPos++;
 182       sindexF[countPos+n] = val;
 183       sindexF[countPos]++;
 184       last =val;
 185     }
 186   }
 187   if (last==val) countPos++;
 188   // sort according frequency
 189   TMath::Sort(countPos, sindexF, sindexS, kTRUE);
 190   for (Int_t i=0;i<countPos;i++){
 191     outlist[2*i  ] = sindexF[sindexS[i]+n];
 192     outlist[2*i+1] = sindexF[sindexS[i]];
 193   }
 194   delete [] sindexS;
 195   delete [] sindexF;
 196
 197   return countPos;
 198
 199 }
 200
 201 //___AliMathBase__________________________________________________________________________
 202 void AliMathBase::TruncatedMean(TH1F * his, TVectorD *param, Float_t down, Float_t up, Bool_t verbose){
 203   //
 204   //
 205   //
 206   Int_t nbins    = his->GetNbinsX();
 207   Float_t nentries = his->GetEntries();
 208   Float_t sum      =0;
 209   Float_t mean   = 0;
 210   Float_t sigma2 = 0;
 211   Float_t ncumul=0;
 212   for (Int_t ibin=1;ibin<nbins; ibin++){
 213     ncumul+= his->GetBinContent(ibin);
 214     Float_t fraction = Float_t(ncumul)/Float_t(nentries);
 215     if (fraction>down && fraction<up){
 216       sum+=his->GetBinContent(ibin);
 217       mean+=his->GetBinCenter(ibin)*his->GetBinContent(ibin);
 218       sigma2+=his->GetBinCenter(ibin)*his->GetBinCenter(ibin)*his->GetBinContent(ibin);
 219     }
 220   }
 221   mean/=sum;
 222   sigma2= TMath::Sqrt(TMath::Abs(sigma2/sum-mean*mean));
 223   if (param){
 224     (*param)[0] = his->GetMaximum();
 225     (*param)[1] = mean;
 226     (*param)[2] = sigma2;
 227
 228   }
 229   if (verbose)  printf("Mean\t%f\t Sigma2\t%f\n", mean,sigma2);
 230 }
 231
 232 void AliMathBase::LTM(TH1F * his, TVectorD *param , Float_t fraction,  Bool_t verbose){
 233   //
 234   // LTM
 235   //
 236   Int_t nbins    = his->GetNbinsX();
 237   Int_t nentries = (Int_t)his->GetEntries();
 238   Double_t *data  = new Double_t[nentries];
 239   Int_t npoints=0;
 240   for (Int_t ibin=1;ibin<nbins; ibin++){
 241     Float_t entriesI = his->GetBinContent(ibin);
 242     Float_t xcenter= his->GetBinCenter(ibin);
 243     for (Int_t ic=0; ic<entriesI; ic++){
 244       if (npoints<nentries){
 245         data[npoints]= xcenter;
 246         npoints++;
 247       }
 248     }
 249   }
 250   Double_t mean, sigma;
 251   Int_t npoints2=TMath::Min(Int_t(fraction*Float_t(npoints)),npoints-1);
 252   npoints2=TMath::Max(Int_t(0.5*Float_t(npoints)),npoints2);
 253   AliMathBase::EvaluateUni(npoints, data, mean,sigma,npoints2);
 254   delete [] data;
 255   if (verbose)  printf("Mean\t%f\t Sigma2\t%f\n", mean,sigma);if (param){
 256     (*param)[0] = his->GetMaximum();
 257     (*param)[1] = mean;
 258     (*param)[2] = sigma;
 259   }
 260 }
 261
 262 Double_t  AliMathBase::FitGaus(TH1F* his, TVectorD *param, TMatrixD *matrix, Float_t xmin, Float_t xmax, Bool_t verbose){
 263   //
 264   //  Fit histogram with gaussian function
 265   //
 266   //  Prameters:
 267   //       return value- chi2 - if negative ( not enough points)
 268   //       his        -  input histogram
 269   //       param      -  vector with parameters
 270   //       xmin, xmax -  range to fit - if xmin=xmax=0 - the full histogram range used
 271   //  Fitting:
 272   //  1. Step - make logarithm
 273   //  2. Linear  fit (parabola) - more robust - always converge
 274   //  3. In case of small statistic bins are averaged
 275   //
 276   static TLinearFitter fitter(3,"pol2");
 277   TVectorD  par(3);
 278   TVectorD  sigma(3);
 279   TMatrixD mat(3,3);
 280   if (his->GetMaximum()<4) return -1;
 281   if (his->GetEntries()<12) return -1;
 282   if (his->GetRMS()<mat.GetTol()) return -1;
 283   Float_t maxEstimate   = his->GetEntries()*his->GetBinWidth(1)/TMath::Sqrt((2.*TMath::Pi()*his->GetRMS()));
 284   Int_t dsmooth = TMath::Nint(6./TMath::Sqrt(maxEstimate));
 285
 286   if (maxEstimate<1) return -1;
 287   Int_t nbins    = his->GetNbinsX();
 288   Int_t npoints=0;
 289   //
 290
 291
 292   if (xmin>=xmax){
 293     xmin = his->GetXaxis()->GetXmin();
 294     xmax = his->GetXaxis()->GetXmax();
 295   }
 296   for (Int_t iter=0; iter<2; iter++){
 297     fitter.ClearPoints();
 298     npoints=0;
 299     for (Int_t ibin=1;ibin<nbins-1; ibin++){
 300       Int_t countB=1;
 301       Float_t entriesI =  his->GetBinContent(ibin);
 302       for (Int_t delta = -dsmooth; delta<=dsmooth; delta++){
 303         if (ibin+delta>1 &&ibin+delta<nbins-1){
 304           entriesI +=  his->GetBinContent(ibin+delta);
 305           countB++;
 306         }
 307       }
 308       entriesI/=countB;
 309       Double_t xcenter= his->GetBinCenter(ibin);
 310       if (xcenter<xmin || xcenter>xmax) continue;
 311       Double_t error=1./TMath::Sqrt(countB);
 312       Float_t   cont=2;
 313       if (iter>0){
 314         if (par[0]+par[1]*xcenter+par[2]*xcenter*xcenter>20) return 0;
 315         cont = TMath::Exp(par[0]+par[1]*xcenter+par[2]*xcenter*xcenter);
 316         if (cont>1.) error = 1./TMath::Sqrt(cont*Float_t(countB));
 317       }
 318       if (entriesI>1&&cont>1){
 319         fitter.AddPoint(&xcenter,TMath::Log(Float_t(entriesI)),error);
 320         npoints++;
 321       }
 322     }
 323     if (npoints>3){
 324       fitter.Eval();
 325       fitter.GetParameters(par);
 326     }else{
 327       break;
 328     }
 329   }
 330   if (npoints<=3){
 331     return -1;
 332   }
 333   fitter.GetParameters(par);
 334   fitter.GetCovarianceMatrix(mat);
 335   if (TMath::Abs(par[1])<mat.GetTol()) return -1;
 336   if (TMath::Abs(par[2])<mat.GetTol()) return -1;
 337   Double_t chi2 = fitter.GetChisquare()/Float_t(npoints);
 338   //fitter.GetParameters();
 339   if (!param)  param  = new TVectorD(3);
 340   if (!matrix) matrix = new TMatrixD(3,3);
 341   (*param)[1] = par[1]/(-2.*par[2]);
 342   (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
 343   (*param)[0] = TMath::Exp(par[0]+ par[1]* (*param)[1] +  par[2]*(*param)[1]*(*param)[1]);
 344   if (verbose){
 345     par.Print();
 346     mat.Print();
 347     param->Print();
 348     printf("Chi2=%f\n",chi2);
 349     TF1 * f1= new TF1("f1","[0]*exp(-(x-[1])^2/(2*[2]*[2]))",his->GetXaxis()->GetXmin(),his->GetXaxis()->GetXmax());
 350     f1->SetParameter(0, (*param)[0]);
 351     f1->SetParameter(1, (*param)[1]);
 352     f1->SetParameter(2, (*param)[2]);
 353     f1->Draw("same");
 354   }
 355   return chi2;
 356 }
 357
 358