STAT/TStatToolkit.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 TStatToolkit
  19 //
  20 // Subset of  matheamtical functions  not included in the TMath
  21 //
  22 //
  23 /////////////////////////////////////////////////////////////////////////
  24 #include "TMath.h"
  25 #include "Riostream.h"
  26 #include "TH1F.h"
  27 #include "TH2F.h"
  28 #include "TH3.h"
  29 #include "TF1.h"
  30 #include "TTree.h"
  31 #include "TChain.h"
  32 #include "TObjString.h"
  33 #include "TLinearFitter.h"
  34 #include "TGraph2D.h"
  35 #include "TGraph.h"
  36 #include "TGraphErrors.h"
  37 #include "TMultiGraph.h"
  38 #include "TCanvas.h"
  39 #include "TLatex.h"
  40 #include "TCut.h"
  41 //
  42 // includes neccessary for test functions
  43 //
  44 #include "TSystem.h"
  45 #include "TRandom.h"
  46 #include "TStopwatch.h"
  47 #include "TTreeStream.h"
  48
  49 #include "TStatToolkit.h"
  50
  51
  52 ClassImp(TStatToolkit) // Class implementation to enable ROOT I/O
  53
  54 TStatToolkit::TStatToolkit() : TObject()
  55 {
  56   //
  57   // Default constructor
  58   //
  59 }
  60 ///////////////////////////////////////////////////////////////////////////
  61 TStatToolkit::~TStatToolkit()
  62 {
  63   //
  64   // Destructor
  65   //
  66 }
  67
  68
  69 //_____________________________________________________________________________
  70 void TStatToolkit::EvaluateUni(Int_t nvectors, Double_t *data, Double_t &mean
  71                            , Double_t &sigma, Int_t hh)
  72 {
  73   //
  74   // Robust estimator in 1D case MI version - (faster than ROOT version)
  75   //
  76   // For the univariate case
  77   // estimates of location and scatter are returned in mean and sigma parameters
  78   // the algorithm works on the same principle as in multivariate case -
  79   // it finds a subset of size hh with smallest sigma, and then returns mean and
  80   // sigma of this subset
  81   //
  82
  83   if (hh==0)
  84     hh=(nvectors+2)/2;
  85   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};
  86   Int_t *index=new Int_t[nvectors];
  87   TMath::Sort(nvectors, data, index, kFALSE);
  88
  89   Int_t    nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11);
  90   Double_t factor = faclts[TMath::Max(0,nquant-1)];
  91
  92   Double_t sumx  =0;
  93   Double_t sumx2 =0;
  94   Int_t    bestindex = -1;
  95   Double_t bestmean  = 0;
  96   Double_t bestsigma = (data[index[nvectors-1]]-data[index[0]]+1.);   // maximal possible sigma
  97   bestsigma *=bestsigma;
  98
  99   for (Int_t i=0; i<hh; i++){
 100     sumx  += data[index[i]];
 101     sumx2 += data[index[i]]*data[index[i]];
 102   }
 103
 104   Double_t norm = 1./Double_t(hh);
 105   Double_t norm2 = (hh-1)>0 ? 1./Double_t(hh-1):1;
 106   for (Int_t i=hh; i<nvectors; i++){
 107     Double_t cmean  = sumx*norm;
 108     Double_t csigma = (sumx2 - hh*cmean*cmean)*norm2;
 109     if (csigma<bestsigma){
 110       bestmean  = cmean;
 111       bestsigma = csigma;
 112       bestindex = i-hh;
 113     }
 114
 115     sumx  += data[index[i]]-data[index[i-hh]];
 116     sumx2 += data[index[i]]*data[index[i]]-data[index[i-hh]]*data[index[i-hh]];
 117   }
 118
 119   Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma));
 120   mean  = bestmean;
 121   sigma = bstd;
 122   delete [] index;
 123
 124 }
 125
 126
 127
 128 void TStatToolkit::EvaluateUniExternal(Int_t nvectors, Double_t *data, Double_t &mean, Double_t &sigma, Int_t hh,  Float_t externalfactor)
 129 {
 130   // Modified version of ROOT robust EvaluateUni
 131   // robust estimator in 1D case MI version
 132   // added external factor to include precision of external measurement
 133   //
 134
 135   if (hh==0)
 136     hh=(nvectors+2)/2;
 137   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};
 138   Int_t *index=new Int_t[nvectors];
 139   TMath::Sort(nvectors, data, index, kFALSE);
 140   //
 141   Int_t    nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11);
 142   Double_t factor = faclts[0];
 143   if (nquant>0){
 144     // fix proper normalization - Anja
 145     factor = faclts[nquant-1];
 146   }
 147
 148   //
 149   //
 150   Double_t sumx  =0;
 151   Double_t sumx2 =0;
 152   Int_t    bestindex = -1;
 153   Double_t bestmean  = 0;
 154   Double_t bestsigma = -1;
 155   for (Int_t i=0; i<hh; i++){
 156     sumx  += data[index[i]];
 157     sumx2 += data[index[i]]*data[index[i]];
 158   }
 159   //
 160   Double_t kfactor = 2.*externalfactor - externalfactor*externalfactor;
 161   Double_t norm = 1./Double_t(hh);
 162   for (Int_t i=hh; i<nvectors; i++){
 163     Double_t cmean  = sumx*norm;
 164     Double_t csigma = (sumx2*norm - cmean*cmean*kfactor);
 165     if (csigma<bestsigma ||  bestsigma<0){
 166       bestmean  = cmean;
 167       bestsigma = csigma;
 168       bestindex = i-hh;
 169     }
 170     //
 171     //
 172     sumx  += data[index[i]]-data[index[i-hh]];
 173     sumx2 += data[index[i]]*data[index[i]]-data[index[i-hh]]*data[index[i-hh]];
 174   }
 175
 176   Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma));
 177   mean  = bestmean;
 178   sigma = bstd;
 179   delete [] index;
 180 }
 181
 182
 183 //_____________________________________________________________________________
 184 Int_t TStatToolkit::Freq(Int_t n, const Int_t *inlist
 185                         , Int_t *outlist, Bool_t down)
 186 {
 187   //
 188   //  Sort eleements according occurancy
 189   //  The size of output array has is 2*n
 190   //
 191
 192   Int_t * sindexS = new Int_t[n];     // temp array for sorting
 193   Int_t * sindexF = new Int_t[2*n];
 194   for (Int_t i=0;i<n;i++) sindexS[i]=0;
 195   for (Int_t i=0;i<2*n;i++) sindexF[i]=0;
 196   //
 197   TMath::Sort(n,inlist, sindexS, down);
 198   Int_t last      = inlist[sindexS[0]];
 199   Int_t val       = last;
 200   sindexF[0]      = 1;
 201   sindexF[0+n]    = last;
 202   Int_t countPos  = 0;
 203   //
 204   //  find frequency
 205   for(Int_t i=1;i<n; i++){
 206     val = inlist[sindexS[i]];
 207     if (last == val)   sindexF[countPos]++;
 208     else{
 209       countPos++;
 210       sindexF[countPos+n] = val;
 211       sindexF[countPos]++;
 212       last =val;
 213     }
 214   }
 215   if (last==val) countPos++;
 216   // sort according frequency
 217   TMath::Sort(countPos, sindexF, sindexS, kTRUE);
 218   for (Int_t i=0;i<countPos;i++){
 219     outlist[2*i  ] = sindexF[sindexS[i]+n];
 220     outlist[2*i+1] = sindexF[sindexS[i]];
 221   }
 222   delete [] sindexS;
 223   delete [] sindexF;
 224
 225   return countPos;
 226
 227 }
 228
 229 //___TStatToolkit__________________________________________________________________________
 230 void TStatToolkit::TruncatedMean(const TH1 * his, TVectorD *param, Float_t down, Float_t up, Bool_t verbose){
 231   //
 232   //
 233   //
 234   Int_t nbins    = his->GetNbinsX();
 235   Float_t nentries = his->GetEntries();
 236   Float_t sum      =0;
 237   Float_t mean   = 0;
 238   Float_t sigma2 = 0;
 239   Float_t ncumul=0;
 240   for (Int_t ibin=1;ibin<nbins; ibin++){
 241     ncumul+= his->GetBinContent(ibin);
 242     Float_t fraction = Float_t(ncumul)/Float_t(nentries);
 243     if (fraction>down && fraction<up){
 244       sum+=his->GetBinContent(ibin);
 245       mean+=his->GetBinCenter(ibin)*his->GetBinContent(ibin);
 246       sigma2+=his->GetBinCenter(ibin)*his->GetBinCenter(ibin)*his->GetBinContent(ibin);
 247     }
 248   }
 249   mean/=sum;
 250   sigma2= TMath::Sqrt(TMath::Abs(sigma2/sum-mean*mean));
 251   if (param){
 252     (*param)[0] = his->GetMaximum();
 253     (*param)[1] = mean;
 254     (*param)[2] = sigma2;
 255
 256   }
 257   if (verbose)  printf("Mean\t%f\t Sigma2\t%f\n", mean,sigma2);
 258 }
 259
 260 void TStatToolkit::LTM(TH1 * his, TVectorD *param , Float_t fraction,  Bool_t verbose){
 261   //
 262   // LTM : Trimmed mean on histogram - Modified version for binned data
 263   //
 264   // Robust statistic to estimate properties of the distribution
 265   // See http://en.wikipedia.org/w/index.php?title=Trimmed_estimator&oldid=582847999
 266   //
 267   // New faster version is under preparation
 268   //
 269   if (!param) return;
 270   (*param)[0]=0;
 271   (*param)[1]=0;
 272   (*param)[2]=0;
 273   Int_t nbins    = his->GetNbinsX();
 274   Int_t nentries = (Int_t)his->GetEntries();
 275   if (nentries<=0) return;
 276   Double_t *data  = new Double_t[nentries];
 277   Int_t npoints=0;
 278   for (Int_t ibin=1;ibin<nbins; ibin++){
 279     Double_t entriesI = his->GetBinContent(ibin);
 280     //Double_t xcenter= his->GetBinCenter(ibin);
 281     Double_t x0 =  his->GetXaxis()->GetBinLowEdge(ibin);
 282     Double_t w  =  his->GetXaxis()->GetBinWidth(ibin);
 283     for (Int_t ic=0; ic<entriesI; ic++){
 284       if (npoints<nentries){
 285         data[npoints]= x0+w*Double_t((ic+0.5)/entriesI);
 286         npoints++;
 287       }
 288     }
 289   }
 290   Double_t mean, sigma;
 291   Int_t npoints2=TMath::Min(Int_t(fraction*Float_t(npoints)),npoints-1);
 292   npoints2=TMath::Max(Int_t(0.5*Float_t(npoints)),npoints2);
 293   TStatToolkit::EvaluateUni(npoints, data, mean,sigma,npoints2);
 294   delete [] data;
 295   if (verbose)  printf("Mean\t%f\t Sigma2\t%f\n", mean,sigma);if (param){
 296     (*param)[0] = his->GetMaximum();
 297     (*param)[1] = mean;
 298     (*param)[2] = sigma;
 299   }
 300 }
 301
 302
 303 void TStatToolkit::MedianFilter(TH1 * his1D, Int_t nmedian){
 304   //
 305   // Algorithm to filter  histogram
 306   // author:  marian.ivanov@cern.ch
 307   // Details of algorithm:
 308   // http://en.wikipedia.org/w/index.php?title=Median_filter&oldid=582191524
 309   // Input parameters:
 310   //    his1D - input histogam - to be modiefied by Medianfilter
 311   //    nmendian - number of bins in median filter
 312   //
 313   Int_t nbins    = his1D->GetNbinsX();
 314   TVectorD vectorH(nbins);
 315   for (Int_t ibin=0; ibin<nbins; ibin++) vectorH[ibin]=his1D->GetBinContent(ibin+1);
 316   for (Int_t ibin=0; ibin<nbins; ibin++) {
 317     Int_t index0=ibin-nmedian;
 318     Int_t index1=ibin+nmedian;
 319     if (index0<0) {index1+=-index0; index0=0;}
 320     if (index1>=nbins) {index0-=index1-nbins+1; index1=nbins-1;}
 321     Double_t value= TMath::Median(index1-index0,&(vectorH.GetMatrixArray()[index0]));
 322     his1D->SetBinContent(ibin+1, value);
 323   }
 324 }
 325
 326 Bool_t TStatToolkit::LTMHisto(TH1 *his1D, TVectorD &params , Float_t fraction){
 327   //
 328   // LTM : Trimmed mean on histogram - Modified version for binned data
 329   //
 330   // Robust statistic to estimate properties of the distribution
 331   // To handle binning error special treatment
 332   // for definition of unbinned data see:
 333   //     http://en.wikipedia.org/w/index.php?title=Trimmed_estimator&oldid=582847999
 334   //
 335   // Function parameters:
 336   //     his1D   - input histogram
 337   //     params  - vector with parameters
 338   //             - 0 - area
 339   //             - 1 - mean
 340   //             - 2 - rms
 341   //             - 3 - error estimate of mean
 342   //             - 4 - error estimate of RMS
 343   //             - 5 - first accepted bin position
 344   //             - 6 - last accepted  bin position
 345   //
 346   Int_t nbins    = his1D->GetNbinsX();
 347   Int_t nentries = (Int_t)his1D->GetEntries();
 348   const Double_t kEpsilon=0.0000000001;
 349
 350   if (nentries<=0) return 0;
 351   if (fraction>1) fraction=0;
 352   if (fraction<0) return 0;
 353   TVectorD vectorX(nbins);
 354   TVectorD vectorMean(nbins);
 355   TVectorD vectorRMS(nbins);
 356   Double_t sumCont=0;
 357   for (Int_t ibin0=1; ibin0<=nbins; ibin0++) sumCont+=his1D->GetBinContent(ibin0);
 358   //
 359   Double_t minRMS=his1D->GetRMS()*10000;
 360   Int_t maxBin=0;
 361   //
 362   for (Int_t ibin0=1; ibin0<nbins; ibin0++){
 363     Double_t sum0=0, sum1=0, sum2=0;
 364     Int_t ibin1=ibin0;
 365     for ( ibin1=ibin0; ibin1<nbins; ibin1++){
 366       Double_t cont=his1D->GetBinContent(ibin1);
 367       Double_t x= his1D->GetBinCenter(ibin1);
 368       sum0+=cont;
 369       sum1+=cont*x;
 370       sum2+=cont*x*x;
 371       if ( (ibin0!=ibin1) && sum0>=fraction*sumCont) break;
 372     }
 373     vectorX[ibin0]=his1D->GetBinCenter(ibin0);
 374     if (sum0<fraction*sumCont) continue;
 375     //
 376     // substract fractions of bin0 and bin1 to keep sum0=fration*sumCont
 377     //
 378     Double_t diff = sum0-fraction*sumCont;
 379     Double_t mean = sum1/sum0;
 380     //
 381     Double_t x0=his1D->GetBinCenter(ibin0);
 382     Double_t x1=his1D->GetBinCenter(ibin1);
 383     Double_t y0=his1D->GetBinContent(ibin0);
 384     Double_t y1=his1D->GetBinContent(ibin1);
 385     //
 386     Double_t d = y0+y1-diff;    //enties to keep
 387     Double_t w0=0,w1=0;
 388     if (y0<=kEpsilon&&y1>kEpsilon){
 389       w1=d/y1;
 390     }
 391     if (y1<=kEpsilon&&y0>kEpsilon){
 392       w0=d/y0;
 393     }
 394     if (y0>kEpsilon && y1>kEpsilon && x1>x0  ){
 395       w0 = (d*(x1-mean))/((x1-x0)*y0);
 396       w1 = (d-y0*w0)/y1;
 397       //
 398       if (w0>1) {w1+=(w0-1)*y0/y1; w0=1;}
 399       if (w1>1) {w0+=(w1-1)*y1/y0; w1=1;}
 400     }
 401     if ( (x1>x0) &&TMath::Abs(y0*w0+y1*w1-d)>kEpsilon*sum0){
 402       printf(" TStatToolkit::LTMHisto error\n");
 403     }
 404     sum0-=y0+y1;
 405     sum1-=y0*x0;
 406     sum1-=y1*x1;
 407     sum2-=y0*x0*x0;
 408     sum2-=y1*x1*x1;
 409     //
 410     Double_t xx0=his1D->GetXaxis()->GetBinUpEdge(ibin0)-0.5*w0*his1D->GetBinWidth(ibin0);
 411     Double_t xx1=his1D->GetXaxis()->GetBinLowEdge(ibin1)+0.5*w1*his1D->GetBinWidth(ibin1);
 412     sum0+=y0*w0+y1*w1;
 413     sum1+=y0*w0*xx0;
 414     sum1+=y1*w1*xx1;
 415     sum2+=y0*w0*xx0*xx0;
 416     sum2+=y1*w1*xx1*xx1;
 417
 418     //
 419     // choose the bin with smallest rms
 420     //
 421     if (sum0>0){
 422       vectorMean[ibin0]=sum1/sum0;
 423       vectorRMS[ibin0]=TMath::Sqrt(TMath::Abs(sum2/sum0-vectorMean[ibin0]*vectorMean[ibin0]));
 424       if (vectorRMS[ibin0]<minRMS){
 425         minRMS=vectorRMS[ibin0];
 426         params[0]=sum0;
 427         params[1]=vectorMean[ibin0];
 428         params[2]=vectorRMS[ibin0];
 429         params[3]=vectorRMS[ibin0]/TMath::Sqrt(sumCont*fraction);
 430         params[4]=0; // what is the formula for error of RMS???
 431         params[5]=ibin0;
 432         params[6]=ibin1;
 433         params[7]=his1D->GetBinCenter(ibin0);
 434         params[8]=his1D->GetBinCenter(ibin1);
 435         maxBin=ibin0;
 436       }
 437     }else{
 438       break;
 439     }
 440   }
 441   return kTRUE;
 442 }
 443
 444
 445 Double_t  TStatToolkit::FitGaus(TH1* his, TVectorD *param, TMatrixD */*matrix*/, Float_t xmin, Float_t xmax, Bool_t verbose){
 446   //
 447   //  Fit histogram with gaussian function
 448   //
 449   //  Prameters:
 450   //       return value- chi2 - if negative ( not enough points)
 451   //       his        -  input histogram
 452   //       param      -  vector with parameters
 453   //       xmin, xmax -  range to fit - if xmin=xmax=0 - the full histogram range used
 454   //  Fitting:
 455   //  1. Step - make logarithm
 456   //  2. Linear  fit (parabola) - more robust - always converge
 457   //  3. In case of small statistic bins are averaged
 458   //
 459   static TLinearFitter fitter(3,"pol2");
 460   TVectorD  par(3);
 461   TVectorD  sigma(3);
 462   TMatrixD mat(3,3);
 463   if (his->GetMaximum()<4) return -1;
 464   if (his->GetEntries()<12) return -1;
 465   if (his->GetRMS()<mat.GetTol()) return -1;
 466   Float_t maxEstimate   = his->GetEntries()*his->GetBinWidth(1)/TMath::Sqrt((TMath::TwoPi()*his->GetRMS()));
 467   Int_t dsmooth = TMath::Nint(6./TMath::Sqrt(maxEstimate));
 468
 469   if (maxEstimate<1) return -1;
 470   Int_t nbins    = his->GetNbinsX();
 471   Int_t npoints=0;
 472   //
 473
 474
 475   if (xmin>=xmax){
 476     xmin = his->GetXaxis()->GetXmin();
 477     xmax = his->GetXaxis()->GetXmax();
 478   }
 479   for (Int_t iter=0; iter<2; iter++){
 480     fitter.ClearPoints();
 481     npoints=0;
 482     for (Int_t ibin=1;ibin<nbins+1; ibin++){
 483       Int_t countB=1;
 484       Float_t entriesI =  his->GetBinContent(ibin);
 485       for (Int_t delta = -dsmooth; delta<=dsmooth; delta++){
 486         if (ibin+delta>1 &&ibin+delta<nbins-1){
 487           entriesI +=  his->GetBinContent(ibin+delta);
 488           countB++;
 489         }
 490       }
 491       entriesI/=countB;
 492       Double_t xcenter= his->GetBinCenter(ibin);
 493       if (xcenter<xmin || xcenter>xmax) continue;
 494       Double_t error=1./TMath::Sqrt(countB);
 495       Float_t   cont=2;
 496       if (iter>0){
 497         if (par[0]+par[1]*xcenter+par[2]*xcenter*xcenter>20) return 0;
 498         cont = TMath::Exp(par[0]+par[1]*xcenter+par[2]*xcenter*xcenter);
 499         if (cont>1.) error = 1./TMath::Sqrt(cont*Float_t(countB));
 500       }
 501       if (entriesI>1&&cont>1){
 502         fitter.AddPoint(&xcenter,TMath::Log(Float_t(entriesI)),error);
 503         npoints++;
 504       }
 505     }
 506     if (npoints>3){
 507       fitter.Eval();
 508       fitter.GetParameters(par);
 509     }else{
 510       break;
 511     }
 512   }
 513   if (npoints<=3){
 514     return -1;
 515   }
 516   fitter.GetParameters(par);
 517   fitter.GetCovarianceMatrix(mat);
 518   if (TMath::Abs(par[1])<mat.GetTol()) return -1;
 519   if (TMath::Abs(par[2])<mat.GetTol()) return -1;
 520   Double_t chi2 = fitter.GetChisquare()/Float_t(npoints);
 521   //fitter.GetParameters();
 522   if (!param)  param  = new TVectorD(3);
 523   // if (!matrix) matrix = new TMatrixD(3,3); // Covariance matrix to be implemented
 524   (*param)[1] = par[1]/(-2.*par[2]);
 525   (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
 526   (*param)[0] = TMath::Exp(par[0]+ par[1]* (*param)[1] +  par[2]*(*param)[1]*(*param)[1]);
 527   if (verbose){
 528     par.Print();
 529     mat.Print();
 530     param->Print();
 531     printf("Chi2=%f\n",chi2);
 532     TF1 * f1= new TF1("f1","[0]*exp(-(x-[1])^2/(2*[2]*[2]))",his->GetXaxis()->GetXmin(),his->GetXaxis()->GetXmax());
 533     f1->SetParameter(0, (*param)[0]);
 534     f1->SetParameter(1, (*param)[1]);
 535     f1->SetParameter(2, (*param)[2]);
 536     f1->Draw("same");
 537   }
 538   return chi2;
 539 }
 540
 541 Double_t  TStatToolkit::FitGaus(Float_t *arr, Int_t nBins, Float_t xMin, Float_t xMax, TVectorD *param, TMatrixD */*matrix*/, Bool_t verbose){
 542   //
 543   //  Fit histogram with gaussian function
 544   //
 545   //  Prameters:
 546   //     nbins: size of the array and number of histogram bins
 547   //     xMin, xMax: histogram range
 548   //     param: paramters of the fit (0-Constant, 1-Mean, 2-Sigma)
 549   //     matrix: covariance matrix -- not implemented yet, pass dummy matrix!!!
 550   //
 551   //  Return values:
 552   //    >0: the chi2 returned by TLinearFitter
 553   //    -3: only three points have been used for the calculation - no fitter was used
 554   //    -2: only two points have been used for the calculation - center of gravity was uesed for calculation
 555   //    -1: only one point has been used for the calculation - center of gravity was uesed for calculation
 556   //    -4: invalid result!!
 557   //
 558   //  Fitting:
 559   //  1. Step - make logarithm
 560   //  2. Linear  fit (parabola) - more robust - always converge
 561   //
 562   static TLinearFitter fitter(3,"pol2");
 563   static TMatrixD mat(3,3);
 564   static Double_t kTol = mat.GetTol();
 565   fitter.StoreData(kFALSE);
 566   fitter.ClearPoints();
 567   TVectorD  par(3);
 568   TVectorD  sigma(3);
 569   TMatrixD matA(3,3);
 570   TMatrixD b(3,1);
 571   Float_t rms = TMath::RMS(nBins,arr);
 572   Float_t max = TMath::MaxElement(nBins,arr);
 573   Float_t binWidth = (xMax-xMin)/(Float_t)nBins;
 574
 575   Float_t meanCOG = 0;
 576   Float_t rms2COG = 0;
 577   Float_t sumCOG  = 0;
 578
 579   Float_t entries = 0;
 580   Int_t nfilled=0;
 581
 582   for (Int_t i=0; i<nBins; i++){
 583       entries+=arr[i];
 584       if (arr[i]>0) nfilled++;
 585   }
 586
 587   if (max<4) return -4;
 588   if (entries<12) return -4;
 589   if (rms<kTol) return -4;
 590
 591   Int_t npoints=0;
 592   //
 593
 594   //
 595   for (Int_t ibin=0;ibin<nBins; ibin++){
 596       Float_t entriesI = arr[ibin];
 597     if (entriesI>1){
 598       Double_t xcenter = xMin+(ibin+0.5)*binWidth;
 599
 600       Float_t error    = 1./TMath::Sqrt(entriesI);
 601       Float_t val = TMath::Log(Float_t(entriesI));
 602       fitter.AddPoint(&xcenter,val,error);
 603       if (npoints<3){
 604           matA(npoints,0)=1;
 605           matA(npoints,1)=xcenter;
 606           matA(npoints,2)=xcenter*xcenter;
 607           b(npoints,0)=val;
 608           meanCOG+=xcenter*entriesI;
 609           rms2COG +=xcenter*entriesI*xcenter;
 610           sumCOG +=entriesI;
 611       }
 612       npoints++;
 613     }
 614   }
 615
 616
 617   Double_t chi2 = 0;
 618   if (npoints>=3){
 619       if ( npoints == 3 ){
 620           //analytic calculation of the parameters for three points
 621           matA.Invert();
 622           TMatrixD res(1,3);
 623           res.Mult(matA,b);
 624           par[0]=res(0,0);
 625           par[1]=res(0,1);
 626           par[2]=res(0,2);
 627           chi2 = -3.;
 628       } else {
 629           // use fitter for more than three points
 630           fitter.Eval();
 631           fitter.GetParameters(par);
 632           fitter.GetCovarianceMatrix(mat);
 633           chi2 = fitter.GetChisquare()/Float_t(npoints);
 634       }
 635       if (TMath::Abs(par[1])<kTol) return -4;
 636       if (TMath::Abs(par[2])<kTol) return -4;
 637
 638       if (!param)  param  = new TVectorD(3);
 639       //if (!matrix) matrix = new TMatrixD(3,3);  // !!!!might be a memory leek. use dummy matrix pointer to call this function! // Covariance matrix to be implemented
 640
 641       (*param)[1] = par[1]/(-2.*par[2]);
 642       (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
 643       Double_t lnparam0 = par[0]+ par[1]* (*param)[1] +  par[2]*(*param)[1]*(*param)[1];
 644       if ( lnparam0>307 ) return -4;
 645       (*param)[0] = TMath::Exp(lnparam0);
 646       if (verbose){
 647           par.Print();
 648           mat.Print();
 649           param->Print();
 650           printf("Chi2=%f\n",chi2);
 651           TF1 * f1= new TF1("f1","[0]*exp(-(x-[1])^2/(2*[2]*[2]))",xMin,xMax);
 652           f1->SetParameter(0, (*param)[0]);
 653           f1->SetParameter(1, (*param)[1]);
 654           f1->SetParameter(2, (*param)[2]);
 655           f1->Draw("same");
 656       }
 657       return chi2;
 658   }
 659
 660   if (npoints == 2){
 661       //use center of gravity for 2 points
 662       meanCOG/=sumCOG;
 663       rms2COG /=sumCOG;
 664       (*param)[0] = max;
 665       (*param)[1] = meanCOG;
 666       (*param)[2] = TMath::Sqrt(TMath::Abs(meanCOG*meanCOG-rms2COG));
 667       chi2=-2.;
 668   }
 669   if ( npoints == 1 ){
 670       meanCOG/=sumCOG;
 671       (*param)[0] = max;
 672       (*param)[1] = meanCOG;
 673       (*param)[2] = binWidth/TMath::Sqrt(12);
 674       chi2=-1.;
 675   }
 676   return chi2;
 677
 678 }
 679
 680
 681 Float_t TStatToolkit::GetCOG(const Short_t *arr, Int_t nBins, Float_t xMin, Float_t xMax, Float_t *rms, Float_t *sum)
 682 {
 683     //
 684     //  calculate center of gravity rms and sum for array 'arr' with nBins an a x range xMin to xMax
 685     //  return COG; in case of failure return xMin
 686     //
 687     Float_t meanCOG = 0;
 688     Float_t rms2COG = 0;
 689     Float_t sumCOG  = 0;
 690     Int_t npoints   = 0;
 691
 692     Float_t binWidth = (xMax-xMin)/(Float_t)nBins;
 693
 694     for (Int_t ibin=0; ibin<nBins; ibin++){
 695         Float_t entriesI = (Float_t)arr[ibin];
 696         Double_t xcenter = xMin+(ibin+0.5)*binWidth;
 697         if ( entriesI>0 ){
 698             meanCOG += xcenter*entriesI;
 699             rms2COG += xcenter*entriesI*xcenter;
 700             sumCOG  += entriesI;
 701             npoints++;
 702         }
 703     }
 704     if ( sumCOG == 0 ) return xMin;
 705     meanCOG/=sumCOG;
 706
 707     if ( rms ){
 708         rms2COG /=sumCOG;
 709         (*rms) = TMath::Sqrt(TMath::Abs(meanCOG*meanCOG-rms2COG));
 710         if ( npoints == 1 ) (*rms) = binWidth/TMath::Sqrt(12);
 711     }
 712
 713     if ( sum )
 714         (*sum) = sumCOG;
 715
 716     return meanCOG;
 717 }
 718
 719
 720
 721 ///////////////////////////////////////////////////////////////
 722 //////////////         TEST functions /////////////////////////
 723 ///////////////////////////////////////////////////////////////
 724
 725
 726
 727
 728
 729 void TStatToolkit::TestGausFit(Int_t nhistos){
 730   //
 731   // Test performance of the parabolic - gaussian fit - compare it with
 732   // ROOT gauss fit
 733   //  nhistos - number of histograms to be used for test
 734   //
 735   TTreeSRedirector *pcstream = new TTreeSRedirector("fitdebug.root","recreate");
 736
 737   Float_t  *xTrue = new Float_t[nhistos];
 738   Float_t  *sTrue = new Float_t[nhistos];
 739   TVectorD **par1  = new TVectorD*[nhistos];
 740   TVectorD **par2  = new TVectorD*[nhistos];
 741   TMatrixD dummy(3,3);
 742
 743
 744   TH1F **h1f = new TH1F*[nhistos];
 745   TF1  *myg = new TF1("myg","gaus");
 746   TF1  *fit = new TF1("fit","gaus");
 747   gRandom->SetSeed(0);
 748
 749   //init
 750   for (Int_t i=0;i<nhistos; i++){
 751     par1[i] = new TVectorD(3);
 752     par2[i] = new TVectorD(3);
 753     h1f[i]  = new TH1F(Form("h1f%d",i),Form("h1f%d",i),20,-10,10);
 754     xTrue[i]= gRandom->Rndm();
 755     gSystem->Sleep(2);
 756     sTrue[i]= .75+gRandom->Rndm()*.5;
 757     myg->SetParameters(1,xTrue[i],sTrue[i]);
 758     h1f[i]->FillRandom("myg");
 759   }
 760
 761   TStopwatch s;
 762   s.Start();
 763   //standard gaus fit
 764   for (Int_t i=0; i<nhistos; i++){
 765     h1f[i]->Fit(fit,"0q");
 766     (*par1[i])(0) = fit->GetParameter(0);
 767     (*par1[i])(1) = fit->GetParameter(1);
 768     (*par1[i])(2) = fit->GetParameter(2);
 769   }
 770   s.Stop();
 771   printf("Gaussian fit\t");
 772   s.Print();
 773
 774   s.Start();
 775   //TStatToolkit gaus fit
 776   for (Int_t i=0; i<nhistos; i++){
 777     TStatToolkit::FitGaus(h1f[i]->GetArray()+1,h1f[i]->GetNbinsX(),h1f[i]->GetXaxis()->GetXmin(),h1f[i]->GetXaxis()->GetXmax(),par2[i],&dummy);
 778   }
 779
 780   s.Stop();
 781   printf("Parabolic fit\t");
 782   s.Print();
 783   //write stream
 784   for (Int_t i=0;i<nhistos; i++){
 785     Float_t xt  = xTrue[i];
 786     Float_t st  = sTrue[i];
 787     (*pcstream)<<"data"
 788                <<"xTrue="<<xt
 789                <<"sTrue="<<st
 790                <<"pg.="<<(par1[i])
 791                <<"pa.="<<(par2[i])
 792                <<"\n";
 793   }
 794   //delete pointers
 795   for (Int_t i=0;i<nhistos; i++){
 796     delete par1[i];
 797     delete par2[i];
 798     delete h1f[i];
 799   }
 800   delete pcstream;
 801   delete []h1f;
 802   delete []xTrue;
 803   delete []sTrue;
 804   //
 805   delete []par1;
 806   delete []par2;
 807
 808 }
 809
 810
 811
 812 TGraph2D * TStatToolkit::MakeStat2D(TH3 * his, Int_t delta0, Int_t delta1, Int_t type){
 813   //
 814   //
 815   //
 816   // delta - number of bins to integrate
 817   // type - 0 - mean value
 818
 819   TAxis * xaxis  = his->GetXaxis();
 820   TAxis * yaxis  = his->GetYaxis();
 821   //  TAxis * zaxis  = his->GetZaxis();
 822   Int_t   nbinx  = xaxis->GetNbins();
 823   Int_t   nbiny  = yaxis->GetNbins();
 824   char name[1000];
 825   Int_t icount=0;
 826   TGraph2D  *graph = new TGraph2D(nbinx*nbiny);
 827   TF1 f1("f1","gaus");
 828   for (Int_t ix=0; ix<nbinx;ix++)
 829     for (Int_t iy=0; iy<nbiny;iy++){
 830       Float_t xcenter = xaxis->GetBinCenter(ix);
 831       Float_t ycenter = yaxis->GetBinCenter(iy);
 832       snprintf(name,1000,"%s_%d_%d",his->GetName(), ix,iy);
 833       TH1 *projection = his->ProjectionZ(name,ix-delta0,ix+delta0,iy-delta1,iy+delta1);
 834       Float_t stat= 0;
 835       if (type==0) stat = projection->GetMean();
 836       if (type==1) stat = projection->GetRMS();
 837       if (type==2 || type==3){
 838         TVectorD vec(10);
 839         TStatToolkit::LTM((TH1F*)projection,&vec,0.7);
 840         if (type==2) stat= vec[1];
 841         if (type==3) stat= vec[0];
 842       }
 843       if (type==4|| type==5){
 844         projection->Fit(&f1);
 845         if (type==4) stat= f1.GetParameter(1);
 846         if (type==5) stat= f1.GetParameter(2);
 847       }
 848       //printf("%d\t%f\t%f\t%f\n", icount,xcenter, ycenter, stat);
 849       graph->SetPoint(icount,xcenter, ycenter, stat);
 850       icount++;
 851     }
 852   return graph;
 853 }
 854
 855 TGraphErrors * TStatToolkit::MakeStat1D(TH2 * his, Int_t deltaBin, Double_t fraction, Int_t returnType, Int_t markerStyle, Int_t markerColor){
 856   //
 857   // function to retrieve the "mean and RMS estimate" of 2D histograms
 858   //
 859   // Robust statistic to estimate properties of the distribution
 860   // See http://en.wikipedia.org/wiki/Trimmed_estimator
 861   //
 862   // deltaBin - number of bins to integrate (bin+-deltaBin)
 863   // fraction - fraction of values for the LTM and for the gauss fit
 864   // returnType -
 865   //        0 - mean value
 866   //        1 - RMS
 867   //        2 - LTM mean
 868   //        3 - LTM sigma
 869   //        4 - Gaus fit mean  - on LTM range
 870   //        5 - Gaus fit sigma - on LTM  range
 871   //
 872   TAxis * xaxis  = his->GetXaxis();
 873   Int_t   nbinx  = xaxis->GetNbins();
 874   char name[1000];
 875   Int_t icount=0;
 876   //
 877   TVectorD vecX(nbinx);
 878   TVectorD vecXErr(nbinx);
 879   TVectorD vecY(nbinx);
 880   TVectorD vecYErr(nbinx);
 881   //
 882   TF1 f1("f1","gaus");
 883   TVectorD vecLTM(10);
 884
 885   for (Int_t jx=1; jx<=nbinx;jx++){
 886     Int_t ix=jx-1;
 887     Float_t xcenter = xaxis->GetBinCenter(jx);
 888     snprintf(name,1000,"%s_%d",his->GetName(), ix);
 889     TH1 *projection = his->ProjectionY(name,TMath::Max(jx-deltaBin,1),TMath::Min(jx+deltaBin,nbinx));
 890     Double_t stat= 0;
 891     Double_t err =0;
 892     TStatToolkit::LTMHisto((TH1F*)projection,vecLTM,fraction);
 893     //
 894     if (returnType==0) {
 895       stat = projection->GetMean();
 896       err  = projection->GetMeanError();
 897     }
 898     if (returnType==1) {
 899       stat = projection->GetRMS();
 900       err = projection->GetRMSError();
 901     }
 902     if (returnType==2 || returnType==3){
 903       if (returnType==2) {stat= vecLTM[1];  err =projection->GetRMSError();}
 904         if (returnType==3) {stat= vecLTM[2];     err =projection->GetRMSError();}
 905     }
 906     if (returnType==4|| returnType==5){
 907       projection->Fit(&f1,"QN","QN", vecLTM[7], vecLTM[8]);
 908       if (returnType==4) {
 909         stat= f1.GetParameter(1);
 910         err=f1.GetParError(1);
 911       }
 912       if (returnType==5) {
 913         stat= f1.GetParameter(2);
 914         err=f1.GetParError(2);
 915       }
 916     }
 917     vecX[icount]=xcenter;
 918     vecY[icount]=stat;
 919     vecYErr[icount]=err;
 920     icount++;
 921     delete projection;
 922   }
 923   TGraphErrors  *graph = new TGraphErrors(icount,vecX.GetMatrixArray(), vecY.GetMatrixArray(),0, vecYErr.GetMatrixArray());
 924   graph->SetMarkerStyle(markerStyle);
 925   graph->SetMarkerColor(markerColor);
 926   return graph;
 927 }
 928
 929
 930
 931
 932
 933 TString* TStatToolkit::FitPlane(TTree *tree, const char* drawCommand, const char* formula, const char* cuts, Double_t & chi2, Int_t &npoints, TVectorD &fitParam, TMatrixD &covMatrix, Float_t frac, Int_t start, Int_t stop,Bool_t fix0){
 934    //
 935    // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts
 936    // returns chi2, fitParam and covMatrix
 937    // returns TString with fitted formula
 938    //
 939
 940    TString formulaStr(formula);
 941    TString drawStr(drawCommand);
 942    TString cutStr(cuts);
 943    TString ferr("1");
 944
 945    TString strVal(drawCommand);
 946    if (strVal.Contains(":")){
 947      TObjArray* valTokens = strVal.Tokenize(":");
 948      drawStr = valTokens->At(0)->GetName();
 949      ferr       = valTokens->At(1)->GetName();
 950      delete valTokens;
 951    }
 952
 953
 954    formulaStr.ReplaceAll("++", "~");
 955    TObjArray* formulaTokens = formulaStr.Tokenize("~");
 956    Int_t dim = formulaTokens->GetEntriesFast();
 957
 958    fitParam.ResizeTo(dim);
 959    covMatrix.ResizeTo(dim,dim);
 960
 961    TLinearFitter* fitter = new TLinearFitter(dim+1, Form("hyp%d",dim));
 962    fitter->StoreData(kTRUE);
 963    fitter->ClearPoints();
 964
 965    Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff",  stop-start, start);
 966    if (entries == -1) {
 967      delete formulaTokens;
 968      return new TString("An ERROR has occured during fitting!");
 969    }
 970    Double_t **values = new Double_t*[dim+1] ;
 971    for (Int_t i=0; i<dim+1; i++) values[i]=NULL;
 972    //
 973    entries = tree->Draw(ferr.Data(), cutStr.Data(), "goff",  stop-start, start);
 974    if (entries == -1) {
 975      delete formulaTokens;
 976      delete []values;
 977      return new TString("An ERROR has occured during fitting!");
 978    }
 979    Double_t *errors = new Double_t[entries];
 980    memcpy(errors,  tree->GetV1(), entries*sizeof(Double_t));
 981
 982    for (Int_t i = 0; i < dim + 1; i++){
 983       Int_t centries = 0;
 984       if (i < dim) centries = tree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start);
 985       else  centries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start);
 986
 987       if (entries != centries) {
 988         delete []errors;
 989         delete []values;
 990         return new TString("An ERROR has occured during fitting!");
 991       }
 992       values[i] = new Double_t[entries];
 993       memcpy(values[i],  tree->GetV1(), entries*sizeof(Double_t));
 994    }
 995
 996    // add points to the fitter
 997    for (Int_t i = 0; i < entries; i++){
 998       Double_t x[1000];
 999       for (Int_t j=0; j<dim;j++) x[j]=values[j][i];
1000       fitter->AddPoint(x, values[dim][i], errors[i]);
1001    }
1002
1003    fitter->Eval();
1004    if (frac>0.5 && frac<1){
1005      fitter->EvalRobust(frac);
1006    }else{
1007      if (fix0) {
1008        fitter->FixParameter(0,0);
1009        fitter->Eval();
1010      }
1011    }
1012    fitter->GetParameters(fitParam);
1013    fitter->GetCovarianceMatrix(covMatrix);
1014    chi2 = fitter->GetChisquare();
1015    npoints = entries;
1016    TString *preturnFormula = new TString(Form("( %f+",fitParam[0])), &returnFormula = *preturnFormula;
1017
1018    for (Int_t iparam = 0; iparam < dim; iparam++) {
1019      returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam+1]));
1020      if (iparam < dim-1) returnFormula.Append("+");
1021    }
1022    returnFormula.Append(" )");
1023
1024
1025    for (Int_t j=0; j<dim+1;j++) delete [] values[j];
1026
1027
1028    delete formulaTokens;
1029    delete fitter;
1030    delete[] values;
1031    delete[] errors;
1032    return preturnFormula;
1033 }
1034
1035 TString* TStatToolkit::FitPlaneConstrain(TTree *tree, const char* drawCommand, const char* formula, const char* cuts, Double_t & chi2, Int_t &npoints, TVectorD &fitParam, TMatrixD &covMatrix, Float_t frac, Int_t start, Int_t stop,Double_t constrain){
1036    //
1037    // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts
1038    // returns chi2, fitParam and covMatrix
1039    // returns TString with fitted formula
1040    //
1041
1042    TString formulaStr(formula);
1043    TString drawStr(drawCommand);
1044    TString cutStr(cuts);
1045    TString ferr("1");
1046
1047    TString strVal(drawCommand);
1048    if (strVal.Contains(":")){
1049      TObjArray* valTokens = strVal.Tokenize(":");
1050      drawStr = valTokens->At(0)->GetName();
1051      ferr       = valTokens->At(1)->GetName();
1052      delete valTokens;
1053    }
1054
1055
1056    formulaStr.ReplaceAll("++", "~");
1057    TObjArray* formulaTokens = formulaStr.Tokenize("~");
1058    Int_t dim = formulaTokens->GetEntriesFast();
1059
1060    fitParam.ResizeTo(dim);
1061    covMatrix.ResizeTo(dim,dim);
1062
1063    TLinearFitter* fitter = new TLinearFitter(dim+1, Form("hyp%d",dim));
1064    fitter->StoreData(kTRUE);
1065    fitter->ClearPoints();
1066
1067    Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff",  stop-start, start);
1068    if (entries == -1) {
1069      delete formulaTokens;
1070      return new TString("An ERROR has occured during fitting!");
1071    }
1072    Double_t **values = new Double_t*[dim+1] ;
1073    for (Int_t i=0; i<dim+1; i++) values[i]=NULL;
1074    //
1075    entries = tree->Draw(ferr.Data(), cutStr.Data(), "goff",  stop-start, start);
1076    if (entries == -1) {
1077      delete formulaTokens;
1078      delete [] values;
1079      return new TString("An ERROR has occured during fitting!");
1080    }
1081    Double_t *errors = new Double_t[entries];
1082    memcpy(errors,  tree->GetV1(), entries*sizeof(Double_t));
1083
1084    for (Int_t i = 0; i < dim + 1; i++){
1085       Int_t centries = 0;
1086       if (i < dim) centries = tree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start);
1087       else  centries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start);
1088
1089       if (entries != centries) {
1090         delete []errors;
1091         delete []values;
1092         delete formulaTokens;
1093         return new TString("An ERROR has occured during fitting!");
1094       }
1095       values[i] = new Double_t[entries];
1096       memcpy(values[i],  tree->GetV1(), entries*sizeof(Double_t));
1097    }
1098
1099    // add points to the fitter
1100    for (Int_t i = 0; i < entries; i++){
1101       Double_t x[1000];
1102       for (Int_t j=0; j<dim;j++) x[j]=values[j][i];
1103       fitter->AddPoint(x, values[dim][i], errors[i]);
1104    }
1105    if (constrain>0){
1106      for (Int_t i = 0; i < dim; i++){
1107        Double_t x[1000];
1108        for (Int_t j=0; j<dim;j++) if (i!=j) x[j]=0;
1109        x[i]=1.;
1110        fitter->AddPoint(x, 0, constrain);
1111      }
1112    }
1113
1114
1115    fitter->Eval();
1116    if (frac>0.5 && frac<1){
1117      fitter->EvalRobust(frac);
1118    }
1119    fitter->GetParameters(fitParam);
1120    fitter->GetCovarianceMatrix(covMatrix);
1121    chi2 = fitter->GetChisquare();
1122    npoints = entries;
1123
1124    TString *preturnFormula = new TString(Form("( %f+",fitParam[0])), &returnFormula = *preturnFormula;
1125
1126    for (Int_t iparam = 0; iparam < dim; iparam++) {
1127      returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam+1]));
1128      if (iparam < dim-1) returnFormula.Append("+");
1129    }
1130    returnFormula.Append(" )");
1131
1132    for (Int_t j=0; j<dim+1;j++) delete [] values[j];
1133
1134
1135
1136    delete formulaTokens;
1137    delete fitter;
1138    delete[] values;
1139    delete[] errors;
1140    return preturnFormula;
1141 }
1142
1143
1144
1145 TString* TStatToolkit::FitPlaneFixed(TTree *tree, const char* drawCommand, const char* formula, const char* cuts, Double_t & chi2, Int_t &npoints, TVectorD &fitParam, TMatrixD &covMatrix, Float_t frac, Int_t start, Int_t stop){
1146    //
1147    // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts
1148    // returns chi2, fitParam and covMatrix
1149    // returns TString with fitted formula
1150    //
1151
1152    TString formulaStr(formula);
1153    TString drawStr(drawCommand);
1154    TString cutStr(cuts);
1155    TString ferr("1");
1156
1157    TString strVal(drawCommand);
1158    if (strVal.Contains(":")){
1159      TObjArray* valTokens = strVal.Tokenize(":");
1160      drawStr = valTokens->At(0)->GetName();
1161      ferr       = valTokens->At(1)->GetName();
1162      delete valTokens;
1163    }
1164
1165
1166    formulaStr.ReplaceAll("++", "~");
1167    TObjArray* formulaTokens = formulaStr.Tokenize("~");
1168    Int_t dim = formulaTokens->GetEntriesFast();
1169
1170    fitParam.ResizeTo(dim);
1171    covMatrix.ResizeTo(dim,dim);
1172    TString fitString="x0";
1173    for (Int_t i=1; i<dim; i++) fitString+=Form("++x%d",i);
1174    TLinearFitter* fitter = new TLinearFitter(dim, fitString.Data());
1175    fitter->StoreData(kTRUE);
1176    fitter->ClearPoints();
1177
1178    Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff",  stop-start, start);
1179    if (entries == -1) {
1180      delete formulaTokens;
1181      return new TString("An ERROR has occured during fitting!");
1182    }
1183    Double_t **values = new Double_t*[dim+1] ;
1184    for (Int_t i=0; i<dim+1; i++) values[i]=NULL;
1185    //
1186    entries = tree->Draw(ferr.Data(), cutStr.Data(), "goff",  stop-start, start);
1187    if (entries == -1) {
1188      delete []values;
1189      delete formulaTokens;
1190      return new TString("An ERROR has occured during fitting!");
1191    }
1192    Double_t *errors = new Double_t[entries];
1193    memcpy(errors,  tree->GetV1(), entries*sizeof(Double_t));
1194
1195    for (Int_t i = 0; i < dim + 1; i++){
1196       Int_t centries = 0;
1197       if (i < dim) centries = tree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start);
1198       else  centries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start);
1199
1200       if (entries != centries) {
1201         delete []errors;
1202         delete []values;
1203         delete formulaTokens;
1204         return new TString("An ERROR has occured during fitting!");
1205       }
1206       values[i] = new Double_t[entries];
1207       memcpy(values[i],  tree->GetV1(), entries*sizeof(Double_t));
1208    }
1209
1210    // add points to the fitter
1211    for (Int_t i = 0; i < entries; i++){
1212       Double_t x[1000];
1213       for (Int_t j=0; j<dim;j++) x[j]=values[j][i];
1214       fitter->AddPoint(x, values[dim][i], errors[i]);
1215    }
1216
1217    fitter->Eval();
1218    if (frac>0.5 && frac<1){
1219      fitter->EvalRobust(frac);
1220    }
1221    fitter->GetParameters(fitParam);
1222    fitter->GetCovarianceMatrix(covMatrix);
1223    chi2 = fitter->GetChisquare();
1224    npoints = entries;
1225
1226    TString *preturnFormula = new TString("("), &returnFormula = *preturnFormula;
1227
1228    for (Int_t iparam = 0; iparam < dim; iparam++) {
1229      returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam]));
1230      if (iparam < dim-1) returnFormula.Append("+");
1231    }
1232    returnFormula.Append(" )");
1233
1234
1235    for (Int_t j=0; j<dim+1;j++) delete [] values[j];
1236
1237    delete formulaTokens;
1238    delete fitter;
1239    delete[] values;
1240    delete[] errors;
1241    return preturnFormula;
1242 }
1243
1244
1245
1246
1247
1248 Int_t TStatToolkit::GetFitIndex(const TString fString, const TString subString){
1249   //
1250   // fitString - ++ separated list of fits
1251   // substring - ++ separated list of the requiered substrings
1252   //
1253   // return the last occurance of substring in fit string
1254   //
1255   TObjArray *arrFit = fString.Tokenize("++");
1256   TObjArray *arrSub = subString.Tokenize("++");
1257   Int_t index=-1;
1258   for (Int_t i=0; i<arrFit->GetEntries(); i++){
1259     Bool_t isOK=kTRUE;
1260     TString str =arrFit->At(i)->GetName();
1261     for (Int_t isub=0; isub<arrSub->GetEntries(); isub++){
1262       if (str.Contains(arrSub->At(isub)->GetName())==0) isOK=kFALSE;
1263     }
1264     if (isOK) index=i;
1265   }
1266   delete arrFit;
1267   delete arrSub;
1268   return index;
1269 }
1270
1271
1272 TString  TStatToolkit::FilterFit(const TString &input, const TString filter, TVectorD &param, TMatrixD & covar){
1273   //
1274   // Filter fit expression make sub-fit
1275   //
1276   TObjArray *array0= input.Tokenize("++");
1277   TObjArray *array1= filter.Tokenize("++");
1278   //TString *presult=new TString("(0");
1279   TString result="(0.0";
1280   for (Int_t i=0; i<array0->GetEntries(); i++){
1281     Bool_t isOK=kTRUE;
1282     TString str(array0->At(i)->GetName());
1283     for (Int_t j=0; j<array1->GetEntries(); j++){
1284       if (str.Contains(array1->At(j)->GetName())==0) isOK=kFALSE;
1285     }
1286     if (isOK) {
1287       result+="+"+str;
1288       result+=Form("*(%f)",param[i+1]);
1289       printf("%f\t%f\t%s\n",param[i+1], TMath::Sqrt(covar(i+1,i+1)),str.Data());
1290     }
1291   }
1292   result+="-0.)";
1293   delete array0;
1294   delete array1;
1295   return result;
1296 }
1297
1298 void TStatToolkit::Update1D(Double_t delta, Double_t sigma, Int_t s1, TMatrixD &vecXk, TMatrixD &covXk){
1299   //
1300   // Update parameters and covariance - with one measurement
1301   // Input:
1302   // vecXk - input vector - Updated in function
1303   // covXk - covariance matrix - Updated in function
1304   // delta, sigma, s1 - new measurement, rms of new measurement and the index of measurement
1305   const Int_t knMeas=1;
1306   Int_t knElem=vecXk.GetNrows();
1307
1308   TMatrixD mat1(knElem,knElem);            // update covariance matrix
1309   TMatrixD matHk(1,knElem);        // vector to mesurement
1310   TMatrixD vecYk(knMeas,1);        // Innovation or measurement residual
1311   TMatrixD matHkT(knElem,knMeas);  // helper matrix Hk transpose
1312   TMatrixD matSk(knMeas,knMeas);   // Innovation (or residual) covariance
1313   TMatrixD matKk(knElem,knMeas);   // Optimal Kalman gain
1314   TMatrixD covXk2(knElem,knElem);  // helper matrix
1315   TMatrixD covXk3(knElem,knElem);  // helper matrix
1316   TMatrixD vecZk(1,1);
1317   TMatrixD measR(1,1);
1318   vecZk(0,0)=delta;
1319   measR(0,0)=sigma*sigma;
1320   //
1321   // reset matHk
1322   for (Int_t iel=0;iel<knElem;iel++)
1323     for (Int_t ip=0;ip<knMeas;ip++) matHk(ip,iel)=0;
1324   //mat1
1325   for (Int_t iel=0;iel<knElem;iel++) {
1326     for (Int_t jel=0;jel<knElem;jel++) mat1(iel,jel)=0;
1327     mat1(iel,iel)=1;
1328   }
1329   //
1330   matHk(0, s1)=1;
1331   vecYk = vecZk-matHk*vecXk;               // Innovation or measurement residual
1332   matHkT=matHk.T(); matHk.T();
1333   matSk = (matHk*(covXk*matHkT))+measR;    // Innovation (or residual) covariance
1334   matSk.Invert();
1335   matKk = (covXk*matHkT)*matSk;            //  Optimal Kalman gain
1336   vecXk += matKk*vecYk;                    //  updated vector
1337   covXk2= (mat1-(matKk*matHk));
1338   covXk3 =  covXk2*covXk;
1339   covXk = covXk3;
1340   Int_t nrows=covXk3.GetNrows();
1341
1342   for (Int_t irow=0; irow<nrows; irow++)
1343     for (Int_t icol=0; icol<nrows; icol++){
1344       // rounding problems - make matrix again symteric
1345       covXk(irow,icol)=(covXk3(irow,icol)+covXk3(icol,irow))*0.5;
1346     }
1347 }
1348
1349
1350
1351 void   TStatToolkit::Constrain1D(const TString &input, const TString filter, TVectorD &param, TMatrixD & covar, Double_t mean, Double_t sigma){
1352   //
1353   // constrain linear fit
1354   // input  - string description of fit function
1355   // filter - string filter to select sub fits
1356   // param,covar - parameters and covariance matrix of the fit
1357   // mean,sigma  - new measurement uning which the fit is updated
1358   //
1359
1360   TObjArray *array0= input.Tokenize("++");
1361   TObjArray *array1= filter.Tokenize("++");
1362   TMatrixD paramM(param.GetNrows(),1);
1363   for (Int_t i=0; i<=array0->GetEntries(); i++){paramM(i,0)=param(i);}
1364
1365   if (filter.Length()==0){
1366     TStatToolkit::Update1D(mean, sigma, 0, paramM, covar);//
1367   }else{
1368     for (Int_t i=0; i<array0->GetEntries(); i++){
1369       Bool_t isOK=kTRUE;
1370       TString str(array0->At(i)->GetName());
1371       for (Int_t j=0; j<array1->GetEntries(); j++){
1372         if (str.Contains(array1->At(j)->GetName())==0) isOK=kFALSE;
1373       }
1374       if (isOK) {
1375         TStatToolkit::Update1D(mean, sigma, i+1, paramM, covar);//
1376       }
1377     }
1378   }
1379   for (Int_t i=0; i<=array0->GetEntries(); i++){
1380     param(i)=paramM(i,0);
1381   }
1382   delete array0;
1383   delete array1;
1384 }
1385
1386 TString  TStatToolkit::MakeFitString(const TString &input, const TVectorD &param, const TMatrixD & covar, Bool_t verbose){
1387   //
1388   //
1389   //
1390   TObjArray *array0= input.Tokenize("++");
1391   TString result=Form("(%f",param[0]);
1392   printf("%f\t%f\t\n", param[0], TMath::Sqrt(covar(0,0)));
1393   for (Int_t i=0; i<array0->GetEntries(); i++){
1394     TString str(array0->At(i)->GetName());
1395     result+="+"+str;
1396     result+=Form("*(%f)",param[i+1]);
1397     if (verbose) printf("%f\t%f\t%s\n", param[i+1], TMath::Sqrt(covar(i+1,i+1)),str.Data());
1398   }
1399   result+="-0.)";
1400   delete array0;
1401   return result;
1402 }
1403
1404 TGraphErrors * TStatToolkit::MakeGraphErrors(TTree * tree, const char * expr, const char * cut,  Int_t mstyle, Int_t mcolor, Float_t msize, Float_t offset){
1405   //
1406   // Query a graph errors
1407   // return TGraphErrors specified by expr and cut
1408   // Example  usage TStatToolkit::MakeGraphError(tree,"Y:X:ErrY","X>0", 25,2,0.4)
1409   // tree   - tree with variable
1410   // expr   - examp
1411   const Int_t entries =  tree->Draw(expr,cut,"goff");
1412   if (entries<=0) {
1413     TStatToolkit t;
1414     t.Error("TStatToolkit::MakeGraphError",Form("Empty or Not valid expression (%s) or cut *%s)", expr,cut));
1415     return 0;
1416   }
1417   if (  tree->GetV2()==0){
1418     TStatToolkit t;
1419     t.Error("TStatToolkit::MakeGraphError",Form("Not valid expression (%s) ", expr));
1420     return 0;
1421   }
1422   TGraphErrors * graph=0;
1423   if ( tree->GetV3()!=0){
1424     graph = new TGraphErrors (entries, tree->GetV2(),tree->GetV1(),0,tree->GetV3());
1425   }else{
1426     graph = new TGraphErrors (entries, tree->GetV2(),tree->GetV1(),0,0);
1427   }
1428   graph->SetMarkerStyle(mstyle);
1429   graph->SetMarkerColor(mcolor);
1430   graph->SetLineColor(mcolor);
1431   if (msize>0) graph->SetMarkerSize(msize);
1432   for(Int_t i=0;i<graph->GetN();i++) graph->GetX()[i]+=offset;
1433   return graph;
1434
1435 }
1436
1437
1438 TGraph * TStatToolkit::MakeGraphSparse(TTree * tree, const char * expr, const char * cut, Int_t mstyle, Int_t mcolor, Float_t msize, Float_t offset){
1439   //
1440   // Make a sparse draw of the variables
1441   // Format of expr : Var:Run or Var:Run:ErrorY or Var:Run:ErrorY:ErrorX
1442   // offset : points can slightly be shifted in x for better visibility with more graphs
1443   //
1444   // Written by Weilin.Yu
1445   // updated & merged with QA-code by Patrick Reichelt
1446   //
1447   const Int_t entries = tree->Draw(expr,cut,"goff");
1448   if (entries<=0) {
1449     TStatToolkit t;
1450     t.Error("TStatToolkit::MakeGraphSparse",Form("Empty or Not valid expression (%s) or cut (%s)", expr, cut));
1451     return 0;
1452   }
1453   //  TGraph * graph = (TGraph*)gPad->GetPrimitive("Graph"); // 2D
1454
1455   Double_t *graphY, *graphX;
1456   graphY = tree->GetV1();
1457   graphX = tree->GetV2();
1458
1459   // sort according to run number
1460   Int_t *index = new Int_t[entries*4];
1461   TMath::Sort(entries,graphX,index,kFALSE);
1462
1463   // define arrays for the new graph
1464   Double_t *unsortedX = new Double_t[entries];
1465   Int_t *runNumber = new Int_t[entries];
1466   Double_t count = 0.5;
1467
1468   // evaluate arrays for the new graph according to the run-number
1469   Int_t icount=0;
1470   //first entry
1471   unsortedX[index[0]] = count;
1472   runNumber[0] = graphX[index[0]];
1473   // loop the rest of entries
1474   for(Int_t i=1;i<entries;i++)
1475   {
1476     if(graphX[index[i]]==graphX[index[i-1]])
1477       unsortedX[index[i]] = count;
1478     else if(graphX[index[i]]!=graphX[index[i-1]]){
1479       count++;
1480       icount++;
1481       unsortedX[index[i]] = count;
1482       runNumber[icount]=graphX[index[i]];
1483     }
1484   }
1485
1486   // count the number of xbins (run-wise) for the new graph
1487   const Int_t newNbins = int(count+0.5);
1488   Double_t *newBins = new Double_t[newNbins+1];
1489   for(Int_t i=0; i<=count+1;i++){
1490     newBins[i] = i;
1491   }
1492
1493   // define and fill the new graph
1494   TGraph *graphNew = 0;
1495   if (tree->GetV3()) {
1496     if (tree->GetV4()) {
1497       graphNew = new TGraphErrors(entries,unsortedX,graphY,tree->GetV4(),tree->GetV3());
1498     }
1499     else { graphNew = new TGraphErrors(entries,unsortedX,graphY,0,tree->GetV3()); }
1500   }
1501   else { graphNew = new TGraphErrors(entries,unsortedX,graphY,0,0); }
1502   // with "Set(...)", the x-axis is being sorted
1503   graphNew->GetXaxis()->Set(newNbins,newBins);
1504
1505   // set the bins for the x-axis, apply shifting of points
1506   Char_t xName[50];
1507   for(Int_t i=0;i<count;i++){
1508     snprintf(xName,50,"%d",runNumber[i]);
1509     graphNew->GetXaxis()->SetBinLabel(i+1,xName);
1510     graphNew->GetX()[i]+=offset;
1511   }
1512
1513   graphNew->GetHistogram()->SetTitle("");
1514   graphNew->SetMarkerStyle(mstyle);
1515   graphNew->SetMarkerColor(mcolor);
1516   if (msize>0) graphNew->SetMarkerSize(msize);
1517   delete [] unsortedX;
1518   delete [] runNumber;
1519   delete [] index;
1520   delete [] newBins;
1521   //
1522   return graphNew;
1523 }
1524
1525
1526
1527 //
1528 // functions used for the trending
1529 //
1530
1531 Int_t  TStatToolkit::MakeStatAlias(TTree * tree, const char * expr, const char * cut, const char * alias)
1532 {
1533   //
1534   // Add alias using statistical values of a given variable.
1535   // (by MI, Patrick Reichelt)
1536   //
1537   // tree - input tree
1538   // expr - variable expression
1539   // cut  - selection criteria
1540   // Output - return number of entries used to define variable
1541   // In addition mean, rms, median, and robust mean and rms (choosing fraction of data with smallest RMS)
1542   //
1543   /* Example usage:
1544      1.) create the robust estimators for variable expr="QA.TPC.CPass1.meanTPCncl" and create a corresponding
1545      aliases with the prefix alias[0]="ncl", calculated using fraction alias[1]="0.90"
1546
1547      TStatToolkit::MakeStatAlias(tree,"QA.TPC.CPass1.meanTPCncl","QA.TPC.CPass1.status>0","ncl:0.9");
1548      root [4] tree->GetListOfAliases().Print()
1549      OBJ: TNamed    ncl_Median      (130.964333+0)
1550      OBJ: TNamed    ncl_Mean        (122.120387+0)
1551      OBJ: TNamed    ncl_RMS         (33.509623+0)
1552      OBJ: TNamed    ncl_Mean90      (131.503862+0)
1553      OBJ: TNamed    ncl_RMS90       (3.738260+0)
1554   */
1555   //
1556   Int_t entries = tree->Draw(expr,cut,"goff");
1557   if (entries<=1){
1558     printf("Expression or cut not valid:\t%s\t%s\n", expr, cut);
1559     return 0;
1560   }
1561   //
1562   TObjArray* oaAlias = TString(alias).Tokenize(":");
1563   if (oaAlias->GetEntries()<2) return 0;
1564   Float_t entryFraction = atof( oaAlias->At(1)->GetName() );
1565   //
1566   Double_t median = TMath::Median(entries,tree->GetV1());
1567   Double_t mean   = TMath::Mean(entries,tree->GetV1());
1568   Double_t rms    = TMath::RMS(entries,tree->GetV1());
1569   Double_t meanEF=0, rmsEF=0;
1570   TStatToolkit::EvaluateUni(entries, tree->GetV1(), meanEF, rmsEF, entries*entryFraction);
1571   //
1572   tree->SetAlias(Form("%s_Median",oaAlias->At(0)->GetName()), Form("(%f+0)",median));
1573   tree->SetAlias(Form("%s_Mean",oaAlias->At(0)->GetName()), Form("(%f+0)",mean));
1574   tree->SetAlias(Form("%s_RMS",oaAlias->At(0)->GetName()), Form("(%f+0)",rms));
1575   tree->SetAlias(Form("%s_Mean%d",oaAlias->At(0)->GetName(),Int_t(entryFraction*100)), Form("(%f+0)",meanEF));
1576   tree->SetAlias(Form("%s_RMS%d",oaAlias->At(0)->GetName(),Int_t(entryFraction*100)), Form("(%f+0)",rmsEF));
1577   delete oaAlias;
1578   return entries;
1579 }
1580
1581 Int_t  TStatToolkit::SetStatusAlias(TTree * tree, const char * expr, const char * cut, const char * alias)
1582 {
1583   //
1584   // Add alias to trending tree using statistical values of a given variable.
1585   // (by MI, Patrick Reichelt)
1586   //
1587   // format of expr :  varname (e.g. meanTPCncl)
1588   // format of cut  :  char like in TCut
1589   // format of alias:  alias:query:entryFraction(EF) (fraction of entries used for uniformity evaluation)
1590   //            e.g.:  varname_Out:(abs(varname-meanEF)>6.*rmsEF):0.8
1591   // available internal variables are: 'varname, Median, Mean, MeanEF, RMS, RMSEF'
1592   // in the alias, 'varname' will be replaced by its content, and 'EF' by the percentage (e.g. MeanEF -> Mean80)
1593   //
1594   /* Example usage:
1595      1.) Define robust mean (possible, but easier done with TStatToolkit::MakeStatAlias(...))
1596      TStatToolkit::SetStatusAlias(tree, "meanTPCnclF", "meanTPCnclF>0", "meanTPCnclF_MeanEF:MeanEF:0.80") ;
1597      root [10] tree->GetListOfAliases()->Print()
1598                Collection name='TList', class='TList', size=1
1599                OBJ: TNamed    meanTPCnclF_Mean80      0.899308
1600      2.) create alias outlyers  - 6 sigma cut
1601      TStatToolkit::SetStatusAlias(tree, "meanTPCnclF", "meanTPCnclF>0", "meanTPCnclF_Out:(abs(meanTPCnclF-MeanEF)>6.*RMSEF):0.8")
1602      meanTPCnclF_Out ==> (abs(meanTPCnclF-0.899308)>6.*0.016590)
1603      3.) the same functionality as in 2.)
1604      TStatToolkit::SetStatusAlias(tree, "meanTPCnclF", "meanTPCnclF>0", "varname_Out2:(abs(varname-MeanEF)>6.*RMSEF):0.8")
1605      meanTPCnclF_Out2 ==> (abs(meanTPCnclF-0.899308)>6.*0.016590)
1606   */
1607   //
1608   Int_t entries = tree->Draw(expr,cut,"goff");
1609   if (entries<=1){
1610     printf("Expression or cut not valid:\t%s\t%s\n", expr, cut);
1611     return 0;
1612   }
1613   //
1614   TObjArray* oaVar = TString(expr).Tokenize(":");
1615   char varname[50];
1616   snprintf(varname,50,"%s", oaVar->At(0)->GetName());
1617   //
1618   TObjArray* oaAlias = TString(alias).Tokenize(":");
1619   if (oaAlias->GetEntries()<3) return 0;
1620   Float_t entryFraction = atof( oaAlias->At(2)->GetName() );
1621   //
1622   Double_t median = TMath::Median(entries,tree->GetV1());
1623   Double_t mean   = TMath::Mean(entries,tree->GetV1());
1624   Double_t rms    = TMath::RMS(entries,tree->GetV1());
1625   Double_t meanEF=0, rmsEF=0;
1626   TStatToolkit::EvaluateUni(entries, tree->GetV1(), meanEF, rmsEF, entries*entryFraction);
1627   //
1628   TString sAlias( oaAlias->At(0)->GetName() );
1629   sAlias.ReplaceAll("varname",varname);
1630   sAlias.ReplaceAll("MeanEF", Form("Mean%1.0f",entryFraction*100) );
1631   sAlias.ReplaceAll("RMSEF",  Form("RMS%1.0f",entryFraction*100) );
1632   TString sQuery( oaAlias->At(1)->GetName() );
1633   sQuery.ReplaceAll("varname",varname);
1634   sQuery.ReplaceAll("MeanEF", Form("%f",meanEF) );
1635   sQuery.ReplaceAll("RMSEF",  Form("%f",rmsEF) ); //make sure to replace 'RMSEF' before 'RMS'...
1636   sQuery.ReplaceAll("Median", Form("%f",median) );
1637   sQuery.ReplaceAll("Mean",   Form("%f",mean) );
1638   sQuery.ReplaceAll("RMS",    Form("%f",rms) );
1639   printf("define alias:\t%s = %s\n", sAlias.Data(), sQuery.Data());
1640   //
1641   char query[200];
1642   char aname[200];
1643   snprintf(query,200,"%s", sQuery.Data());
1644   snprintf(aname,200,"%s", sAlias.Data());
1645   tree->SetAlias(aname, query);
1646   delete oaVar;
1647   delete oaAlias;
1648   return entries;
1649 }
1650
1651 TMultiGraph*  TStatToolkit::MakeStatusMultGr(TTree * tree, const char * expr, const char * cut, const char * alias, Int_t igr)
1652 {
1653   //
1654   // Compute a trending multigraph that shows for which runs a variable has outliers.
1655   // (by MI, Patrick Reichelt)
1656   //
1657   // format of expr :  varname:xaxis (e.g. meanTPCncl:run)
1658   // format of cut  :  char like in TCut
1659   // format of alias:  (1):(varname_Out==0):(varname_Out)[:(varname_Warning):...]
1660   // in the alias, 'varname' will be replaced by its content (e.g. varname_Out -> meanTPCncl_Out)
1661   // note: the aliases 'varname_Out' etc have to be defined by function TStatToolkit::SetStatusAlias(...)
1662   // counter igr is used to shift the multigraph in y when filling a TObjArray.
1663   //
1664   TObjArray* oaVar = TString(expr).Tokenize(":");
1665   if (oaVar->GetEntries()<2) return 0;
1666   char varname[50];
1667   char var_x[50];
1668   snprintf(varname,50,"%s", oaVar->At(0)->GetName());
1669   snprintf(var_x  ,50,"%s", oaVar->At(1)->GetName());
1670   //
1671   TString sAlias(alias);
1672   sAlias.ReplaceAll("varname",varname);
1673   TObjArray* oaAlias = TString(sAlias.Data()).Tokenize(":");
1674   if (oaAlias->GetEntries()<3) return 0;
1675   //
1676   char query[200];
1677   TMultiGraph* multGr = new TMultiGraph();
1678   Int_t marArr[6]    = {24+igr%2, 20+igr%2, 20+igr%2, 20+igr%2, 22, 23};
1679   Int_t colArr[6]    = {kBlack, kBlack, kRed, kOrange, kMagenta, kViolet};
1680   Double_t sizArr[6] = {1.2, 1.1, 1.0, 1.0, 1, 1};
1681   const Int_t ngr = oaAlias->GetEntriesFast();
1682   for (Int_t i=0; i<ngr; i++){
1683     if (i==2) continue; // the Fatal(Out) graph will be added in the end to be plotted on top!
1684     snprintf(query,200, "%f*(%s-0.5):%s", 1.+igr, oaAlias->At(i)->GetName(), var_x);
1685     multGr->Add( (TGraphErrors*) TStatToolkit::MakeGraphSparse(tree,query,cut,marArr[i],colArr[i],sizArr[i]) );
1686   }
1687   snprintf(query,200, "%f*(%s-0.5):%s", 1.+igr, oaAlias->At(2)->GetName(), var_x);
1688   multGr->Add( (TGraphErrors*) TStatToolkit::MakeGraphSparse(tree,query,cut,marArr[2],colArr[2],sizArr[2]) );
1689   //
1690   multGr->SetName(varname);
1691   multGr->SetTitle(varname); // used for y-axis labels. // details to be included!
1692   delete oaVar;
1693   delete oaAlias;
1694   return multGr;
1695 }
1696
1697
1698 void  TStatToolkit::AddStatusPad(TCanvas* c1, Float_t padratio, Float_t bottommargin)
1699 {
1700   //
1701   // add pad to bottom of canvas for Status graphs (by Patrick Reichelt)
1702   // call function "DrawStatusGraphs(...)" afterwards
1703   //
1704   TCanvas* c1_clone = (TCanvas*) c1->Clone("c1_clone");
1705   c1->Clear();
1706   // produce new pads
1707   c1->cd();
1708   TPad* pad1 = new TPad("pad1", "pad1", 0., padratio, 1., 1.);
1709   pad1->Draw();
1710   pad1->SetNumber(1); // so it can be called via "c1->cd(1);"
1711   c1->cd();
1712   TPad* pad2 = new TPad("pad2", "pad2", 0., 0., 1., padratio);
1713   pad2->Draw();
1714   pad2->SetNumber(2);
1715   // draw original canvas into first pad
1716   c1->cd(1);
1717   c1_clone->DrawClonePad();
1718   pad1->SetBottomMargin(0.001);
1719   pad1->SetRightMargin(0.01);
1720   // set up second pad
1721   c1->cd(2);
1722   pad2->SetGrid(3);
1723   pad2->SetTopMargin(0);
1724   pad2->SetBottomMargin(bottommargin); // for the long x-axis labels (runnumbers)
1725   pad2->SetRightMargin(0.01);
1726 }
1727
1728
1729 void  TStatToolkit::DrawStatusGraphs(TObjArray* oaMultGr)
1730 {
1731   //
1732   // draw Status graphs into active pad of canvas (by MI, Patrick Reichelt)
1733   // ...into bottom pad, if called after "AddStatusPad(...)"
1734   //
1735   const Int_t nvars = oaMultGr->GetEntriesFast();
1736   TGraph* grAxis = (TGraph*) ((TMultiGraph*) oaMultGr->At(0))->GetListOfGraphs()->At(0);
1737   grAxis->SetMaximum(0.5*nvars+0.5);
1738   grAxis->SetMinimum(0);
1739   grAxis->GetYaxis()->SetLabelSize(0);
1740   Int_t entries = grAxis->GetN();
1741   printf("entries (via GetN()) = %d\n",entries);
1742   grAxis->GetXaxis()->SetLabelSize(5.7*TMath::Min(TMath::Max(5./entries,0.01),0.03));
1743   grAxis->GetXaxis()->LabelsOption("v");
1744   grAxis->Draw("ap");
1745   //
1746   // draw multigraphs & names of status variables on the y axis
1747   for (Int_t i=0; i<nvars; i++){
1748     ((TMultiGraph*) oaMultGr->At(i))->Draw("p");
1749     TLatex* ylabel = new TLatex(-0.1, 0.5*i+0.5, ((TMultiGraph*) oaMultGr->At(i))->GetTitle());
1750     ylabel->SetTextAlign(32); //hor:right & vert:centered
1751     ylabel->SetTextSize(0.025/gPad->GetHNDC());
1752     ylabel->Draw();
1753   }
1754 }
1755
1756
1757 void TStatToolkit::DrawHistogram(TTree * tree, const char* drawCommand, const char* cuts, const char* histoname, const char* histotitle, Int_t nsigma, Float_t fraction )
1758 {
1759   //
1760   // Draw histogram from TTree with robust range
1761   // Only for 1D so far!
1762   //
1763   // Parameters:
1764   // - histoname:  name of histogram
1765   // - histotitle: title of histgram
1766   // - fraction:   fraction of data to define the robust mean
1767   // - nsigma:     nsigma value for range
1768   //
1769
1770    TString drawStr(drawCommand);
1771    TString cutStr(cuts);
1772    Int_t dim = 1;
1773
1774    if(!tree) {
1775      cerr<<" Tree pointer is NULL!"<<endl;
1776      return;
1777    }
1778
1779    // get entries
1780    Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff");
1781    if (entries == -1) {
1782      cerr<<"TTree draw returns -1"<<endl;
1783      return;
1784    }
1785
1786    // get dimension
1787    if(tree->GetV1()) dim = 1;
1788    if(tree->GetV2()) dim = 2;
1789    if(tree->GetV3()) dim = 3;
1790    if(dim > 2){
1791      cerr<<"TTree has more than 2 dimensions (not yet supported)"<<endl;
1792      return;
1793    }
1794
1795    // draw robust
1796    Double_t meanX, rmsX=0;
1797    Double_t meanY, rmsY=0;
1798    TStatToolkit::EvaluateUni(entries, tree->GetV1(),meanX,rmsX, fraction*entries);
1799    if(dim==2){
1800      TStatToolkit::EvaluateUni(entries, tree->GetV1(),meanY,rmsY, fraction*entries);
1801      TStatToolkit::EvaluateUni(entries, tree->GetV2(),meanX,rmsX, fraction*entries);
1802    }
1803    TH1* hOut;
1804    if(dim==1){
1805      hOut = new TH1F(histoname, histotitle, 200, meanX-nsigma*rmsX, meanX+nsigma*rmsX);
1806      for (Int_t i=0; i<entries; i++) hOut->Fill(tree->GetV1()[i]);
1807      hOut->GetXaxis()->SetTitle(tree->GetHistogram()->GetXaxis()->GetTitle());
1808      hOut->Draw();
1809    }
1810    else if(dim==2){
1811      hOut = new TH2F(histoname, histotitle, 200, meanX-nsigma*rmsX, meanX+nsigma*rmsX,200, meanY-nsigma*rmsY, meanY+nsigma*rmsY);
1812      for (Int_t i=0; i<entries; i++) hOut->Fill(tree->GetV2()[i],tree->GetV1()[i]);
1813      hOut->GetXaxis()->SetTitle(tree->GetHistogram()->GetXaxis()->GetTitle());
1814      hOut->GetYaxis()->SetTitle(tree->GetHistogram()->GetYaxis()->GetTitle());
1815      hOut->Draw("colz");
1816    }
1817
1818 }