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
1218 TGraph * TStatToolkit::MakeGraphSparse(TTree * tree, const char * expr, const char * cut, Int_t mstyle, Int_t mcolor, Float_t msize, Float_t offset){
1219   //
1220   // Make a sparse draw of the variables
1221   // Writen by Weilin.Yu
1222   const Int_t entries =  tree->Draw(expr,cut,"goff");
1223   if (entries<=0) {
1224     TStatToolkit t;
1225     t.Error("TStatToolkit::MakeGraphSparse",Form("Empty or Not valid expression (%s) or cut *%s)", expr,cut));
1226     return 0;
1227   }
1228   //  TGraph * graph = (TGraph*)gPad->GetPrimitive("Graph"); // 2D
1229   TGraph * graph = 0;
1230   if (tree->GetV3()) graph = new TGraphErrors (entries, tree->GetV2(),tree->GetV1(),0,tree->GetV3());
1231   graph =  new TGraphErrors (entries, tree->GetV2(),tree->GetV1(),0,0);
1232   graph->SetMarkerStyle(mstyle);
1233   graph->SetMarkerColor(mcolor);
1234   //
1235   Int_t *index = new Int_t[entries*4];
1236   TMath::Sort(entries,graph->GetX(),index,kFALSE);
1237
1238   Double_t *tempArray = new Double_t[entries];
1239
1240   Double_t count = 0.5;
1241   Double_t  *vrun = new Double_t[entries];
1242   Int_t icount=0;
1243   //
1244   tempArray[index[0]] = count;
1245   vrun[0] = graph->GetX()[index[0]];
1246   for(Int_t i=1;i<entries;i++){
1247     if(graph->GetX()[index[i]]==graph->GetX()[index[i-1]])
1248       tempArray[index[i]] = count;
1249     else if(graph->GetX()[index[i]]!=graph->GetX()[index[i-1]]){
1250       count++;
1251       icount++;
1252       tempArray[index[i]] = count;
1253       vrun[icount]=graph->GetX()[index[i]];
1254     }
1255   }
1256
1257   const Int_t newNbins = int(count+0.5);
1258   Double_t *newBins = new Double_t[newNbins+1];
1259   for(Int_t i=0; i<=count+1;i++){
1260     newBins[i] = i;
1261   }
1262
1263   TGraph *graphNew = 0;
1264   if (tree->GetV3()) graphNew = new TGraphErrors(entries,tempArray,graph->GetY(),0,tree->GetV3());
1265   else
1266     graphNew = new TGraphErrors(entries,tempArray,graph->GetY(),0,0);
1267   graphNew->GetXaxis()->Set(newNbins,newBins);
1268
1269   Char_t xName[50];
1270   for(Int_t i=0;i<count;i++){
1271     snprintf(xName,50,"%d",Int_t(vrun[i]));
1272     graphNew->GetXaxis()->SetBinLabel(i+1,xName);
1273   }
1274   graphNew->GetHistogram()->SetTitle("");
1275   graphNew->SetMarkerStyle(mstyle);
1276   graphNew->SetMarkerColor(mcolor);
1277   if (msize>0) graphNew->SetMarkerSize(msize);
1278   for(Int_t i=0;i<graphNew->GetN();i++) graphNew->GetX()[i]+=offset;
1279   delete [] tempArray;
1280   delete [] index;
1281   delete [] newBins;
1282   delete [] vrun;
1283   return graphNew;
1284 }
1285
1286
1287
1288 //
1289 // function used for the trending
1290 //
1291
1292 Int_t  TStatToolkit::MakeStatAlias(TTree * tree, const char * expr, const char * cut, const char * alias)
1293 {
1294   //
1295   // Add alias using statistical values of a given variable.
1296   // (by MI, Patrick Reichelt)
1297   //
1298   // tree - input tree
1299   // expr - variable expression
1300   // cut  - selection criteria
1301   // Output - return number of entries used to define variable
1302   // In addition mean, rms, median, and robust mean and rms (choosing fraction of data with smallest RMS)
1303   //
1304   // Example usage:
1305   /*
1306     Example usage to create the robust estimators for  variable expr="QA.TPC.CPass1.meanTPCncl" and create a corresponding
1307     aliases with the prefix alias[0]="ncl", calculated using fraction alias[1]="0.90"
1308
1309     TStatToolkit::MakeStatAlias(tree,"QA.TPC.CPass1.meanTPCncl","QA.TPC.CPass1.status>0","ncl:0.9");
1310     root [4] tree->GetListOfAliases().Print()
1311     OBJ: TNamed    ncl_Mean        (122.120387+0)
1312     OBJ: TNamed    ncl_RMS (33.509623+0)
1313     OBJ: TNamed    ncl_Median      (130.964333+0)
1314     OBJ: TNamed    ncl_Mean90      (131.503862+0)
1315     OBJ: TNamed    ncl_RMS90       (3.738260+0)
1316    */
1317   //
1318   Int_t entries= tree->Draw(expr,cut,"goff");
1319   if (entries<=1){
1320     printf("Expression or cut not valid:\t%s\t%s\n", expr, cut);
1321     return 0;
1322   }
1323   //
1324   TObjArray* oaAlias = TString(alias).Tokenize(":");
1325   if (oaAlias->GetEntries()<2) return 0;
1326   Float_t entryFraction = atof( oaAlias->At(1)->GetName() );
1327   //
1328   Double_t median = TMath::Median(entries,tree->GetV1());
1329   Double_t mean = TMath::Mean(entries,tree->GetV1());
1330   Double_t rms    = TMath::RMS(entries,tree->GetV1());
1331   Double_t meanEF=0, rmsEF=0;
1332   TStatToolkit::EvaluateUni(entries, tree->GetV1(), meanEF, rmsEF, entries*entryFraction);
1333   //
1334   tree->SetAlias(Form("%s_Mean",oaAlias->At(0)->GetName()), Form("(%f+0)",mean));
1335   tree->SetAlias(Form("%s_RMS",oaAlias->At(0)->GetName()), Form("(%f+0)",rms));
1336   tree->SetAlias(Form("%s_Median",oaAlias->At(0)->GetName()), Form("(%f+0)",median));
1337   tree->SetAlias(Form("%s_Mean%d",oaAlias->At(0)->GetName(),Int_t(entryFraction*100)), Form("(%f+0)",meanEF));
1338   tree->SetAlias(Form("%s_RMS%d",oaAlias->At(0)->GetName(),Int_t(entryFraction*100)), Form("(%f+0)",rmsEF));
1339   delete oaAlias;
1340   return entries;
1341 }
1342
1343 Int_t  TStatToolkit::SetStatusAlias(TTree * tree, const char * expr, const char * cut, const char * alias)
1344 {
1345   //
1346   // Add alias to trending tree using statistical values of a given variable.
1347   // (by MI, Patrick Reichelt)
1348   //
1349   // format of expr :  varname (e.g. meanTPCncl)
1350   // format of cut  :  char like in TCut
1351   // format of alias:  alias:query:entryFraction(EF) (fraction of entries used for uniformity evaluation)
1352   //            e.g.:  varname_Out:(abs(varname-meanEF)>6.*rmsEF):0.8
1353   // available internal variables are: 'varname, median, meanEF, rms, rmsEF'
1354   // in the alias, 'varname' will be replaced by its content, and 'EF' by the percentage (e.g. meanEF -> mean80)
1355   //
1356   /* Example usage:
1357      1.) Define robust  mean
1358
1359      TStatToolkit::SetStatusAlias(tree, "meanTPCnclF", "meanTPCnclF>0", "meanTPCnclF_meanEF:meanEF:0.80") ;
1360      -->
1361      root [10] tree->GetListOfAliases()->Print()
1362                Collection name='TList', class='TList', size=1
1363                OBJ: TNamed    meanTPCnclF_mean80      0.899308
1364      2.) create alias outlyers  - 6 sigma cut
1365       TStatToolkit::SetStatusAlias(tree, "meanTPCnclF", "meanTPCnclF>0", "meanTPCnclF_Out:(abs(meanTPCnclF-meanEF)>6.*rmsEF):0.8")          meanTPCnclF_Out ==> (abs(meanTPCnclF-0.899308)>6.*0.016590)
1366      3.) the same functionality as in 2.)
1367          TStatToolkit::SetStatusAlias(tree, "meanTPCnclF", "meanTPCnclF>0", "varname_Out2:(abs(varname-meanEF)>6.*rmsEF):0.8")
1368          ->
1369          meanTPCnclF_Out2 ==> (abs(meanTPCnclF-0.899308)>6.*0.016590)
1370   */
1371   //
1372   TObjArray* oaVar = TString(expr).Tokenize(":");
1373   char varname[50];
1374   //char var_x[50];
1375   snprintf(varname,50,"%s", oaVar->At(0)->GetName());
1376   //snprintf(var_x  ,50,"%s", oaVar->At(1)->GetName());
1377   TCut userCut(cut);
1378   TObjArray* oaAlias = TString(alias).Tokenize(":");
1379   Float_t entryFraction = atof( oaAlias->At(2)->GetName() );
1380   //
1381   Int_t entries = tree->Draw(expr, userCut, "goff");
1382   if (entries<=1){
1383     printf("Expression or cut not valid:\t%s\t%s\n", expr, cut);
1384     return 0;
1385   }
1386   //printf("  entries (via tree->Draw(...)) = %d\n",entries);
1387   Double_t mean = TMath::Mean(entries,tree->GetV1());
1388   Double_t median = TMath::Median(entries,tree->GetV1());
1389   Double_t rms    = TMath::RMS(entries,tree->GetV1());
1390   Double_t meanEF=0, rmsEF=0;
1391   TStatToolkit::EvaluateUni(entries, tree->GetV1(), meanEF, rmsEF, entries*entryFraction);
1392   //printf("%s\t%f\t%f\t%f\t%f\n",varname, median, meanEF, rms, rmsEF);
1393   //
1394   TString sAlias( oaAlias->At(0)->GetName() );
1395   sAlias.ReplaceAll("varname",varname);
1396   sAlias.ReplaceAll("MeanEF", Form("mean%1.0f",entryFraction*100) );
1397   sAlias.ReplaceAll("RMSEF",  Form("rms%1.0f",entryFraction*100) );
1398   TString sQuery( oaAlias->At(1)->GetName() );
1399   sQuery.ReplaceAll("varname",varname);
1400   sQuery.ReplaceAll("MeanEF", Form("%f",meanEF) );
1401   sQuery.ReplaceAll("RMSEF",  Form("%f",rmsEF) ); //make sure to replace 'rmsEF' before 'rms'...
1402   sQuery.ReplaceAll("Median", Form("%f",median) );
1403   sQuery.ReplaceAll("RMS",    Form("%f",rms) );
1404   sQuery.ReplaceAll("Mean",    Form("%f",mean) );
1405   printf("%s\n", sQuery.Data());
1406   //
1407   char query[200];
1408   char aname[200];
1409   snprintf(query,200,"%s", sQuery.Data());
1410   snprintf(aname,200,"%s", sAlias.Data());
1411   tree->SetAlias(aname, query);
1412   return entries;
1413 }
1414
1415 TMultiGraph*  TStatToolkit::MakeStatusMultGr(TTree * tree, const char * expr, const char * cut, const char * alias, Int_t igr)
1416 {
1417   //
1418   // Compute a trending multigraph that shows for which runs a variable has outliers.
1419   // (by MI, Patrick Reichelt)
1420   //
1421   // format of expr :  varname:xaxis (e.g. meanTPCncl:run)
1422   // format of cut  :  char like in TCut
1423   // format of alias:  (1):(varname_Out==0):(varname_Out)[:(varname_Warning):...]
1424   // in the alias, 'varname' will be replaced by its content (e.g. varname_Out -> meanTPCncl_Out)
1425   // note: the aliases 'varname_Out' etc have to be defined by function 'SetStatisticAlias(...)'
1426   // counter igr is used to shift the multigraph in y when filling a TObjArray.
1427   //
1428   TObjArray* oaVar = TString(expr).Tokenize(":");
1429   char varname[50];
1430   char var_x[50];
1431   snprintf(varname,50,"%s", oaVar->At(0)->GetName());
1432   snprintf(var_x  ,50,"%s", oaVar->At(1)->GetName());
1433   TCut userCut(cut);
1434   TString sAlias(alias);
1435   sAlias.ReplaceAll("varname",varname);
1436   TObjArray* oaAlias = TString(sAlias.Data()).Tokenize(":");
1437   //
1438   char query[200];
1439   TMultiGraph* multGr = new TMultiGraph();
1440   Int_t marArr[6]    = {24+igr%2, 20+igr%2, 20+igr%2, 20+igr%2, 22, 23};
1441   Int_t colArr[6]    = {kBlack, kBlack, kRed, kOrange, kMagenta, kViolet};
1442   Double_t sizArr[6] = {1.2, 1.1, 1.0, 1.0, 1, 1};
1443   const Int_t ngr = oaAlias->GetEntriesFast();
1444   for (Int_t i=0; i<ngr; i++){
1445     if (i==2) continue; // the Fatal(Out) graph will be added in the end to be plotted on top!
1446     snprintf(query,200, "%f*(%s-0.5):%s", 1.+igr, oaAlias->At(i)->GetName(), var_x);
1447     multGr->Add( (TGraphErrors*) TStatToolkit::MakeGraphSparse(tree,query,userCut,marArr[i],colArr[i],sizArr[i]) );
1448   }
1449   snprintf(query,200, "%f*(%s-0.5):%s", 1.+igr, oaAlias->At(2)->GetName(), var_x);
1450   multGr->Add( (TGraphErrors*) TStatToolkit::MakeGraphSparse(tree,query,userCut,marArr[2],colArr[2],sizArr[2]) );
1451   //
1452   multGr->SetName(varname);
1453   multGr->SetTitle(varname); // used for y-axis labels. // details to be included!
1454   return multGr;
1455 }
1456
1457
1458 void  TStatToolkit::AddStatusPad(TCanvas* c1, Float_t padratio, Float_t bottommargin)
1459 {
1460   //
1461   // add pad to bottom of canvas for Status graphs (by Patrick Reichelt)
1462   // call function "DrawStatusGraphs(...)" afterwards
1463   //
1464   TCanvas* c1_clone = (TCanvas*) c1->Clone("c1_clone");
1465   c1->Clear();
1466   // produce new pads
1467   c1->cd();
1468   TPad* pad1 = new TPad("pad1", "pad1", 0., padratio, 1., 1.);
1469   pad1->Draw();
1470   pad1->SetNumber(1); // so it can be called via "c1->cd(1);"
1471   c1->cd();
1472   TPad* pad2 = new TPad("pad2", "pad2", 0., 0., 1., padratio);
1473   pad2->Draw();
1474   pad2->SetNumber(2);
1475   // draw original canvas into first pad
1476   c1->cd(1);
1477   c1_clone->DrawClonePad();
1478   pad1->SetBottomMargin(0.001);
1479   pad1->SetRightMargin(0.01);
1480   // set up second pad
1481   c1->cd(2);
1482   pad2->SetGrid(3);
1483   pad2->SetTopMargin(0);
1484   pad2->SetBottomMargin(bottommargin); // for the long x-axis labels (runnumbers)
1485   pad2->SetRightMargin(0.01);
1486 }
1487
1488
1489 void  TStatToolkit::DrawStatusGraphs(TObjArray* oaMultGr)
1490 {
1491   //
1492   // draw Status graphs into active pad of canvas (by MI, Patrick Reichelt)
1493   // ...into bottom pad, if called after "AddStatusPad(...)"
1494   //
1495   const Int_t nvars = oaMultGr->GetEntriesFast();
1496   TGraph* grAxis = (TGraph*) ((TMultiGraph*) oaMultGr->At(0))->GetListOfGraphs()->At(0);
1497   grAxis->SetMaximum(0.5*nvars+0.5);
1498   grAxis->SetMinimum(0);
1499   grAxis->GetYaxis()->SetLabelSize(0);
1500   Int_t entries = grAxis->GetN();
1501   printf("entries (via GetN()) = %d\n",entries);
1502   grAxis->GetXaxis()->SetLabelSize(5.7*TMath::Min(TMath::Max(5./entries,0.01),0.03));
1503   grAxis->GetXaxis()->LabelsOption("v");
1504   grAxis->Draw("ap");
1505   //
1506   // draw multigraphs & names of status variables on the y axis
1507   for (Int_t i=0; i<nvars; i++){
1508     ((TMultiGraph*) oaMultGr->At(i))->Draw("p");
1509     TLatex* ylabel = new TLatex(-0.1, 0.5*i+0.5, ((TMultiGraph*) oaMultGr->At(i))->GetTitle());
1510     ylabel->SetTextAlign(32); //hor:right & vert:centered
1511     ylabel->SetTextSize(0.025/gPad->GetHNDC());
1512     ylabel->Draw();
1513   }
1514 }