]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - STAT/TStatToolkit.cxx
coverty
[u/mrichter/AliRoot.git] / STAT / TStatToolkit.cxx
index 36a9d08688eab3e731497c8d7be8101f8b168498..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
@@ -220,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){
   //
   //
   //
@@ -280,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
   //  
@@ -404,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);
@@ -439,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;
@@ -456,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);
@@ -516,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
@@ -1031,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
@@ -1053,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
   //
@@ -1130,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
@@ -1138,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++){
@@ -1158,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;
+}
+