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    }
 760
 761
 762    formulaStr.ReplaceAll("++", "~");
 763    TObjArray* formulaTokens = formulaStr.Tokenize("~");
 764    Int_t dim = formulaTokens->GetEntriesFast();
 765
 766    fitParam.ResizeTo(dim);
 767    covMatrix.ResizeTo(dim,dim);
 768
 769    TLinearFitter* fitter = new TLinearFitter(dim+1, Form("hyp%d",dim));
 770    fitter->StoreData(kTRUE);
 771    fitter->ClearPoints();
 772
 773    Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff",  stop-start, start);
 774    if (entries == -1) return new TString("An ERROR has occured during fitting!");
 775    Double_t **values = new Double_t*[dim+1] ;
 776    for (Int_t i=0; i<dim+1; i++) values[i]=NULL;
 777    //
 778    entries = tree->Draw(ferr.Data(), cutStr.Data(), "goff",  stop-start, start);
 779    if (entries == -1) {
 780      delete []values;
 781      return new TString("An ERROR has occured during fitting!");
 782    }
 783    Double_t *errors = new Double_t[entries];
 784    memcpy(errors,  tree->GetV1(), entries*sizeof(Double_t));
 785
 786    for (Int_t i = 0; i < dim + 1; i++){
 787       Int_t centries = 0;
 788       if (i < dim) centries = tree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start);
 789       else  centries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start);
 790
 791       if (entries != centries) {
 792         delete []errors;
 793         delete []values;
 794         return new TString("An ERROR has occured during fitting!");
 795       }
 796       values[i] = new Double_t[entries];
 797       memcpy(values[i],  tree->GetV1(), entries*sizeof(Double_t));
 798    }
 799
 800    // add points to the fitter
 801    for (Int_t i = 0; i < entries; i++){
 802       Double_t x[1000];
 803       for (Int_t j=0; j<dim;j++) x[j]=values[j][i];
 804       fitter->AddPoint(x, values[dim][i], errors[i]);
 805    }
 806
 807    fitter->Eval();
 808    if (frac>0.5 && frac<1){
 809      fitter->EvalRobust(frac);
 810    }else{
 811      if (fix0) {
 812        fitter->FixParameter(0,0);
 813        fitter->Eval();
 814      }
 815    }
 816    fitter->GetParameters(fitParam);
 817    fitter->GetCovarianceMatrix(covMatrix);
 818    chi2 = fitter->GetChisquare();
 819    npoints = entries;
 820    TString *preturnFormula = new TString(Form("( %f+",fitParam[0])), &returnFormula = *preturnFormula;
 821
 822    for (Int_t iparam = 0; iparam < dim; iparam++) {
 823      returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam+1]));
 824      if (iparam < dim-1) returnFormula.Append("+");
 825    }
 826    returnFormula.Append(" )");
 827
 828
 829    for (Int_t j=0; j<dim+1;j++) delete [] values[j];
 830
 831
 832    delete formulaTokens;
 833    delete fitter;
 834    delete[] values;
 835    delete[] errors;
 836    return preturnFormula;
 837 }
 838
 839 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){
 840    //
 841    // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts
 842    // returns chi2, fitParam and covMatrix
 843    // returns TString with fitted formula
 844    //
 845
 846    TString formulaStr(formula);
 847    TString drawStr(drawCommand);
 848    TString cutStr(cuts);
 849    TString ferr("1");
 850
 851    TString strVal(drawCommand);
 852    if (strVal.Contains(":")){
 853      TObjArray* valTokens = strVal.Tokenize(":");
 854      drawStr = valTokens->At(0)->GetName();
 855      ferr       = valTokens->At(1)->GetName();
 856    }
 857
 858
 859    formulaStr.ReplaceAll("++", "~");
 860    TObjArray* formulaTokens = formulaStr.Tokenize("~");
 861    Int_t dim = formulaTokens->GetEntriesFast();
 862
 863    fitParam.ResizeTo(dim);
 864    covMatrix.ResizeTo(dim,dim);
 865
 866    TLinearFitter* fitter = new TLinearFitter(dim+1, Form("hyp%d",dim));
 867    fitter->StoreData(kTRUE);
 868    fitter->ClearPoints();
 869
 870    Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff",  stop-start, start);
 871    if (entries == -1) return new TString("An ERROR has occured during fitting!");
 872    Double_t **values = new Double_t*[dim+1] ;
 873    for (Int_t i=0; i<dim+1; i++) values[i]=NULL;
 874    //
 875    entries = tree->Draw(ferr.Data(), cutStr.Data(), "goff",  stop-start, start);
 876    if (entries == -1) {
 877      delete [] values;
 878      return new TString("An ERROR has occured during fitting!");
 879    }
 880    Double_t *errors = new Double_t[entries];
 881    memcpy(errors,  tree->GetV1(), entries*sizeof(Double_t));
 882
 883    for (Int_t i = 0; i < dim + 1; i++){
 884       Int_t centries = 0;
 885       if (i < dim) centries = tree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start);
 886       else  centries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start);
 887
 888       if (entries != centries) {
 889         delete []errors;
 890         delete []values;
 891         return new TString("An ERROR has occured during fitting!");
 892       }
 893       values[i] = new Double_t[entries];
 894       memcpy(values[i],  tree->GetV1(), entries*sizeof(Double_t));
 895    }
 896
 897    // add points to the fitter
 898    for (Int_t i = 0; i < entries; i++){
 899       Double_t x[1000];
 900       for (Int_t j=0; j<dim;j++) x[j]=values[j][i];
 901       fitter->AddPoint(x, values[dim][i], errors[i]);
 902    }
 903    if (constrain>0){
 904      for (Int_t i = 0; i < dim; i++){
 905        Double_t x[1000];
 906        for (Int_t j=0; j<dim;j++) if (i!=j) x[j]=0;
 907        x[i]=1.;
 908        fitter->AddPoint(x, 0, constrain);
 909      }
 910    }
 911
 912
 913    fitter->Eval();
 914    if (frac>0.5 && frac<1){
 915      fitter->EvalRobust(frac);
 916    }
 917    fitter->GetParameters(fitParam);
 918    fitter->GetCovarianceMatrix(covMatrix);
 919    chi2 = fitter->GetChisquare();
 920    npoints = entries;
 921
 922    TString *preturnFormula = new TString(Form("( %f+",fitParam[0])), &returnFormula = *preturnFormula;
 923
 924    for (Int_t iparam = 0; iparam < dim; iparam++) {
 925      returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam+1]));
 926      if (iparam < dim-1) returnFormula.Append("+");
 927    }
 928    returnFormula.Append(" )");
 929
 930    for (Int_t j=0; j<dim+1;j++) delete [] values[j];
 931
 932
 933
 934    delete formulaTokens;
 935    delete fitter;
 936    delete[] values;
 937    delete[] errors;
 938    return preturnFormula;
 939 }
 940
 941
 942
 943 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){
 944    //
 945    // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts
 946    // returns chi2, fitParam and covMatrix
 947    // returns TString with fitted formula
 948    //
 949
 950    TString formulaStr(formula);
 951    TString drawStr(drawCommand);
 952    TString cutStr(cuts);
 953    TString ferr("1");
 954
 955    TString strVal(drawCommand);
 956    if (strVal.Contains(":")){
 957      TObjArray* valTokens = strVal.Tokenize(":");
 958      drawStr = valTokens->At(0)->GetName();
 959      ferr       = valTokens->At(1)->GetName();
 960    }
 961
 962
 963    formulaStr.ReplaceAll("++", "~");
 964    TObjArray* formulaTokens = formulaStr.Tokenize("~");
 965    Int_t dim = formulaTokens->GetEntriesFast();
 966
 967    fitParam.ResizeTo(dim);
 968    covMatrix.ResizeTo(dim,dim);
 969    TString fitString="x0";
 970    for (Int_t i=1; i<dim; i++) fitString+=Form("++x%d",i);
 971    TLinearFitter* fitter = new TLinearFitter(dim, fitString.Data());
 972    fitter->StoreData(kTRUE);
 973    fitter->ClearPoints();
 974
 975    Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff",  stop-start, start);
 976    if (entries == -1) return new TString("An ERROR has occured during fitting!");
 977    Double_t **values = new Double_t*[dim+1] ;
 978    for (Int_t i=0; i<dim+1; i++) values[i]=NULL;
 979    //
 980    entries = tree->Draw(ferr.Data(), cutStr.Data(), "goff",  stop-start, start);
 981    if (entries == -1) {
 982      delete []values;
 983      return new TString("An ERROR has occured during fitting!");
 984    }
 985    Double_t *errors = new Double_t[entries];
 986    memcpy(errors,  tree->GetV1(), entries*sizeof(Double_t));
 987
 988    for (Int_t i = 0; i < dim + 1; i++){
 989       Int_t centries = 0;
 990       if (i < dim) centries = tree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start);
 991       else  centries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start);
 992
 993       if (entries != centries) {
 994         delete []errors;
 995         delete []values;
 996         return new TString("An ERROR has occured during fitting!");
 997       }
 998       values[i] = new Double_t[entries];
 999       memcpy(values[i],  tree->GetV1(), entries*sizeof(Double_t));
