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