X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STAT%2FTStatToolkit.cxx;h=1fb67062b0e109a2cd3e0134c30d5338e144055c;hb=f9d25f75bb6a84264f5594c513e5b3b108ff8beb;hp=24493eaf58a0bd813620e09685523b0bb2d6fdbd;hpb=7c9cf6e41f59dcf1deb5503952f74ac54d9ac47c;p=u%2Fmrichter%2FAliRoot.git

diff --git a/STAT/TStatToolkit.cxx b/STAT/TStatToolkit.cxx
index 24493eaf58a..1fb67062b0e 100644
--- a/STAT/TStatToolkit.cxx
+++ b/STAT/TStatToolkit.cxx
@@ -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
   //  
@@ -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);
@@ -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
@@ -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)); 
    }
@@ -804,15 +814,7 @@ TString* TStatToolkit::FitPlane(TTree *tree, const char* drawCommand, const char
    fitter->GetParameters(fitParam);
    fitter->GetCovarianceMatrix(covMatrix);
    chi2 = fitter->GetChisquare();
-   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]/fitParam[0]));
-//      if (iparam < dim-1) returnFormula.Append("+");
-//    }
-//    returnFormula.Append(" )");
-   
+   npoints = entries;   
    TString *preturnFormula = new TString(Form("( %f+",fitParam[0])), &returnFormula = *preturnFormula; 
    
    for (Int_t iparam = 0; iparam < dim; iparam++) {
@@ -822,11 +824,13 @@ 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;
 }
 
@@ -866,7 +870,10 @@ TString* TStatToolkit::FitPlaneConstrain(TTree *tree, const char* drawCommand, c
    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));
    
@@ -875,7 +882,11 @@ TString* TStatToolkit::FitPlaneConstrain(TTree *tree, const char* drawCommand, c
       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)); 
    }
@@ -904,13 +915,6 @@ TString* TStatToolkit::FitPlaneConstrain(TTree *tree, const char* drawCommand, c
    fitter->GetCovarianceMatrix(covMatrix);
    chi2 = fitter->GetChisquare();
    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]/fitParam[0]));
-//      if (iparam < dim-1) returnFormula.Append("+");
-//    }
-//    returnFormula.Append(" )");
    
    TString *preturnFormula = new TString(Form("( %f+",fitParam[0])), &returnFormula = *preturnFormula; 
    
@@ -920,12 +924,14 @@ TString* TStatToolkit::FitPlaneConstrain(TTree *tree, const char* drawCommand, c
    }
    returnFormula.Append(" )");
    
+   for (Int_t j=0; j<dim+1;j++) delete [] values[j];
    
 
 
    delete formulaTokens;
    delete fitter;
    delete[] values;
+   delete[] errors;
    return preturnFormula;
 }
 
@@ -968,7 +974,10 @@ TString* TStatToolkit::FitPlaneFixed(TTree *tree, const char* drawCommand, const
    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));
    
@@ -977,7 +986,11 @@ TString* TStatToolkit::FitPlaneFixed(TTree *tree, const char* drawCommand, const
       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)); 
    }
@@ -997,13 +1010,6 @@ TString* TStatToolkit::FitPlaneFixed(TTree *tree, const char* drawCommand, const
    fitter->GetCovarianceMatrix(covMatrix);
    chi2 = fitter->GetChisquare();
    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]/fitParam[0]));
-//      if (iparam < dim-1) returnFormula.Append("+");
-//    }
-//    returnFormula.Append(" )");
    
    TString *preturnFormula = new TString("("), &returnFormula = *preturnFormula; 
    
@@ -1014,11 +1020,12 @@ TString* TStatToolkit::FitPlaneFixed(TTree *tree, const char* drawCommand, const
    returnFormula.Append(" )");
    
    
-   
+   for (Int_t j=0; j<dim+1;j++) delete [] values[j];
    
    delete formulaTokens;
    delete fitter;
    delete[] values;
+   delete[] errors;
    return preturnFormula;
 }
 
@@ -1026,7 +1033,7 @@ TString* TStatToolkit::FitPlaneFixed(TTree *tree, const char* drawCommand, const
 
 
 
-Int_t TStatToolkit::GetFitIndex(TString fString, TString subString){
+Int_t TStatToolkit::GetFitIndex(const TString fString, const TString subString){
   //
   // fitString - ++ separated list of fits
   // substring - ++ separated list of the requiered substrings
@@ -1048,7 +1055,7 @@ Int_t TStatToolkit::GetFitIndex(TString fString, TString subString){
 }
 
 
-TString  TStatToolkit::FilterFit(TString &input, TString filter, TVectorD &param, TMatrixD & covar){
+TString  TStatToolkit::FilterFit(const TString &input, const TString filter, TVectorD &param, TMatrixD & covar){
   //
   // Filter fit expression make sub-fit
   //
@@ -1125,7 +1132,7 @@ void TStatToolkit::Update1D(Double_t delta, Double_t sigma, Int_t s1, TMatrixD &
 
 
 
-void   TStatToolkit::Constrain1D(TString &input, TString filter, TVectorD &param, TMatrixD & covar, Double_t mean, Double_t sigma){
+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
@@ -1133,19 +1140,24 @@ void   TStatToolkit::Constrain1D(TString &input, TString filter, TVectorD &param
   // 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);}
   
-  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);//
+  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++){
@@ -1153,18 +1165,74 @@ void   TStatToolkit::Constrain1D(TString &input, TString filter, TVectorD &param
   }
 }
 
-TString  TStatToolkit::MakeFitString(TString &input, TVectorD &param, TMatrixD & covar){
+TString  TStatToolkit::MakeFitString(const TString &input, const TVectorD &param, const TMatrixD & covar, Bool_t verbose){
   //
   //
   //
   TObjArray *array0= input.Tokenize("++");
-  TString result="(0.0";
+  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]);
-    printf("%f\t%f\t%s\n", param[i+1], TMath::Sqrt(covar(i+1,i+1)),str.Data());    
+    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;
+}
+