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
  37 //
  38 // includes neccessary for test functions
  39 //
  40 #include "TSystem.h"
  41 #include "TRandom.h"
  42 #include "TStopwatch.h"
  43 #include "TTreeStream.h"
  44
  45 #include "TStatToolkit.h"
  46
  47
  48 ClassImp(TStatToolkit) // Class implementation to enable ROOT I/O
  49
  50 TStatToolkit::TStatToolkit() : TObject()
  51 {
  52   //
  53   // Default constructor
  54   //
  55 }
  56 ///////////////////////////////////////////////////////////////////////////
  57 TStatToolkit::~TStatToolkit()
  58 {
  59   //
  60   // Destructor
  61   //
  62 }
  63
  64
  65 //_____________________________________________________________________________
  66 void TStatToolkit::EvaluateUni(Int_t nvectors, Double_t *data, Double_t &mean
  67                            , Double_t &sigma, Int_t hh)
  68 {
  69   //
  70   // Robust estimator in 1D case MI version - (faster than ROOT version)
  71   //
  72   // For the univariate case
  73   // estimates of location and scatter are returned in mean and sigma parameters
  74   // the algorithm works on the same principle as in multivariate case -
  75   // it finds a subset of size hh with smallest sigma, and then returns mean and
  76   // sigma of this subset
  77   //
  78
  79   if (hh==0)
  80     hh=(nvectors+2)/2;
  81   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};
  82   Int_t *index=new Int_t[nvectors];
  83   TMath::Sort(nvectors, data, index, kFALSE);
  84
  85   Int_t    nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11);
  86   Double_t factor = faclts[TMath::Max(0,nquant-1)];
  87
  88   Double_t sumx  =0;
  89   Double_t sumx2 =0;
  90   Int_t    bestindex = -1;
  91   Double_t bestmean  = 0;
  92   Double_t bestsigma = (data[index[nvectors-1]]-data[index[0]]+1.);   // maximal possible sigma
  93   bestsigma *=bestsigma;
  94
  95   for (Int_t i=0; i<hh; i++){
  96     sumx  += data[index[i]];
  97     sumx2 += data[index[i]]*data[index[i]];
  98   }
  99
 100   Double_t norm = 1./Double_t(hh);
 101   Double_t norm2 = (hh-1)>0 ? 1./Double_t(hh-1):1;
 102   for (Int_t i=hh; i<nvectors; i++){
 103     Double_t cmean  = sumx*norm;
 104     Double_t csigma = (sumx2 - hh*cmean*cmean)*norm2;
 105     if (csigma<bestsigma){
 106       bestmean  = cmean;
 107       bestsigma = csigma;
 108       bestindex = i-hh;
 109     }
 110
 111     sumx  += data[index[i]]-data[index[i-hh]];
 112     sumx2 += data[index[i]]*data[index[i]]-data[index[i-hh]]*data[index[i-hh]];
 113   }
 114
 115   Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma));
 116   mean  = bestmean;
 117   sigma = bstd;
 118   delete [] index;
 119
 120 }
 121
 122
 123
 124 void TStatToolkit::EvaluateUniExternal(Int_t nvectors, Double_t *data, Double_t &mean, Double_t &sigma, Int_t hh,  Float_t externalfactor)
 125 {
 126   // Modified version of ROOT robust EvaluateUni
 127   // robust estimator in 1D case MI version
 128   // added external factor to include precision of external measurement
 129   //
 130
 131   if (hh==0)
 132     hh=(nvectors+2)/2;
 133   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};
 134   Int_t *index=new Int_t[nvectors];
 135   TMath::Sort(nvectors, data, index, kFALSE);
 136   //
 137   Int_t    nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11);
 138   Double_t factor = faclts[0];
 139   if (nquant>0){
 140     // fix proper normalization - Anja
 141     factor = faclts[nquant-1];
 142   }
 143
 144   //
 145   //
 146   Double_t sumx  =0;
 147   Double_t sumx2 =0;
 148   Int_t    bestindex = -1;
 149   Double_t bestmean  = 0;
 150   Double_t bestsigma = -1;
 151   for (Int_t i=0; i<hh; i++){
 152     sumx  += data[index[i]];
 153     sumx2 += data[index[i]]*data[index[i]];
 154   }
 155   //
 156   Double_t kfactor = 2.*externalfactor - externalfactor*externalfactor;
 157   Double_t norm = 1./Double_t(hh);
 158   for (Int_t i=hh; i<nvectors; i++){
 159     Double_t cmean  = sumx*norm;
 160     Double_t csigma = (sumx2*norm - cmean*cmean*kfactor);
 161     if (csigma<bestsigma ||  bestsigma<0){
 162       bestmean  = cmean;
 163       bestsigma = csigma;
 164       bestindex = i-hh;
 165     }
 166     //
 167     //
 168     sumx  += data[index[i]]-data[index[i-hh]];
 169     sumx2 += data[index[i]]*data[index[i]]-data[index[i-hh]]*data[index[i-hh]];
 170   }
 171
 172   Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma));
 173   mean  = bestmean;
 174   sigma = bstd;
 175   delete [] index;
 176 }
 177
 178
 179 //_____________________________________________________________________________
 180 Int_t TStatToolkit::Freq(Int_t n, const Int_t *inlist
 181                         , Int_t *outlist, Bool_t down)
 182 {
 183   //
 184   //  Sort eleements according occurancy
 185   //  The size of output array has is 2*n
 186   //
 187
 188   Int_t * sindexS = new Int_t[n];     // temp array for sorting
 189   Int_t * sindexF = new Int_t[2*n];
 190   for (Int_t i=0;i<n;i++) sindexS[i]=0;
 191   for (Int_t i=0;i<2*n;i++) sindexF[i]=0;
 192   //
 193   TMath::Sort(n,inlist, sindexS, down);
 194   Int_t last      = inlist[sindexS[0]];
 195   Int_t val       = last;
 196   sindexF[0]      = 1;
 197   sindexF[0+n]    = last;
 198   Int_t countPos  = 0;
 199   //
 200   //  find frequency
 201   for(Int_t i=1;i<n; i++){
 202     val = inlist[sindexS[i]];
 203     if (last == val)   sindexF[countPos]++;
 204     else{
 205       countPos++;
 206       sindexF[countPos+n] = val;
 207       sindexF[countPos]++;
 208       last =val;
 209     }
 210   }
 211   if (last==val) countPos++;
 212   // sort according frequency
 213   TMath::Sort(countPos, sindexF, sindexS, kTRUE);
 214   for (Int_t i=0;i<countPos;i++){
 215     outlist[2*i  ] = sindexF[sindexS[i]+n];
 216     outlist[2*i+1] = sindexF[sindexS[i]];
 217   }
 218   delete [] sindexS;
 219   delete [] sindexF;
 220
 221   return countPos;
 222
 223 }
 224
 225 //___TStatToolkit__________________________________________________________________________
 226 void TStatToolkit::TruncatedMean(const TH1 * his, TVectorD *param, Float_t down, Float_t up, Bool_t verbose){
 227   //
 228   //
 229   //
 230   Int_t nbins    = his->GetNbinsX();
 231   Float_t nentries = his->GetEntries();
 232   Float_t sum      =0;
 233   Float_t mean   = 0;
 234   Float_t sigma2 = 0;
 235   Float_t ncumul=0;
 236   for (Int_t ibin=1;ibin<nbins; ibin++){
 237     ncumul+= his->GetBinContent(ibin);
 238     Float_t fraction = Float_t(ncumul)/Float_t(nentries);
 239     if (fraction>down && fraction<up){
 240       sum+=his->GetBinContent(ibin);
 241       mean+=his->GetBinCenter(ibin)*his->GetBinContent(ibin);
 242       sigma2+=his->GetBinCenter(ibin)*his->GetBinCenter(ibin)*his->GetBinContent(ibin);
 243     }
 244   }
 245   mean/=sum;
 246   sigma2= TMath::Sqrt(TMath::Abs(sigma2/sum-mean*mean));
 247   if (param){
 248     (*param)[0] = his->GetMaximum();
 249     (*param)[1] = mean;
 250     (*param)[2] = sigma2;
 251
 252   }
 253   if (verbose)  printf("Mean\t%f\t Sigma2\t%f\n", mean,sigma2);
 254 }
 255
 256 void TStatToolkit::LTM(TH1F * his, TVectorD *param , Float_t fraction,  Bool_t verbose){
 257   //
 258   // LTM
 259   //
 260   Int_t nbins    = his->GetNbinsX();
 261   Int_t nentries = (Int_t)his->GetEntries();
 262   Double_t *data  = new Double_t[nentries];
 263   Int_t npoints=0;
 264   for (Int_t ibin=1;ibin<nbins; ibin++){
 265     Float_t entriesI = his->GetBinContent(ibin);
 266     Float_t xcenter= his->GetBinCenter(ibin);
 267     for (Int_t ic=0; ic<entriesI; ic++){
 268       if (npoints<nentries){
 269         data[npoints]= xcenter;
 270         npoints++;
 271       }
 272     }
 273   }
 274   Double_t mean, sigma;
 275   Int_t npoints2=TMath::Min(Int_t(fraction*Float_t(npoints)),npoints-1);
 276   npoints2=TMath::Max(Int_t(0.5*Float_t(npoints)),npoints2);
 277   TStatToolkit::EvaluateUni(npoints, data, mean,sigma,npoints2);
 278   delete [] data;
 279   if (verbose)  printf("Mean\t%f\t Sigma2\t%f\n", mean,sigma);if (param){
 280     (*param)[0] = his->GetMaximum();
 281     (*param)[1] = mean;
 282     (*param)[2] = sigma;
 283   }
 284 }
 285
 286 Double_t  TStatToolkit::FitGaus(TH1* his, TVectorD *param, TMatrixD */*matrix*/, Float_t xmin, Float_t xmax, Bool_t verbose){
 287   //
 288   //  Fit histogram with gaussian function
 289   //
 290   //  Prameters:
 291   //       return value- chi2 - if negative ( not enough points)
 292   //       his        -  input histogram
 293   //       param      -  vector with parameters
 294   //       xmin, xmax -  range to fit - if xmin=xmax=0 - the full histogram range used
 295   //  Fitting:
 296   //  1. Step - make logarithm
 297   //  2. Linear  fit (parabola) - more robust - always converge
 298   //  3. In case of small statistic bins are averaged
 299   //
 300   static TLinearFitter fitter(3,"pol2");
 301   TVectorD  par(3);
 302   TVectorD  sigma(3);
 303   TMatrixD mat(3,3);
 304   if (his->GetMaximum()<4) return -1;
 305   if (his->GetEntries()<12) return -1;
 306   if (his->GetRMS()<mat.GetTol()) return -1;
 307   Float_t maxEstimate   = his->GetEntries()*his->GetBinWidth(1)/TMath::Sqrt((TMath::TwoPi()*his->GetRMS()));
 308   Int_t dsmooth = TMath::Nint(6./TMath::Sqrt(maxEstimate));
 309
 310   if (maxEstimate<1) return -1;
 311   Int_t nbins    = his->GetNbinsX();
 312   Int_t npoints=0;
 313   //
 314
 315
 316   if (xmin>=xmax){
 317     xmin = his->GetXaxis()->GetXmin();
 318     xmax = his->GetXaxis()->GetXmax();
 319   }
 320   for (Int_t iter=0; iter<2; iter++){
 321     fitter.ClearPoints();
 322     npoints=0;
 323     for (Int_t ibin=1;ibin<nbins+1; ibin++){
 324       Int_t countB=1;
 325       Float_t entriesI =  his->GetBinContent(ibin);
 326       for (Int_t delta = -dsmooth; delta<=dsmooth; delta++){
 327         if (ibin+delta>1 &&ibin+delta<nbins-1){
 328           entriesI +=  his->GetBinContent(ibin+delta);
 329           countB++;
 330         }
 331       }
 332       entriesI/=countB;
 333       Double_t xcenter= his->GetBinCenter(ibin);
 334       if (xcenter<xmin || xcenter>xmax) continue;
 335       Double_t error=1./TMath::Sqrt(countB);
 336       Float_t   cont=2;
 337       if (iter>0){
 338         if (par[0]+par[1]*xcenter+par[2]*xcenter*xcenter>20) return 0;
 339         cont = TMath::Exp(par[0]+par[1]*xcenter+par[2]*xcenter*xcenter);
 340         if (cont>1.) error = 1./TMath::Sqrt(cont*Float_t(countB));
 341       }
 342       if (entriesI>1&&cont>1){
 343         fitter.AddPoint(&xcenter,TMath::Log(Float_t(entriesI)),error);
 344         npoints++;
 345       }
 346     }
 347     if (npoints>3){
 348       fitter.Eval();
 349       fitter.GetParameters(par);
 350     }else{
 351       break;
 352     }
 353   }
 354   if (npoints<=3){
 355     return -1;
 356   }
 357   fitter.GetParameters(par);
 358   fitter.GetCovarianceMatrix(mat);
 359   if (TMath::Abs(par[1])<mat.GetTol()) return -1;
 360   if (TMath::Abs(par[2])<mat.GetTol()) return -1;
 361   Double_t chi2 = fitter.GetChisquare()/Float_t(npoints);
 362   //fitter.GetParameters();
 363   if (!param)  param  = new TVectorD(3);
 364   // if (!matrix) matrix = new TMatrixD(3,3); // Covariance matrix to be implemented
 365   (*param)[1] = par[1]/(-2.*par[2]);
 366   (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
 367   (*param)[0] = TMath::Exp(par[0]+ par[1]* (*param)[1] +  par[2]*(*param)[1]*(*param)[1]);
 368   if (verbose){
 369     par.Print();
 370     mat.Print();
 371     param->Print();
 372     printf("Chi2=%f\n",chi2);
 373     TF1 * f1= new TF1("f1","[0]*exp(-(x-[1])^2/(2*[2]*[2]))",his->GetXaxis()->GetXmin(),his->GetXaxis()->GetXmax());
 374     f1->SetParameter(0, (*param)[0]);
 375     f1->SetParameter(1, (*param)[1]);
 376     f1->SetParameter(2, (*param)[2]);
 377     f1->Draw("same");
 378   }
 379   return chi2;
 380 }
 381
 382 Double_t  TStatToolkit::FitGaus(Float_t *arr, Int_t nBins, Float_t xMin, Float_t xMax, TVectorD *param, TMatrixD */*matrix*/, Bool_t verbose){
 383   //
 384   //  Fit histogram with gaussian function
 385   //
 386   //  Prameters:
 387   //     nbins: size of the array and number of histogram bins
 388   //     xMin, xMax: histogram range
 389   //     param: paramters of the fit (0-Constant, 1-Mean, 2-Sigma)
 390   //     matrix: covariance matrix -- not implemented yet, pass dummy matrix!!!
 391   //
 392   //  Return values:
 393   //    >0: the chi2 returned by TLinearFitter
 394   //    -3: only three points have been used for the calculation - no fitter was used
 395   //    -2: only two points have been used for the calculation - center of gravity was uesed for calculation
 396   //    -1: only one point has been used for the calculation - center of gravity was uesed for calculation
 397   //    -4: invalid result!!
 398   //
 399   //  Fitting:
 400   //  1. Step - make logarithm
 401   //  2. Linear  fit (parabola) - more robust - always converge
 402   //
 403   static TLinearFitter fitter(3,"pol2");
 404   static TMatrixD mat(3,3);
 405   static Double_t kTol = mat.GetTol();
 406   fitter.StoreData(kFALSE);
 407   fitter.ClearPoints();
 408   TVectorD  par(3);
 409   TVectorD  sigma(3);
 410   TMatrixD matA(3,3);
 411   TMatrixD b(3,1);
 412   Float_t rms = TMath::RMS(nBins,arr);
 413   Float_t max = TMath::MaxElement(nBins,arr);
 414   Float_t binWidth = (xMax-xMin)/(Float_t)nBins;
 415
 416   Float_t meanCOG = 0;
 417   Float_t rms2COG = 0;
 418   Float_t sumCOG  = 0;
 419
 420   Float_t entries = 0;
 421   Int_t nfilled=0;
 422
 423   for (Int_t i=0; i<nBins; i++){
 424       entries+=arr[i];
 425       if (arr[i]>0) nfilled++;
 426   }
 427
 428   if (max<4) return -4;
 429   if (entries<12) return -4;
 430   if (rms<kTol) return -4;
 431
 432   Int_t npoints=0;
 433   //
 434
 435   //
 436   for (Int_t ibin=0;ibin<nBins; ibin++){
 437       Float_t entriesI = arr[ibin];
 438     if (entriesI>1){
 439       Double_t xcenter = xMin+(ibin+0.5)*binWidth;
 440
 441       Float_t error    = 1./TMath::Sqrt(entriesI);
 442       Float_t val = TMath::Log(Float_t(entriesI));
 443       fitter.AddPoint(&xcenter,val,error);
 444       if (npoints<3){
 445           matA(npoints,0)=1;
 446           matA(npoints,1)=xcenter;
 447           matA(npoints,2)=xcenter*xcenter;
 448           b(npoints,0)=val;
 449           meanCOG+=xcenter*entriesI;
 450           rms2COG +=xcenter*entriesI*xcenter;
 451           sumCOG +=entriesI;
 452       }
 453       npoints++;
 454     }
 455   }
 456
 457
 458   Double_t chi2 = 0;
 459   if (npoints>=3){
 460       if ( npoints == 3 ){
 461           //analytic calculation of the parameters for three points
 462           matA.Invert();
 463           TMatrixD res(1,3);
 464           res.Mult(matA,b);
 465           par[0]=res(0,0);
 466           par[1]=res(0,1);
 467           par[2]=res(0,2);
 468           chi2 = -3.;
 469       } else {
 470           // use fitter for more than three points
 471           fitter.Eval();
 472           fitter.GetParameters(par);
 473           fitter.GetCovarianceMatrix(mat);
 474           chi2 = fitter.GetChisquare()/Float_t(npoints);
 475       }
 476       if (TMath::Abs(par[1])<kTol) return -4;
 477       if (TMath::Abs(par[2])<kTol) return -4;
 478
 479       if (!param)  param  = new TVectorD(3);
 480       //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
 481
 482       (*param)[1] = par[1]/(-2.*par[2]);
 483       (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
 484       Double_t lnparam0 = par[0]+ par[1]* (*param)[1] +  par[2]*(*param)[1]*(*param)[1];
 485       if ( lnparam0>307 ) return -4;
 486       (*param)[0] = TMath::Exp(lnparam0);
 487       if (verbose){
 488           par.Print();
 489           mat.Print();
 490           param->Print();
 491           printf("Chi2=%f\n",chi2);
 492           TF1 * f1= new TF1("f1","[0]*exp(-(x-[1])^2/(2*[2]*[2]))",xMin,xMax);
 493           f1->SetParameter(0, (*param)[0]);
 494           f1->SetParameter(1, (*param)[1]);
 495           f1->SetParameter(2, (*param)[2]);
 496           f1->Draw("same");
 497       }
 498       return chi2;
 499   }
 500
 501   if (npoints == 2){
 502       //use center of gravity for 2 points
 503       meanCOG/=sumCOG;
 504       rms2COG /=sumCOG;
 505       (*param)[0] = max;
 506       (*param)[1] = meanCOG;
 507       (*param)[2] = TMath::Sqrt(TMath::Abs(meanCOG*meanCOG-rms2COG));
 508       chi2=-2.;
 509   }
 510   if ( npoints == 1 ){
 511       meanCOG/=sumCOG;
 512       (*param)[0] = max;
 513       (*param)[1] = meanCOG;
 514       (*param)[2] = binWidth/TMath::Sqrt(12);
 515       chi2=-1.;
 516   }
 517   return chi2;
 518
 519 }
 520
 521
 522 Float_t TStatToolkit::GetCOG(const Short_t *arr, Int_t nBins, Float_t xMin, Float_t xMax, Float_t *rms, Float_t *sum)
 523 {
 524     //
 525     //  calculate center of gravity rms and sum for array 'arr' with nBins an a x range xMin to xMax
 526     //  return COG; in case of failure return xMin
 527     //
 528     Float_t meanCOG = 0;
 529     Float_t rms2COG = 0;
 530     Float_t sumCOG  = 0;
 531     Int_t npoints   = 0;
 532
 533     Float_t binWidth = (xMax-xMin)/(Float_t)nBins;
 534
 535     for (Int_t ibin=0; ibin<nBins; ibin++){
 536         Float_t entriesI = (Float_t)arr[ibin];
 537         Double_t xcenter = xMin+(ibin+0.5)*binWidth;
 538         if ( entriesI>0 ){
 539             meanCOG += xcenter*entriesI;
 540             rms2COG += xcenter*entriesI*xcenter;
 541             sumCOG  += entriesI;
 542             npoints++;
 543         }
 544     }
 545     if ( sumCOG == 0 ) return xMin;
 546     meanCOG/=sumCOG;
 547
 548     if ( rms ){
 549         rms2COG /=sumCOG;
 550         (*rms) = TMath::Sqrt(TMath::Abs(meanCOG*meanCOG-rms2COG));
 551         if ( npoints == 1 ) (*rms) = binWidth/TMath::Sqrt(12);
 552     }
 553
 554     if ( sum )
 555         (*sum) = sumCOG;
 556
 557     return meanCOG;
 558 }
 559
 560
 561
 562 ///////////////////////////////////////////////////////////////
 563 //////////////         TEST functions /////////////////////////
 564 ///////////////////////////////////////////////////////////////
 565
 566
 567
 568
 569
 570 void TStatToolkit::TestGausFit(Int_t nhistos){
 571   //
 572   // Test performance of the parabolic - gaussian fit - compare it with
 573   // ROOT gauss fit
 574   //  nhistos - number of histograms to be used for test
 575   //
 576   TTreeSRedirector *pcstream = new TTreeSRedirector("fitdebug.root");
 577
 578   Float_t  *xTrue = new Float_t[nhistos];
 579   Float_t  *sTrue = new Float_t[nhistos];
 580   TVectorD **par1  = new TVectorD*[nhistos];
 581   TVectorD **par2  = new TVectorD*[nhistos];
 582   TMatrixD dummy(3,3);
 583
 584
 585   TH1F **h1f = new TH1F*[nhistos];
 586   TF1  *myg = new TF1("myg","gaus");
 587   TF1  *fit = new TF1("fit","gaus");
 588   gRandom->SetSeed(0);
 589
 590   //init
 591   for (Int_t i=0;i<nhistos; i++){
 592     par1[i] = new TVectorD(3);
 593     par2[i] = new TVectorD(3);
 594     h1f[i]  = new TH1F(Form("h1f%d",i),Form("h1f%d",i),20,-10,10);
 595     xTrue[i]= gRandom->Rndm();
 596     gSystem->Sleep(2);
 597     sTrue[i]= .75+gRandom->Rndm()*.5;
 598     myg->SetParameters(1,xTrue[i],sTrue[i]);
 599     h1f[i]->FillRandom("myg");
 600   }
 601
 602   TStopwatch s;
 603   s.Start();
 604   //standard gaus fit
 605   for (Int_t i=0; i<nhistos; i++){
 606     h1f[i]->Fit(fit,"0q");
 607     (*par1[i])(0) = fit->GetParameter(0);
 608     (*par1[i])(1) = fit->GetParameter(1);
 609     (*par1[i])(2) = fit->GetParameter(2);
 610   }
 611   s.Stop();
 612   printf("Gaussian fit\t");
 613   s.Print();
 614
 615   s.Start();
 616   //TStatToolkit gaus fit
 617   for (Int_t i=0; i<nhistos; i++){
 618     TStatToolkit::FitGaus(h1f[i]->GetArray()+1,h1f[i]->GetNbinsX(),h1f[i]->GetXaxis()->GetXmin(),h1f[i]->GetXaxis()->GetXmax(),par2[i],&dummy);
 619   }
 620
 621   s.Stop();
 622   printf("Parabolic fit\t");
 623   s.Print();
 624   //write stream
 625   for (Int_t i=0;i<nhistos; i++){
 626     Float_t xt  = xTrue[i];
 627     Float_t st  = sTrue[i];
 628     (*pcstream)<<"data"
 629                <<"xTrue="<<xt
 630                <<"sTrue="<<st
 631                <<"pg.="<<(par1[i])
 632                <<"pa.="<<(par2[i])
 633                <<"\n";
 634   }
 635   //delete pointers
 636   for (Int_t i=0;i<nhistos; i++){
 637     delete par1[i];
 638     delete par2[i];
 639     delete h1f[i];
 640   }
 641   delete pcstream;
 642   delete []h1f;
 643   delete []xTrue;
 644   delete []sTrue;
 645   //
 646   delete []par1;
 647   delete []par2;
 648
 649 }
 650
 651
 652
 653 TGraph2D * TStatToolkit::MakeStat2D(TH3 * his, Int_t delta0, Int_t delta1, Int_t type){
 654   //
 655   //
 656   //
 657   // delta - number of bins to integrate
 658   // type - 0 - mean value
 659
 660   TAxis * xaxis  = his->GetXaxis();
 661   TAxis * yaxis  = his->GetYaxis();
 662   //  TAxis * zaxis  = his->GetZaxis();
 663   Int_t   nbinx  = xaxis->GetNbins();
 664   Int_t   nbiny  = yaxis->GetNbins();
 665   char name[1000];
 666   Int_t icount=0;
 667   TGraph2D  *graph = new TGraph2D(nbinx*nbiny);
 668   TF1 f1("f1","gaus");
 669   for (Int_t ix=0; ix<nbinx;ix++)
 670     for (Int_t iy=0; iy<nbiny;iy++){
 671       Float_t xcenter = xaxis->GetBinCenter(ix);
 672       Float_t ycenter = yaxis->GetBinCenter(iy);
 673       snprintf(name,1000,"%s_%d_%d",his->GetName(), ix,iy);
 674       TH1 *projection = his->ProjectionZ(name,ix-delta0,ix+delta0,iy-delta1,iy+delta1);
 675       Float_t stat= 0;
 676       if (type==0) stat = projection->GetMean();
 677       if (type==1) stat = projection->GetRMS();
 678       if (type==2 || type==3){
 679         TVectorD vec(3);
 680         TStatToolkit::LTM((TH1F*)projection,&vec,0.7);
 681         if (type==2) stat= vec[1];
 682         if (type==3) stat= vec[0];
 683       }
 684       if (type==4|| type==5){
 685         projection->Fit(&f1);
 686         if (type==4) stat= f1.GetParameter(1);
 687         if (type==5) stat= f1.GetParameter(2);
 688       }
 689       //printf("%d\t%f\t%f\t%f\n", icount,xcenter, ycenter, stat);
 690       graph->SetPoint(icount,xcenter, ycenter, stat);
 691       icount++;
 692     }
 693   return graph;
 694 }
 695
 696 TGraph * TStatToolkit::MakeStat1D(TH3 * his, Int_t delta1, Int_t type){
 697   //
 698   //
 699   //
 700   // delta - number of bins to integrate
 701   // type - 0 - mean value
 702
 703   TAxis * xaxis  = his->GetXaxis();
 704   TAxis * yaxis  = his->GetYaxis();
 705   //  TAxis * zaxis  = his->GetZaxis();
 706   Int_t   nbinx  = xaxis->GetNbins();
 707   Int_t   nbiny  = yaxis->GetNbins();
 708   char name[1000];
 709   Int_t icount=0;
 710   TGraph  *graph = new TGraph(nbinx);
 711   TF1 f1("f1","gaus");
 712   for (Int_t ix=0; ix<nbinx;ix++){
 713     Float_t xcenter = xaxis->GetBinCenter(ix);
 714     //    Float_t ycenter = yaxis->GetBinCenter(iy);
 715     snprintf(name,1000,"%s_%d",his->GetName(), ix);
 716     TH1 *projection = his->ProjectionZ(name,ix-delta1,ix+delta1,0,nbiny);
 717     Float_t stat= 0;
 718     if (type==0) stat = projection->GetMean();
 719     if (type==1) stat = projection->GetRMS();
 720     if (type==2 || type==3){
 721       TVectorD vec(3);
 722         TStatToolkit::LTM((TH1F*)projection,&vec,0.7);
 723         if (type==2) stat= vec[1];
 724         if (type==3) stat= vec[0];
 725     }
 726     if (type==4|| type==5){
 727       projection->Fit(&f1);
 728       if (type==4) stat= f1.GetParameter(1);
 729       if (type==5) stat= f1.GetParameter(2);
 730     }
 731       //printf("%d\t%f\t%f\t%f\n", icount,xcenter, ycenter, stat);
 732     graph->SetPoint(icount,xcenter, stat);
 733     icount++;
 734   }
 735   return graph;
 736 }
 737
 738
 739
 740
 741
 742 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){
 743    //
 744    // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts
 745    // returns chi2, fitParam and covMatrix
 746    // returns TString with fitted formula
 747    //
 748
 749    TString formulaStr(formula);
 750    TString drawStr(drawCommand);
 751    TString cutStr(cuts);
 752    TString ferr("1");
 753
 754    TString strVal(drawCommand);
 755    if (strVal.Contains(":")){
 756      TObjArray* valTokens = strVal.Tokenize(":");
 757      drawStr = valTokens->At(0)->GetName();
 758      ferr       = valTokens->At(1)->GetName();
 759      delete valTokens;
 760    }
 761
 762
 763    formulaStr.ReplaceAll("++", "~");
 764    TObjArray* formulaTokens = formulaStr.Tokenize("~");
 765    Int_t dim = formulaTokens->GetEntriesFast();
 766
 767    fitParam.ResizeTo(dim);
 768    covMatrix.ResizeTo(dim,dim);
 769
 770    TLinearFitter* fitter = new TLinearFitter(dim+1, Form("hyp%d",dim));
 771    fitter->StoreData(kTRUE);
 772    fitter->ClearPoints();
 773
 774    Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff",  stop-start, start);
 775    if (entries == -1) {
 776      delete formulaTokens;
 777      return new TString("An ERROR has occured during fitting!");
 778    }
 779    Double_t **values = new Double_t*[dim+1] ;
 780    for (Int_t i=0; i<dim+1; i++) values[i]=NULL;
 781    //
 782    entries = tree->Draw(ferr.Data(), cutStr.Data(), "goff",  stop-start, start);
 783    if (entries == -1) {
 784      delete formulaTokens;
 785      delete []values;
 786      return new TString("An ERROR has occured during fitting!");
 787    }
 788    Double_t *errors = new Double_t[entries];
 789    memcpy(errors,  tree->GetV1(), entries*sizeof(Double_t));
 790
 791    for (Int_t i = 0; i < dim + 1; i++){
 792       Int_t centries = 0;
 793       if (i < dim) centries = tree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start);
 794       else  centries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start);
 795
 796       if (entries != centries) {
 797         delete []errors;
 798         delete []values;
 799         return new TString("An ERROR has occured during fitting!");
 800       }
 801       values[i] = new Double_t[entries];
 802       memcpy(values[i],  tree->GetV1(), entries*sizeof(Double_t));
 803    }
 804
 805    // add points to the fitter
 806    for (Int_t i = 0; i < entries; i++){
 807       Double_t x[1000];
 808       for (Int_t j=0; j<dim;j++) x[j]=values[j][i];
 809       fitter->AddPoint(x, values[dim][i], errors[i]);
 810    }
 811
 812    fitter->Eval();
 813    if (frac>0.5 && frac<1){
 814      fitter->EvalRobust(frac);
 815    }else{
 816      if (fix0) {
 817        fitter->FixParameter(0,0);
 818        fitter->Eval();
 819      }
 820    }
 821    fitter->GetParameters(fitParam);
 822    fitter->GetCovarianceMatrix(covMatrix);
 823    chi2 = fitter->GetChisquare();
 824    npoints = entries;
 825    TString *preturnFormula = new TString(Form("( %f+",fitParam[0])), &returnFormula = *preturnFormula;
 826
 827    for (Int_t iparam = 0; iparam < dim; iparam++) {
 828      returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam+1]));
 829      if (iparam < dim-1) returnFormula.Append("+");
 830    }
 831    returnFormula.Append(" )");
 832
 833
 834    for (Int_t j=0; j<dim+1;j++) delete [] values[j];
 835
 836
 837    delete formulaTokens;
 838    delete fitter;
 839    delete[] values;
 840    delete[] errors;
 841    return preturnFormula;
 842 }
 843
 844 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){
 845    //
 846    // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts
 847    // returns chi2, fitParam and covMatrix
 848    // returns TString with fitted formula
 849    //
 850
 851    TString formulaStr(formula);
 852    TString drawStr(drawCommand);
 853    TString cutStr(cuts);
 854    TString ferr("1");
 855
 856    TString strVal(drawCommand);
 857    if (strVal.Contains(":")){
 858      TObjArray* valTokens = strVal.Tokenize(":");
 859      drawStr = valTokens->At(0)->GetName();
 860      ferr       = valTokens->At(1)->GetName();
 861      delete valTokens;
 862    }
 863
 864
 865    formulaStr.ReplaceAll("++", "~");
 866    TObjArray* formulaTokens = formulaStr.Tokenize("~");
 867    Int_t dim = formulaTokens->GetEntriesFast();
 868
 869    fitParam.ResizeTo(dim);
 870    covMatrix.ResizeTo(dim,dim);
 871
 872    TLinearFitter* fitter = new TLinearFitter(dim+1, Form("hyp%d",dim));
 873    fitter->StoreData(kTRUE);
 874    fitter->ClearPoints();
 875
 876    Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff",  stop-start, start);
 877    if (entries == -1) {
 878      delete formulaTokens;
 879      return new TString("An ERROR has occured during fitting!");
 880    }
 881    Double_t **values = new Double_t*[dim+1] ;
 882    for (Int_t i=0; i<dim+1; i++) values[i]=NULL;
 883    //
 884    entries = tree->Draw(ferr.Data(), cutStr.Data(), "goff",  stop-start, start);
 885    if (entries == -1) {
 886      delete formulaTokens;
 887      delete [] values;
 888      return new TString("An ERROR has occured during fitting!");
 889    }
 890    Double_t *errors = new Double_t[entries];
 891    memcpy(errors,  tree->GetV1(), entries*sizeof(Double_t));
 892
 893    for (Int_t i = 0; i < dim + 1; i++){
 894       Int_t centries = 0;
 895       if (i < dim) centries = tree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start);
 896       else  centries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start);
 897
 898       if (entries != centries) {
 899         delete []errors;
 900         delete []values;
 901         delete formulaTokens;
 902         return new TString("An ERROR has occured during fitting!");
 903       }
 904       values[i] = new Double_t[entries];
 905       memcpy(values[i],  tree->GetV1(), entries*sizeof(Double_t));
 906    }
 907
 908    // add points to the fitter
 909    for (Int_t i = 0; i < entries; i++){
 910       Double_t x[1000];
 911       for (Int_t j=0; j<dim;j++) x[j]=values[j][i];
 912       fitter->AddPoint(x, values[dim][i], errors[i]);
 913    }
 914    if (constrain>0){
 915      for (Int_t i = 0; i < dim; i++){
 916        Double_t x[1000];
 917        for (Int_t j=0; j<dim;j++) if (i!=j) x[j]=0;
 918        x[i]=1.;
 919        fitter->AddPoint(x, 0, constrain);
 920      }
 921    }
 922
 923
 924    fitter->Eval();
 925    if (frac>0.5 && frac<1){
 926      fitter->EvalRobust(frac);
 927    }
 928    fitter->GetParameters(fitParam);
 929    fitter->GetCovarianceMatrix(covMatrix);
 930    chi2 = fitter->GetChisquare();
 931    npoints = entries;
 932
 933    TString *preturnFormula = new TString(Form("( %f+",fitParam[0])), &returnFormula = *preturnFormula;
 934
 935    for (Int_t iparam = 0; iparam < dim; iparam++) {
 936      returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam+1]));
 937      if (iparam < dim-1) returnFormula.Append("+");
 938    }
 939    returnFormula.Append(" )");
 940
 941    for (Int_t j=0; j<dim+1;j++) delete [] values[j];
 942
 943
 944
 945    delete formulaTokens;
 946    delete fitter;
 947    delete[] values;
 948    delete[] errors;
 949    return preturnFormula;
 950 }
 951
 952
 953
 954 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){
 955    //
 956    // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts
 957    // returns chi2, fitParam and covMatrix
 958    // returns TString with fitted formula
 959    //
 960
 961    TString formulaStr(formula);
 962    TString drawStr(drawCommand);
 963    TString cutStr(cuts);
 964    TString ferr("1");
 965
 966    TString strVal(drawCommand);
 967    if (strVal.Contains(":")){
 968      TObjArray* valTokens = strVal.Tokenize(":");
 969      drawStr = valTokens->At(0)->GetName();
 970      ferr       = valTokens->At(1)->GetName();
 971      delete valTokens;
 972    }
 973
 974
 975    formulaStr.ReplaceAll("++", "~");
 976    TObjArray* formulaTokens = formulaStr.Tokenize("~");
 977    Int_t dim = formulaTokens->GetEntriesFast();
 978
 979    fitParam.ResizeTo(dim);
 980    covMatrix.ResizeTo(dim,dim);
 981    TString fitString="x0";
 982    for (Int_t i=1; i<dim; i++) fitString+=Form("++x%d",i);
 983    TLinearFitter* fitter = new TLinearFitter(dim, fitString.Data());
 984    fitter->StoreData(kTRUE);
 985    fitter->ClearPoints();
 986
 987    Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff",  stop-start, start);
 988    if (entries == -1) {
 989      delete formulaTokens;
 990      return new TString("An ERROR has occured during fitting!");
 991    }
 992    Double_t **values = new Double_t*[dim+1] ;
 993    for (Int_t i=0; i<dim+1; i++) values[i]=NULL;
 994    //
 995    entries = tree->Draw(ferr.Data(), cutStr.Data(), "goff",  stop-start, start);
 996    if (entries == -1) {
 997      delete []values;
 998      delete formulaTokens;
 999      return new TString("An ERROR has occured during fitting!");
1000    }
1001    Double_t *errors = new Double_t[entries];
1002    memcpy(errors,  tree->GetV1(), entries*sizeof(Double_t));
1003
1004    for (Int_t i = 0; i < dim + 1; i++){
1005       Int_t centries = 0;
1006       if (i < dim) centries = tree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start);
1007       else  centries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start);
1008
1009       if (entries != centries) {
1010         delete []errors;
1011         delete []values;
1012         delete formulaTokens;
1013         return new TString("An ERROR has occured during fitting!");
1014       }
1015       values[i] = new Double_t[entries];
1016       memcpy(values[i],  tree->GetV1(), entries*sizeof(Double_t));
1017    }
1018
1019    // add points to the fitter
1020    for (Int_t i = 0; i < entries; i++){
1021       Double_t x[1000];
1022       for (Int_t j=0; j<dim;j++) x[j]=values[j][i];
1023       fitter->AddPoint(x, values[dim][i], errors[i]);
1024    }
1025
1026    fitter->Eval();
1027    if (frac>0.5 && frac<1){
1028      fitter->EvalRobust(frac);
1029    }
1030    fitter->GetParameters(fitParam);
1031    fitter->GetCovarianceMatrix(covMatrix);
1032    chi2 = fitter->GetChisquare();
1033    npoints = entries;
1034
1035    TString *preturnFormula = new TString("("), &returnFormula = *preturnFormula;
1036
1037    for (Int_t iparam = 0; iparam < dim; iparam++) {
1038      returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam]));
1039      if (iparam < dim-1) returnFormula.Append("+");
1040    }
1041    returnFormula.Append(" )");
1042
1043
1044    for (Int_t j=0; j<dim+1;j++) delete [] values[j];
1045
1046    delete formulaTokens;
1047    delete fitter;
1048    delete[] values;
1049    delete[] errors;
1050    return preturnFormula;
1051 }
1052
1053
1054
1055
1056
1057 Int_t TStatToolkit::GetFitIndex(const TString fString, const TString subString){
1058   //
1059   // fitString - ++ separated list of fits
1060   // substring - ++ separated list of the requiered substrings
1061   //
1062   // return the last occurance of substring in fit string
1063   //
1064   TObjArray *arrFit = fString.Tokenize("++");
1065   TObjArray *arrSub = subString.Tokenize("++");
1066   Int_t index=-1;
1067   for (Int_t i=0; i<arrFit->GetEntries(); i++){
1068     Bool_t isOK=kTRUE;
1069     TString str =arrFit->At(i)->GetName();
1070     for (Int_t isub=0; isub<arrSub->GetEntries(); isub++){
1071       if (str.Contains(arrSub->At(isub)->GetName())==0) isOK=kFALSE;
1072     }
1073     if (isOK) index=i;
1074   }
1075   delete arrFit;
1076   delete arrSub;
1077   return index;
1078 }
1079
1080
1081 TString  TStatToolkit::FilterFit(const TString &input, const TString filter, TVectorD &param, TMatrixD & covar){
1082   //
1083   // Filter fit expression make sub-fit
1084   //
1085   TObjArray *array0= input.Tokenize("++");
1086   TObjArray *array1= filter.Tokenize("++");
1087   //TString *presult=new TString("(0");
1088   TString result="(0.0";
1089   for (Int_t i=0; i<array0->GetEntries(); i++){
1090     Bool_t isOK=kTRUE;
1091     TString str(array0->At(i)->GetName());
1092     for (Int_t j=0; j<array1->GetEntries(); j++){
1093       if (str.Contains(array1->At(j)->GetName())==0) isOK=kFALSE;
1094     }
1095     if (isOK) {
1096       result+="+"+str;
1097       result+=Form("*(%f)",param[i+1]);
1098       printf("%f\t%f\t%s\n",param[i+1], TMath::Sqrt(covar(i+1,i+1)),str.Data());
1099     }
1100   }
1101   result+="-0.)";
1102   delete array0;
1103   delete array1;
1104   return result;
1105 }
1106
1107 void TStatToolkit::Update1D(Double_t delta, Double_t sigma, Int_t s1, TMatrixD &vecXk, TMatrixD &covXk){
1108   //
1109   // Update parameters and covariance - with one measurement
1110   // Input:
1111   // vecXk - input vector - Updated in function
1112   // covXk - covariance matrix - Updated in function
1113   // delta, sigma, s1 - new measurement, rms of new measurement and the index of measurement
1114   const Int_t knMeas=1;
1115   Int_t knElem=vecXk.GetNrows();
1116
1117   TMatrixD mat1(knElem,knElem);            // update covariance matrix
1118   TMatrixD matHk(1,knElem);        // vector to mesurement
1119   TMatrixD vecYk(knMeas,1);        // Innovation or measurement residual
1120   TMatrixD matHkT(knElem,knMeas);  // helper matrix Hk transpose
1121   TMatrixD matSk(knMeas,knMeas);   // Innovation (or residual) covariance
1122   TMatrixD matKk(knElem,knMeas);   // Optimal Kalman gain
1123   TMatrixD covXk2(knElem,knElem);  // helper matrix
1124   TMatrixD covXk3(knElem,knElem);  // helper matrix
1125   TMatrixD vecZk(1,1);
1126   TMatrixD measR(1,1);
1127   vecZk(0,0)=delta;
1128   measR(0,0)=sigma*sigma;
1129   //
1130   // reset matHk
1131   for (Int_t iel=0;iel<knElem;iel++)
1132     for (Int_t ip=0;ip<knMeas;ip++) matHk(ip,iel)=0;
1133   //mat1
1134   for (Int_t iel=0;iel<knElem;iel++) {
1135     for (Int_t jel=0;jel<knElem;jel++) mat1(iel,jel)=0;
1136     mat1(iel,iel)=1;
1137   }
1138   //
1139   matHk(0, s1)=1;
1140   vecYk = vecZk-matHk*vecXk;               // Innovation or measurement residual
1141   matHkT=matHk.T(); matHk.T();
1142   matSk = (matHk*(covXk*matHkT))+measR;    // Innovation (or residual) covariance
1143   matSk.Invert();
1144   matKk = (covXk*matHkT)*matSk;            //  Optimal Kalman gain
1145   vecXk += matKk*vecYk;                    //  updated vector
1146   covXk2= (mat1-(matKk*matHk));
1147   covXk3 =  covXk2*covXk;
1148   covXk = covXk3;
1149   Int_t nrows=covXk3.GetNrows();
1150
1151   for (Int_t irow=0; irow<nrows; irow++)
1152     for (Int_t icol=0; icol<nrows; icol++){
1153       // rounding problems - make matrix again symteric
1154       covXk(irow,icol)=(covXk3(irow,icol)+covXk3(icol,irow))*0.5;
1155     }
1156 }
1157
1158
1159
1160 void   TStatToolkit::Constrain1D(const TString &input, const TString filter, TVectorD &param, TMatrixD & covar, Double_t mean, Double_t sigma){
1161   //
1162   // constrain linear fit
1163   // input  - string description of fit function
1164   // filter - string filter to select sub fits
1165   // param,covar - parameters and covariance matrix of the fit
1166   // mean,sigma  - new measurement uning which the fit is updated
1167   //
1168
1169   TObjArray *array0= input.Tokenize("++");
1170   TObjArray *array1= filter.Tokenize("++");
1171   TMatrixD paramM(param.GetNrows(),1);
1172   for (Int_t i=0; i<=array0->GetEntries(); i++){paramM(i,0)=param(i);}
1173
1174   if (filter.Length()==0){
1175     TStatToolkit::Update1D(mean, sigma, 0, paramM, covar);//
1176   }else{
1177     for (Int_t i=0; i<array0->GetEntries(); i++){
1178       Bool_t isOK=kTRUE;
1179       TString str(array0->At(i)->GetName());
1180       for (Int_t j=0; j<array1->GetEntries(); j++){
1181         if (str.Contains(array1->At(j)->GetName())==0) isOK=kFALSE;
1182       }
1183       if (isOK) {
1184         TStatToolkit::Update1D(mean, sigma, i+1, paramM, covar);//
1185       }
1186     }
1187   }
1188   for (Int_t i=0; i<=array0->GetEntries(); i++){
1189     param(i)=paramM(i,0);
1190   }
1191   delete array0;
1192   delete array1;
1193 }
1194
1195 TString  TStatToolkit::MakeFitString(const TString &input, const TVectorD &param, const TMatrixD & covar, Bool_t verbose){
1196   //
1197   //
1198   //
1199   TObjArray *array0= input.Tokenize("++");
1200   TString result=Form("(%f",param[0]);
1201   printf("%f\t%f\t\n", param[0], TMath::Sqrt(covar(0,0)));
1202   for (Int_t i=0; i<array0->GetEntries(); i++){
1203     TString str(array0->At(i)->GetName());
1204     result+="+"+str;
1205     result+=Form("*(%f)",param[i+1]);
1206     if (verbose) printf("%f\t%f\t%s\n", param[i+1], TMath::Sqrt(covar(i+1,i+1)),str.Data());
1207   }
1208   result+="-0.)";
1209   delete array0;
1210   return result;
1211 }
1212
1213
1214 TGraph * TStatToolkit::MakeGraphSparse(TTree * tree, const char * expr, const char * cut, Int_t mstyle, Int_t mcolor, Float_t msize){
1215   //
1216   // Make a sparse draw of the variables
1217   // Writen by Weilin.Yu
1218   const Int_t entries =  tree->Draw(expr,cut,"goff");
1219   if (entries<=0) {
1220     TStatToolkit t;
1221     t.Error("TStatToolkit::MakeGraphSparse",Form("Empty or Not valid expression (%s) or cut *%s)", expr,cut));
1222     return 0;
1223   }
1224   //  TGraph * graph = (TGraph*)gPad->GetPrimitive("Graph"); // 2D
1225   TGraph * graph = 0;
1226   if (tree->GetV3()) graph = new TGraphErrors (entries, tree->GetV2(),tree->GetV1(),0,tree->GetV3());
1227   graph =  new TGraphErrors (entries, tree->GetV2(),tree->GetV1(),0,0);
1228   graph->SetMarkerStyle(mstyle);
1229   graph->SetMarkerColor(mcolor);
1230   //
1231   Int_t *index = new Int_t[entries*4];
1232   TMath::Sort(entries,graph->GetX(),index,kFALSE);
1233
1234   Double_t *tempArray = new Double_t[entries];
1235
1236   Double_t count = 0.5;
1237   Double_t  *vrun = new Double_t[entries];
1238   Int_t icount=0;
1239   //
1240   tempArray[index[0]] = count;
1241   vrun[0] = graph->GetX()[index[0]];
1242   for(Int_t i=1;i<entries;i++){
1243     if(graph->GetX()[index[i]]==graph->GetX()[index[i-1]])
1244       tempArray[index[i]] = count;
1245     else if(graph->GetX()[index[i]]!=graph->GetX()[index[i-1]]){
1246       count++;
1247       icount++;
1248       tempArray[index[i]] = count;
1249       vrun[icount]=graph->GetX()[index[i]];
1250     }
1251   }
1252
1253   const Int_t newNbins = int(count+0.5);
1254   Double_t *newBins = new Double_t[newNbins+1];
1255   for(Int_t i=0; i<=count+1;i++){
1256     newBins[i] = i;
1257   }
1258
1259   TGraph *graphNew = 0;
1260   if (tree->GetV3()) graphNew = new TGraphErrors(entries,tempArray,graph->GetY(),0,tree->GetV3());
1261   else
1262     graphNew = new TGraphErrors(entries,tempArray,graph->GetY(),0,0);
1263   graphNew->GetXaxis()->Set(newNbins,newBins);
1264
1265   Char_t xName[50];
1266   for(Int_t i=0;i<count;i++){
1267     snprintf(xName,50,"%d",Int_t(vrun[i]));
1268     graphNew->GetXaxis()->SetBinLabel(i+1,xName);
1269   }
1270   graphNew->GetHistogram()->SetTitle("");
1271   graphNew->SetMarkerStyle(mstyle);
1272   graphNew->SetMarkerColor(mcolor);
1273   if (msize>0) graphNew->SetMarkerSize(msize);
1274   delete [] tempArray;
1275   delete [] index;
1276   delete [] newBins;
1277   delete [] vrun;
1278   return graphNew;
1279 }
1280