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