]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - STAT/TStatToolkit.cxx
fix coverity
[u/mrichter/AliRoot.git] / STAT / TStatToolkit.cxx
index a5981a759a438cfd60ccfa030a60007eeddd69e7..1fb67062b0e109a2cd3e0134c30d5338e144055c 100644 (file)
@@ -30,6 +30,8 @@
 #include "TChain.h"
 #include "TObjString.h"
 #include "TLinearFitter.h"
+#include "TGraph2D.h"
+#include "TGraph.h"
 
 //
 // includes neccessary for test functions
@@ -184,7 +186,8 @@ Int_t TStatToolkit::Freq(Int_t n, const Int_t *inlist
 
   Int_t * sindexS = new Int_t[n];     // temp array for sorting
   Int_t * sindexF = new Int_t[2*n];   
-  for (Int_t i=0;i<n;i++) sindexF[i]=0;
+  for (Int_t i=0;i<n;i++) sindexS[i]=0;
+  for (Int_t i=0;i<2*n;i++) sindexF[i]=0;
   //
   TMath::Sort(n,inlist, sindexS, down);  
   Int_t last      = inlist[sindexS[0]];
@@ -219,7 +222,7 @@ Int_t TStatToolkit::Freq(Int_t n, const Int_t *inlist
 }
 
 //___TStatToolkit__________________________________________________________________________
-void TStatToolkit::TruncatedMean(TH1F * his, TVectorD *param, Float_t down, Float_t up, Bool_t verbose){
+void TStatToolkit::TruncatedMean(const TH1 * his, TVectorD *param, Float_t down, Float_t up, Bool_t verbose){
   //
   //
   //
@@ -279,7 +282,7 @@ void TStatToolkit::LTM(TH1F * his, TVectorD *param , Float_t fraction,  Bool_t v
   }
 }
 
-Double_t  TStatToolkit::FitGaus(TH1F* his, TVectorD *param, TMatrixD *matrix, Float_t xmin, Float_t xmax, Bool_t verbose){
+Double_t  TStatToolkit::FitGaus(TH1* his, TVectorD *param, TMatrixD */*matrix*/, Float_t xmin, Float_t xmax, Bool_t verbose){
   //
   //  Fit histogram with gaussian function
   //  
@@ -357,7 +360,7 @@ Double_t  TStatToolkit::FitGaus(TH1F* his, TVectorD *param, TMatrixD *matrix, Fl
   Double_t chi2 = fitter.GetChisquare()/Float_t(npoints);
   //fitter.GetParameters();
   if (!param)  param  = new TVectorD(3);
-  if (!matrix) matrix = new TMatrixD(3,3);
+  // if (!matrix) matrix = new TMatrixD(3,3); // Covariance matrix to be implemented
   (*param)[1] = par[1]/(-2.*par[2]);
   (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
   (*param)[0] = TMath::Exp(par[0]+ par[1]* (*param)[1] +  par[2]*(*param)[1]*(*param)[1]);
@@ -375,7 +378,7 @@ Double_t  TStatToolkit::FitGaus(TH1F* his, TVectorD *param, TMatrixD *matrix, Fl
   return chi2;
 }
 
-Double_t  TStatToolkit::FitGaus(Float_t *arr, Int_t nBins, Float_t xMin, Float_t xMax, TVectorD *param, TMatrixD *matrix, Bool_t verbose){
+Double_t  TStatToolkit::FitGaus(Float_t *arr, Int_t nBins, Float_t xMin, Float_t xMax, TVectorD *param, TMatrixD */*matrix*/, Bool_t verbose){
   //
   //  Fit histogram with gaussian function
   //  
@@ -403,7 +406,7 @@ Double_t  TStatToolkit::FitGaus(Float_t *arr, Int_t nBins, Float_t xMin, Float_t
   fitter.ClearPoints();
   TVectorD  par(3);
   TVectorD  sigma(3);
-  TMatrixD A(3,3);
+  TMatrixD matA(3,3);
   TMatrixD b(3,1);
   Float_t rms = TMath::RMS(nBins,arr);
   Float_t max = TMath::MaxElement(nBins,arr);
@@ -438,9 +441,9 @@ Double_t  TStatToolkit::FitGaus(Float_t *arr, Int_t nBins, Float_t xMin, Float_t
       Float_t val = TMath::Log(Float_t(entriesI));
       fitter.AddPoint(&xcenter,val,error);
       if (npoints<3){
-         A(npoints,0)=1;
-         A(npoints,1)=xcenter;
-         A(npoints,2)=xcenter*xcenter;
+         matA(npoints,0)=1;
+         matA(npoints,1)=xcenter;
+         matA(npoints,2)=xcenter*xcenter;
          b(npoints,0)=val;
          meanCOG+=xcenter*entriesI;
          rms2COG +=xcenter*entriesI*xcenter;
@@ -455,9 +458,9 @@ Double_t  TStatToolkit::FitGaus(Float_t *arr, Int_t nBins, Float_t xMin, Float_t
   if (npoints>=3){
       if ( npoints == 3 ){
          //analytic calculation of the parameters for three points
-         A.Invert();
+         matA.Invert();
          TMatrixD res(1,3);
-         res.Mult(A,b);
+         res.Mult(matA,b);
          par[0]=res(0,0);
          par[1]=res(0,1);
          par[2]=res(0,2);
@@ -473,7 +476,7 @@ Double_t  TStatToolkit::FitGaus(Float_t *arr, Int_t nBins, Float_t xMin, Float_t
       if (TMath::Abs(par[2])<kTol) return -4;
 
       if (!param)  param  = new TVectorD(3);
-      if (!matrix) matrix = new TMatrixD(3,3);  // !!!!might be a memory leek. use dummy matrix pointer to call this function!
+      //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
 
       (*param)[1] = par[1]/(-2.*par[2]);
       (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
@@ -515,7 +518,7 @@ Double_t  TStatToolkit::FitGaus(Float_t *arr, Int_t nBins, Float_t xMin, Float_t
 }
 
 
-Float_t TStatToolkit::GetCOG(Short_t *arr, Int_t nBins, Float_t xMin, Float_t xMax, Float_t *rms, Float_t *sum)
+Float_t TStatToolkit::GetCOG(const Short_t *arr, Int_t nBins, Float_t xMin, Float_t xMax, Float_t *rms, Float_t *sum)
 {
     //
     //  calculate center of gravity rms and sum for array 'arr' with nBins an a x range xMin to xMax
@@ -666,7 +669,7 @@ TGraph2D * TStatToolkit::MakeStat2D(TH3 * his, Int_t delta0, Int_t delta1, Int_t
     for (Int_t iy=0; iy<nbiny;iy++){
       Float_t xcenter = xaxis->GetBinCenter(ix); 
       Float_t ycenter = yaxis->GetBinCenter(iy); 
-      sprintf(name,"%s_%d_%d",his->GetName(), ix,iy);
+      snprintf(name,1000,"%s_%d_%d",his->GetName(), ix,iy);
       TH1 *projection = his->ProjectionZ(name,ix-delta0,ix+delta0,iy-delta1,iy+delta1);
       Float_t stat= 0;
       if (type==0) stat = projection->GetMean();
@@ -708,7 +711,7 @@ TGraph * TStatToolkit::MakeStat1D(TH3 * his, Int_t delta1, Int_t type){
   for (Int_t ix=0; ix<nbinx;ix++){
     Float_t xcenter = xaxis->GetBinCenter(ix); 
     //    Float_t ycenter = yaxis->GetBinCenter(iy); 
-    sprintf(name,"%s_%d",his->GetName(), ix);
+    snprintf(name,1000,"%s_%d",his->GetName(), ix);
     TH1 *projection = his->ProjectionZ(name,ix-delta1,ix+delta1,0,nbiny);
     Float_t stat= 0;
     if (type==0) stat = projection->GetMean();
@@ -735,7 +738,7 @@ TGraph * TStatToolkit::MakeStat1D(TH3 * his, Int_t delta1, Int_t type){
 
 
 
-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){
+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){
    //
    // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts
    // returns chi2, fitParam and covMatrix
@@ -771,7 +774,10 @@ TString* TStatToolkit::FitPlane(TTree *tree, const char* drawCommand, const char
    Double_t **values = new Double_t*[dim+1] ; 
    //
    entries = tree->Draw(ferr.Data(), cutStr.Data(), "goff",  stop-start, start);
-   if (entries == -1) return new TString("An ERROR has occured during fitting!");
+   if (entries == -1) {
+     delete []values;
+     return new TString("An ERROR has occured during fitting!");
+   }
    Double_t *errors = new Double_t[entries];
    memcpy(errors,  tree->GetV1(), entries*sizeof(Double_t));
    
@@ -780,7 +786,11 @@ TString* TStatToolkit::FitPlane(TTree *tree, const char* drawCommand, const char
       if (i < dim) centries = tree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start);
       else  centries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start);
       
-      if (entries != centries) return new TString("An ERROR has occured during fitting!");
+      if (entries != centries) {
+       delete []errors;
+       delete []values;
+       return new TString("An ERROR has occured during fitting!");
+      }
       values[i] = new Double_t[entries];
       memcpy(values[i],  tree->GetV1(), entries*sizeof(Double_t)); 
    }
@@ -795,19 +805,116 @@ TString* TStatToolkit::FitPlane(TTree *tree, const char* drawCommand, const char
    fitter->Eval();
    if (frac>0.5 && frac<1){
      fitter->EvalRobust(frac);
+   }else{
+     if (fix0) {
+       fitter->FixParameter(0,0);
+       fitter->Eval();     
+     }
    }
    fitter->GetParameters(fitParam);
    fitter->GetCovarianceMatrix(covMatrix);
    chi2 = fitter->GetChisquare();
-   chi2 = chi2;
-   npoints = entries;
-//    TString *preturnFormula = new TString(Form("%f*(",fitParam[0])), &returnFormula = *preturnFormula; 
+   npoints = entries;   
+   TString *preturnFormula = new TString(Form("( %f+",fitParam[0])), &returnFormula = *preturnFormula; 
+   
+   for (Int_t iparam = 0; iparam < dim; iparam++) {
+     returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam+1]));
+     if (iparam < dim-1) returnFormula.Append("+");
+   }
+   returnFormula.Append(" )");
+   
+   
+   for (Int_t j=0; j<dim+1;j++) delete [] values[j];
+
+
+   delete formulaTokens;
+   delete fitter;
+   delete[] values;
+   delete[] errors;
+   return preturnFormula;
+}
+
+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){
+   //
+   // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts
+   // returns chi2, fitParam and covMatrix
+   // returns TString with fitted formula
+   //
+
+   TString formulaStr(formula); 
+   TString drawStr(drawCommand);
+   TString cutStr(cuts);
+   TString ferr("1");
+
+   TString strVal(drawCommand);
+   if (strVal.Contains(":")){
+     TObjArray* valTokens = strVal.Tokenize(":");
+     drawStr = valTokens->At(0)->GetName();
+     ferr       = valTokens->At(1)->GetName();     
+   }
+
+      
+   formulaStr.ReplaceAll("++", "~");
+   TObjArray* formulaTokens = formulaStr.Tokenize("~"); 
+   Int_t dim = formulaTokens->GetEntriesFast();
+   
+   fitParam.ResizeTo(dim);
+   covMatrix.ResizeTo(dim,dim);
+   
+   TLinearFitter* fitter = new TLinearFitter(dim+1, Form("hyp%d",dim));
+   fitter->StoreData(kTRUE);   
+   fitter->ClearPoints();
+   
+   Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff",  stop-start, start);
+   if (entries == -1) return new TString("An ERROR has occured during fitting!");
+   Double_t **values = new Double_t*[dim+1] ; 
+   //
+   entries = tree->Draw(ferr.Data(), cutStr.Data(), "goff",  stop-start, start);
+   if (entries == -1) {
+     delete [] values;
+     return new TString("An ERROR has occured during fitting!");
+   }
+   Double_t *errors = new Double_t[entries];
+   memcpy(errors,  tree->GetV1(), entries*sizeof(Double_t));
    
-//    for (Int_t iparam = 0; iparam < dim; iparam++) {
-//      returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam+1]/fitParam[0]));
-//      if (iparam < dim-1) returnFormula.Append("+");
-//    }
-//    returnFormula.Append(" )");
+   for (Int_t i = 0; i < dim + 1; i++){
+      Int_t centries = 0;
+      if (i < dim) centries = tree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start);
+      else  centries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start);
+      
+      if (entries != centries) {
+       delete []errors;
+       delete []values;
+       return new TString("An ERROR has occured during fitting!");
+      }
+      values[i] = new Double_t[entries];
+      memcpy(values[i],  tree->GetV1(), entries*sizeof(Double_t)); 
+   }
+   
+   // add points to the fitter
+   for (Int_t i = 0; i < entries; i++){
+      Double_t x[1000];
+      for (Int_t j=0; j<dim;j++) x[j]=values[j][i];
+      fitter->AddPoint(x, values[dim][i], errors[i]);
+   }
+   if (constrain>0){
+     for (Int_t i = 0; i < dim; i++){
+       Double_t x[1000];
+       for (Int_t j=0; j<dim;j++) if (i!=j) x[j]=0;
+       x[i]=1.;
+       fitter->AddPoint(x, 0, constrain);
+     }
+   }
+
+
+   fitter->Eval();
+   if (frac>0.5 && frac<1){
+     fitter->EvalRobust(frac);   
+   }
+   fitter->GetParameters(fitParam);
+   fitter->GetCovarianceMatrix(covMatrix);
+   chi2 = fitter->GetChisquare();
+   npoints = entries;
    
    TString *preturnFormula = new TString(Form("( %f+",fitParam[0])), &returnFormula = *preturnFormula; 
    
@@ -817,11 +924,315 @@ TString* TStatToolkit::FitPlane(TTree *tree, const char* drawCommand, const char
    }
    returnFormula.Append(" )");
    
+   for (Int_t j=0; j<dim+1;j++) delete [] values[j];
    
 
 
    delete formulaTokens;
    delete fitter;
    delete[] values;
+   delete[] errors;
+   return preturnFormula;
+}
+
+
+
+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){
+   //
+   // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts
+   // returns chi2, fitParam and covMatrix
+   // returns TString with fitted formula
+   //
+
+   TString formulaStr(formula); 
+   TString drawStr(drawCommand);
+   TString cutStr(cuts);
+   TString ferr("1");
+
+   TString strVal(drawCommand);
+   if (strVal.Contains(":")){
+     TObjArray* valTokens = strVal.Tokenize(":");
+     drawStr = valTokens->At(0)->GetName();
+     ferr       = valTokens->At(1)->GetName();     
+   }
+
+      
+   formulaStr.ReplaceAll("++", "~");
+   TObjArray* formulaTokens = formulaStr.Tokenize("~"); 
+   Int_t dim = formulaTokens->GetEntriesFast();
+   
+   fitParam.ResizeTo(dim);
+   covMatrix.ResizeTo(dim,dim);
+   TString fitString="x0";
+   for (Int_t i=1; i<dim; i++) fitString+=Form("++x%d",i);     
+   TLinearFitter* fitter = new TLinearFitter(dim, fitString.Data());
+   fitter->StoreData(kTRUE);   
+   fitter->ClearPoints();
+   
+   Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff",  stop-start, start);
+   if (entries == -1) return new TString("An ERROR has occured during fitting!");
+   Double_t **values = new Double_t*[dim+1] ; 
+   //
+   entries = tree->Draw(ferr.Data(), cutStr.Data(), "goff",  stop-start, start);
+   if (entries == -1) {
+     delete []values;
+     return new TString("An ERROR has occured during fitting!");
+   }
+   Double_t *errors = new Double_t[entries];
+   memcpy(errors,  tree->GetV1(), entries*sizeof(Double_t));
+   
+   for (Int_t i = 0; i < dim + 1; i++){
+      Int_t centries = 0;
+      if (i < dim) centries = tree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start);
+      else  centries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start);
+      
+      if (entries != centries) {
+       delete []errors;
+       delete []values;
+       return new TString("An ERROR has occured during fitting!");
+      }
+      values[i] = new Double_t[entries];
+      memcpy(values[i],  tree->GetV1(), entries*sizeof(Double_t)); 
+   }
+   
+   // add points to the fitter
+   for (Int_t i = 0; i < entries; i++){
+      Double_t x[1000];
+      for (Int_t j=0; j<dim;j++) x[j]=values[j][i];
+      fitter->AddPoint(x, values[dim][i], errors[i]);
+   }
+
+   fitter->Eval();
+   if (frac>0.5 && frac<1){
+     fitter->EvalRobust(frac);
+   }
+   fitter->GetParameters(fitParam);
+   fitter->GetCovarianceMatrix(covMatrix);
+   chi2 = fitter->GetChisquare();
+   npoints = entries;
+   
+   TString *preturnFormula = new TString("("), &returnFormula = *preturnFormula; 
+   
+   for (Int_t iparam = 0; iparam < dim; iparam++) {
+     returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam]));
+     if (iparam < dim-1) returnFormula.Append("+");
+   }
+   returnFormula.Append(" )");
+   
+   
+   for (Int_t j=0; j<dim+1;j++) delete [] values[j];
+   
+   delete formulaTokens;
+   delete fitter;
+   delete[] values;
+   delete[] errors;
    return preturnFormula;
 }
+
+
+
+
+
+Int_t TStatToolkit::GetFitIndex(const TString fString, const TString subString){
+  //
+  // fitString - ++ separated list of fits
+  // substring - ++ separated list of the requiered substrings
+  //
+  // return the last occurance of substring in fit string
+  // 
+  TObjArray *arrFit = fString.Tokenize("++");
+  TObjArray *arrSub = subString.Tokenize("++");
+  Int_t index=-1;
+  for (Int_t i=0; i<arrFit->GetEntries(); i++){
+    Bool_t isOK=kTRUE;
+    TString str =arrFit->At(i)->GetName();
+    for (Int_t isub=0; isub<arrSub->GetEntries(); isub++){
+      if (str.Contains(arrSub->At(isub)->GetName())==0) isOK=kFALSE;
+    }
+    if (isOK) index=i;
+  }
+  return index;
+}
+
+
+TString  TStatToolkit::FilterFit(const TString &input, const TString filter, TVectorD &param, TMatrixD & covar){
+  //
+  // Filter fit expression make sub-fit
+  //
+  TObjArray *array0= input.Tokenize("++");
+  TObjArray *array1= filter.Tokenize("++");
+  //TString *presult=new TString("(0");
+  TString result="(0.0";
+  for (Int_t i=0; i<array0->GetEntries(); i++){
+    Bool_t isOK=kTRUE;
+    TString str(array0->At(i)->GetName());
+    for (Int_t j=0; j<array1->GetEntries(); j++){
+      if (str.Contains(array1->At(j)->GetName())==0) isOK=kFALSE;      
+    }
+    if (isOK) {
+      result+="+"+str;
+      result+=Form("*(%f)",param[i+1]);
+      printf("%f\t%f\t%s\n",param[i+1], TMath::Sqrt(covar(i+1,i+1)),str.Data());    
+    }
+  }
+  result+="-0.)";
+  return result;
+}
+
+void TStatToolkit::Update1D(Double_t delta, Double_t sigma, Int_t s1, TMatrixD &vecXk, TMatrixD &covXk){
+  //
+  // Update parameters and covariance - with one measurement
+  // Input:
+  // vecXk - input vector - Updated in function 
+  // covXk - covariance matrix - Updated in function
+  // delta, sigma, s1 - new measurement, rms of new measurement and the index of measurement
+  const Int_t knMeas=1;
+  Int_t knElem=vecXk.GetNrows();
+  TMatrixD mat1(knElem,knElem);            // update covariance matrix
+  TMatrixD matHk(1,knElem);        // vector to mesurement
+  TMatrixD vecYk(knMeas,1);        // Innovation or measurement residual
+  TMatrixD matHkT(knElem,knMeas);  // helper matrix Hk transpose
+  TMatrixD matSk(knMeas,knMeas);   // Innovation (or residual) covariance
+  TMatrixD matKk(knElem,knMeas);   // Optimal Kalman gain
+  TMatrixD covXk2(knElem,knElem);  // helper matrix
+  TMatrixD covXk3(knElem,knElem);  // helper matrix
+  TMatrixD vecZk(1,1);
+  TMatrixD measR(1,1);
+  vecZk(0,0)=delta;
+  measR(0,0)=sigma*sigma;
+  //
+  // reset matHk
+  for (Int_t iel=0;iel<knElem;iel++) 
+    for (Int_t ip=0;ip<knMeas;ip++) matHk(ip,iel)=0; 
+  //mat1
+  for (Int_t iel=0;iel<knElem;iel++) {
+    for (Int_t jel=0;jel<knElem;jel++) mat1(iel,jel)=0;
+    mat1(iel,iel)=1;
+  }
+  //
+  matHk(0, s1)=1;
+  vecYk = vecZk-matHk*vecXk;               // Innovation or measurement residual
+  matHkT=matHk.T(); matHk.T();
+  matSk = (matHk*(covXk*matHkT))+measR;    // Innovation (or residual) covariance
+  matSk.Invert();
+  matKk = (covXk*matHkT)*matSk;            //  Optimal Kalman gain
+  vecXk += matKk*vecYk;                    //  updated vector 
+  covXk2= (mat1-(matKk*matHk));
+  covXk3 =  covXk2*covXk;          
+  covXk = covXk3;  
+  Int_t nrows=covXk3.GetNrows();
+  
+  for (Int_t irow=0; irow<nrows; irow++)
+    for (Int_t icol=0; icol<nrows; icol++){
+      // rounding problems - make matrix again symteric
+      covXk(irow,icol)=(covXk3(irow,icol)+covXk3(icol,irow))*0.5; 
+    }
+}
+
+
+
+void   TStatToolkit::Constrain1D(const TString &input, const TString filter, TVectorD &param, TMatrixD & covar, Double_t mean, Double_t sigma){
+  //
+  // constrain linear fit
+  // input  - string description of fit function
+  // filter - string filter to select sub fits
+  // param,covar - parameters and covariance matrix of the fit
+  // mean,sigma  - new measurement uning which the fit is updated
+  //
+  
+  TObjArray *array0= input.Tokenize("++");
+  TObjArray *array1= filter.Tokenize("++");
+  TMatrixD paramM(param.GetNrows(),1);
+  for (Int_t i=0; i<=array0->GetEntries(); i++){paramM(i,0)=param(i);}
+  
+  if (filter.Length()==0){
+    TStatToolkit::Update1D(mean, sigma, 0, paramM, covar);//
+  }else{  
+    for (Int_t i=0; i<array0->GetEntries(); i++){
+      Bool_t isOK=kTRUE;
+      TString str(array0->At(i)->GetName());
+      for (Int_t j=0; j<array1->GetEntries(); j++){
+       if (str.Contains(array1->At(j)->GetName())==0) isOK=kFALSE;      
+      }
+      if (isOK) {
+       TStatToolkit::Update1D(mean, sigma, i+1, paramM, covar);//
+      }
+    }
+  }
+  for (Int_t i=0; i<=array0->GetEntries(); i++){
+    param(i)=paramM(i,0);
+  }
+}
+
+TString  TStatToolkit::MakeFitString(const TString &input, const TVectorD &param, const TMatrixD & covar, Bool_t verbose){
+  //
+  //
+  //
+  TObjArray *array0= input.Tokenize("++");
+  TString result=Form("(%f",param[0]);
+  printf("%f\t%f\t\n", param[0], TMath::Sqrt(covar(0,0))); 
+  for (Int_t i=0; i<array0->GetEntries(); i++){
+    TString str(array0->At(i)->GetName());
+    result+="+"+str;
+    result+=Form("*(%f)",param[i+1]);
+    if (verbose) printf("%f\t%f\t%s\n", param[i+1], TMath::Sqrt(covar(i+1,i+1)),str.Data());    
+  }
+  result+="-0.)";
+  return result;
+}
+
+
+TGraph * TStatToolkit::MakeGraphSparse(TTree * tree, const char * expr, const char * cut){
+  //
+  // Make a sparse draw of the variables
+  //
+  const Int_t entries =  tree->Draw(expr,cut,"goff");
+  //  TGraph * graph = (TGraph*)gPad->GetPrimitive("Graph"); // 2D
+  TGraph * graph = new TGraph (entries, tree->GetV2(),tree->GetV1());
+  //
+  Int_t *index = new Int_t[entries];
+  TMath::Sort(entries,graph->GetX(),index,kFALSE);
+  
+  Double_t *tempArray = new Double_t[entries];
+
+  Double_t count = 0.5;
+  Double_t  *vrun = new Double_t[entries];
+  Int_t icount=0;
+  //
+  tempArray[index[0]] = count;
+  vrun[0] = graph->GetX()[index[0]];
+  for(Int_t i=1;i<entries;i++){
+    if(graph->GetX()[index[i]]==graph->GetX()[index[i-1]])
+      tempArray[index[i]] = count; 
+    else if(graph->GetX()[index[i]]!=graph->GetX()[index[i-1]]){
+      count++;
+      icount++;
+      tempArray[index[i]] = count;
+      vrun[icount]=graph->GetX()[index[i]];
+    }
+  }
+  
+  const Int_t newNbins = int(count+0.5);
+  Double_t *newBins = new Double_t[newNbins+1];
+  for(Int_t i=0; i<=count+1;i++){
+    newBins[i] = i;
+  }
+  
+  TGraph *graphNew = new TGraph(entries,tempArray,graph->GetY());
+  graphNew->GetXaxis()->Set(newNbins,newBins);
+  
+  Char_t xName[50];
+  for(Int_t i=0;i<count;i++){
+    snprintf(xName,50,"%d",Int_t(vrun[i]));
+    graphNew->GetXaxis()->SetBinLabel(i+1,xName);
+  }
+  graphNew->GetHistogram()->SetTitle("");
+  
+  delete [] tempArray;
+  delete [] index;
+  delete [] newBins;
+  delete [] vrun;
+  return graphNew;
+}
+