1000    }
1001
1002    // add points to the fitter
1003    for (Int_t i = 0; i < entries; i++){
1004       Double_t x[1000];
1005       for (Int_t j=0; j<dim;j++) x[j]=values[j][i];
1006       fitter->AddPoint(x, values[dim][i], errors[i]);
1007    }
1008
1009    fitter->Eval();
1010    if (frac>0.5 && frac<1){
1011      fitter->EvalRobust(frac);
1012    }
1013    fitter->GetParameters(fitParam);
1014    fitter->GetCovarianceMatrix(covMatrix);
1015    chi2 = fitter->GetChisquare();
1016    npoints = entries;
1017
1018    TString *preturnFormula = new TString("("), &returnFormula = *preturnFormula;
1019
1020    for (Int_t iparam = 0; iparam < dim; iparam++) {
1021      returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam]));
1022      if (iparam < dim-1) returnFormula.Append("+");
1023    }
1024    returnFormula.Append(" )");
1025
1026
1027    for (Int_t j=0; j<dim+1;j++) delete [] values[j];
1028
1029    delete formulaTokens;
1030    delete fitter;
1031    delete[] values;
1032    delete[] errors;
1033    return preturnFormula;
1034 }
1035
1036
1037
1038
1039
1040 Int_t TStatToolkit::GetFitIndex(const TString fString, const TString subString){
1041   //
1042   // fitString - ++ separated list of fits
1043   // substring - ++ separated list of the requiered substrings
1044   //
1045   // return the last occurance of substring in fit string
1046   //
1047   TObjArray *arrFit = fString.Tokenize("++");
1048   TObjArray *arrSub = subString.Tokenize("++");
1049   Int_t index=-1;
1050   for (Int_t i=0; i<arrFit->GetEntries(); i++){
1051     Bool_t isOK=kTRUE;
1052     TString str =arrFit->At(i)->GetName();
1053     for (Int_t isub=0; isub<arrSub->GetEntries(); isub++){
1054       if (str.Contains(arrSub->At(isub)->GetName())==0) isOK=kFALSE;
1055     }
1056     if (isOK) index=i;
1057   }
1058   return index;
1059 }
1060
1061
1062 TString  TStatToolkit::FilterFit(const TString &input, const TString filter, TVectorD &param, TMatrixD & covar){
1063   //
1064   // Filter fit expression make sub-fit
1065   //
1066   TObjArray *array0= input.Tokenize("++");
1067   TObjArray *array1= filter.Tokenize("++");
1068   //TString *presult=new TString("(0");
1069   TString result="(0.0";
1070   for (Int_t i=0; i<array0->GetEntries(); i++){
1071     Bool_t isOK=kTRUE;
1072     TString str(array0->At(i)->GetName());
1073     for (Int_t j=0; j<array1->GetEntries(); j++){
1074       if (str.Contains(array1->At(j)->GetName())==0) isOK=kFALSE;
1075     }
1076     if (isOK) {
1077       result+="+"+str;
1078       result+=Form("*(%f)",param[i+1]);
1079       printf("%f\t%f\t%s\n",param[i+1], TMath::Sqrt(covar(i+1,i+1)),str.Data());
1080     }
1081   }
1082   result+="-0.)";
1083   return result;
1084 }
1085
1086 void TStatToolkit::Update1D(Double_t delta, Double_t sigma, Int_t s1, TMatrixD &vecXk, TMatrixD &covXk){
1087   //
1088   // Update parameters and covariance - with one measurement
1089   // Input:
1090   // vecXk - input vector - Updated in function
1091   // covXk - covariance matrix - Updated in function
1092   // delta, sigma, s1 - new measurement, rms of new measurement and the index of measurement
1093   const Int_t knMeas=1;
1094   Int_t knElem=vecXk.GetNrows();
1095
1096   TMatrixD mat1(knElem,knElem);            // update covariance matrix
1097   TMatrixD matHk(1,knElem);        // vector to mesurement
1098   TMatrixD vecYk(knMeas,1);        // Innovation or measurement residual
1099   TMatrixD matHkT(knElem,knMeas);  // helper matrix Hk transpose
1100   TMatrixD matSk(knMeas,knMeas);   // Innovation (or residual) covariance
1101   TMatrixD matKk(knElem,knMeas);   // Optimal Kalman gain
1102   TMatrixD covXk2(knElem,knElem);  // helper matrix
1103   TMatrixD covXk3(knElem,knElem);  // helper matrix
1104   TMatrixD vecZk(1,1);
1105   TMatrixD measR(1,1);
1106   vecZk(0,0)=delta;
1107   measR(0,0)=sigma*sigma;
1108   //
1109   // reset matHk
1110   for (Int_t iel=0;iel<knElem;iel++)
1111     for (Int_t ip=0;ip<knMeas;ip++) matHk(ip,iel)=0;
1112   //mat1
1113   for (Int_t iel=0;iel<knElem;iel++) {
1114     for (Int_t jel=0;jel<knElem;jel++) mat1(iel,jel)=0;
1115     mat1(iel,iel)=1;
1116   }
1117   //
1118   matHk(0, s1)=1;
1119   vecYk = vecZk-matHk*vecXk;               // Innovation or measurement residual
1120   matHkT=matHk.T(); matHk.T();
1121   matSk = (matHk*(covXk*matHkT))+measR;    // Innovation (or residual) covariance
1122   matSk.Invert();
1123   matKk = (covXk*matHkT)*matSk;            //  Optimal Kalman gain
1124   vecXk += matKk*vecYk;                    //  updated vector
1125   covXk2= (mat1-(matKk*matHk));
1126   covXk3 =  covXk2*covXk;
1127   covXk = covXk3;
1128   Int_t nrows=covXk3.GetNrows();
1129
1130   for (Int_t irow=0; irow<nrows; irow++)
1131     for (Int_t icol=0; icol<nrows; icol++){
1132       // rounding problems - make matrix again symteric
1133       covXk(irow,icol)=(covXk3(irow,icol)+covXk3(icol,irow))*0.5;
1134     }
1135 }
1136
1137
1138
1139 void   TStatToolkit::Constrain1D(const TString &input, const TString filter, TVectorD &param, TMatrixD & covar, Double_t mean, Double_t sigma){
1140   //
1141   // constrain linear fit
1142   // input  - string description of fit function
1143   // filter - string filter to select sub fits
1144   // param,covar - parameters and covariance matrix of the fit
1145   // mean,sigma  - new measurement uning which the fit is updated
1146   //
1147
1148   TObjArray *array0= input.Tokenize("++");
1149   TObjArray *array1= filter.Tokenize("++");
1150   TMatrixD paramM(param.GetNrows(),1);
1151   for (Int_t i=0; i<=array0->GetEntries(); i++){paramM(i,0)=param(i);}
1152
1153   if (filter.Length()==0){
1154     TStatToolkit::Update1D(mean, sigma, 0, paramM, covar);//
1155   }else{
1156     for (Int_t i=0; i<array0->GetEntries(); i++){
1157       Bool_t isOK=kTRUE;
1158       TString str(array0->At(i)->GetName());
1159       for (Int_t j=0; j<array1->GetEntries(); j++){
1160         if (str.Contains(array1->At(j)->GetName())==0) isOK=kFALSE;
1161       }
1162       if (isOK) {
1163         TStatToolkit::Update1D(mean, sigma, i+1, paramM, covar);//
1164       }
1165     }
1166   }
1167   for (Int_t i=0; i<=array0->GetEntries(); i++){
1168     param(i)=paramM(i,0);
1169   }
1170 }
1171
1172 TString  TStatToolkit::MakeFitString(const TString &input, const TVectorD &param, const TMatrixD & covar, Bool_t verbose){
1173   //
1174   //
1175   //
1176   TObjArray *array0= input.Tokenize("++");
1177   TString result=Form("(%f",param[0]);
1178   printf("%f\t%f\t\n", param[0], TMath::Sqrt(covar(0,0)));
1179   for (Int_t i=0; i<array0->GetEntries(); i++){
1180     TString str(array0->At(i)->GetName());
1181     result+="+"+str;
1182     result+=Form("*(%f)",param[i+1]);
1183     if (verbose) printf("%f\t%f\t%s\n", param[i+1], TMath::Sqrt(covar(i+1,i+1)),str.Data());
1184   }
1185   result+="-0.)";
1186   return result;
1187 }
1188
1189
1190 TGraph * TStatToolkit::MakeGraphSparse(TTree * tree, const char * expr, const char * cut, Int_t mstyle, Int_t mcolor){
1191   //
1192   // Make a sparse draw of the variables
1193   // Writen by Weilin.Yu
1194   const Int_t entries =  tree->Draw(expr,cut,"goff");
1195   //  TGraph * graph = (TGraph*)gPad->GetPrimitive("Graph"); // 2D
1196   TGraph * graph = 0;
1197   if (tree->GetV3()) graph = new TGraphErrors (entries, tree->GetV2(),tree->GetV1(),0,tree->GetV3());
1198   graph =  new TGraphErrors (entries, tree->GetV2(),tree->GetV1(),0,0);
1199   graph->SetMarkerStyle(mstyle);
1200   graph->SetMarkerColor(mcolor);
1201   //
1202   Int_t *index = new Int_t[entries*4];
1203   TMath::Sort(entries,graph->GetX(),index,kFALSE);
1204
1205   Double_t *tempArray = new Double_t[entries];
1206
1207   Double_t count = 0.5;
1208   Double_t  *vrun = new Double_t[entries];
1209   Int_t icount=0;
1210   //
1211   tempArray[index[0]] = count;
1212   vrun[0] = graph->GetX()[index[0]];
1213   for(Int_t i=1;i<entries;i++){
1214     if(graph->GetX()[index[i]]==graph->GetX()[index[i-1]])
1215       tempArray[index[i]] = count;
1216     else if(graph->GetX()[index[i]]!=graph->GetX()[index[i-1]]){
1217       count++;
1218       icount++;
1219       tempArray[index[i]] = count;
1220       vrun[icount]=graph->GetX()[index[i]];
1221     }
1222   }
1223
1224   const Int_t newNbins = int(count+0.5);
1225   Double_t *newBins = new Double_t[newNbins+1];
1226   for(Int_t i=0; i<=count+1;i++){
1227     newBins[i] = i;
1228   }
1229
1230   TGraph *graphNew = 0;
1231   if (tree->GetV3()) graphNew = new TGraphErrors(entries,tempArray,graph->GetY(),0,tree->GetV3());
1232   else
1233     graphNew = new TGraphErrors(entries,tempArray,graph->GetY(),0,0);
1234   graphNew->GetXaxis()->Set(newNbins,newBins);
1235
1236   Char_t xName[50];
1237   for(Int_t i=0;i<count;i++){
1238     snprintf(xName,50,"%d",Int_t(vrun[i]));
1239     graphNew->GetXaxis()->SetBinLabel(i+1,xName);
1240   }
1241   graphNew->GetHistogram()->SetTitle("");
1242   graphNew->SetMarkerStyle(mstyle);
1243   graphNew->SetMarkerColor(mcolor);
1244   delete [] tempArray;
1245   delete [] index;
1246   delete [] newBins;
1247   delete [] vrun;
1248   return graphNew;
1249 }
1250