]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - STAT/TStatToolkit.cxx
Updated treatment of TOF PID in QA task (Francesco+Pietro)
[u/mrichter/AliRoot.git] / STAT / TStatToolkit.cxx
index 3382cb8605f6898c56b23017a058e79a108da1f1..8ea3c3c6c53089799e425b1b4deefddd3de40c20 100644 (file)
@@ -30,6 +30,9 @@
 #include "TChain.h"
 #include "TObjString.h"
 #include "TLinearFitter.h"
+#include "TGraph2D.h"
+#include "TGraph.h"
+#include "TGraphErrors.h"
 
 //
 // includes neccessary for test functions
@@ -95,7 +98,7 @@ void TStatToolkit::EvaluateUni(Int_t nvectors, Double_t *data, Double_t &mean
   }
   
   Double_t norm = 1./Double_t(hh);
-  Double_t norm2 = 1./Double_t(hh-1);
+  Double_t norm2 = (hh-1)>0 ? 1./Double_t(hh-1):1;
   for (Int_t i=hh; i<nvectors; i++){
     Double_t cmean  = sumx*norm;
     Double_t csigma = (sumx2 - hh*cmean*cmean)*norm2;
@@ -220,7 +223,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 +283,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 +407,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 +442,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 +459,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 +519,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
@@ -769,7 +772,8 @@ TString* TStatToolkit::FitPlane(TTree *tree, const char* drawCommand, const char
    
    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] ; 
+   Double_t **values = new Double_t*[dim+1] ;
+   for (Int_t i=0; i<dim+1; i++) values[i]=NULL; 
    //
    entries = tree->Draw(ferr.Data(), cutStr.Data(), "goff",  stop-start, start);
    if (entries == -1) {
@@ -866,6 +870,7 @@ TString* TStatToolkit::FitPlaneConstrain(TTree *tree, const char* drawCommand, c
    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] ; 
+   for (Int_t i=0; i<dim+1; i++) values[i]=NULL; 
    //
    entries = tree->Draw(ferr.Data(), cutStr.Data(), "goff",  stop-start, start);
    if (entries == -1) {
@@ -970,6 +975,7 @@ TString* TStatToolkit::FitPlaneFixed(TTree *tree, const char* drawCommand, const
    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] ; 
+   for (Int_t i=0; i<dim+1; i++) values[i]=NULL; 
    //
    entries = tree->Draw(ferr.Data(), cutStr.Data(), "goff",  stop-start, start);
    if (entries == -1) {
@@ -1031,7 +1037,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 +1059,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 +1136,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 +1144,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,47 +1169,60 @@ 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){
+TGraph * TStatToolkit::MakeGraphSparse(TTree * tree, const char * expr, const char * cut, Int_t mstyle, Int_t mcolor){
   //
   // Make a sparse draw of the variables
-  //
+  // Writen by Weilin.Yu
   const Int_t entries =  tree->Draw(expr,cut,"goff");
+  if (entries<=0) {
+    TStatToolkit t;
+    t.Error("TStatToolkit::MakeGraphSparse",Form("Empty or Not valid expression (%s) or cut *%s)", expr,cut));
+    return 0;
+  }
   //  TGraph * graph = (TGraph*)gPad->GetPrimitive("Graph"); // 2D
-  TGraph * graph = new TGraph (entries, tree->GetV2(),tree->GetV1());
+  TGraph * graph = 0;
+  if (tree->GetV3()) graph = new TGraphErrors (entries, tree->GetV2(),tree->GetV1(),0,tree->GetV3());
+  graph =  new TGraphErrors (entries, tree->GetV2(),tree->GetV1(),0,0);
+  graph->SetMarkerStyle(mstyle); 
+  graph->SetMarkerColor(mcolor);
   //
-  Int_t *index = new Int_t[entries];
+  Int_t *index = new Int_t[entries*4];
   TMath::Sort(entries,graph->GetX(),index,kFALSE);
   
-  Double_t *VV = new Double_t[entries];
+  Double_t *tempArray = new Double_t[entries];
 
   Double_t count = 0.5;
-  vector<Int_t> vrun;
-  VV[index[0]] = count;
-  vrun.push_back(graph->GetX()[index[0]]);
+  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]])
-      VV[index[i]] = count; 
+      tempArray[index[i]] = count; 
     else if(graph->GetX()[index[i]]!=graph->GetX()[index[i-1]]){
       count++;
-      VV[index[i]] = count;
-      vrun.push_back(graph->GetX()[index[i]]);
+      icount++;
+      tempArray[index[i]] = count;
+      vrun[icount]=graph->GetX()[index[i]];
     }
   }
   
@@ -1208,20 +1232,24 @@ TGraph * TStatToolkit::MakeGraphSparse(TTree * tree, const char * expr, const ch
     newBins[i] = i;
   }
   
-  TGraph *graphNew = new TGraph(entries,VV,graph->GetY());
+  TGraph *graphNew = 0;
+  if (tree->GetV3()) graphNew = new TGraphErrors(entries,tempArray,graph->GetY(),0,tree->GetV3());
+  else
+    graphNew = new TGraphErrors(entries,tempArray,graph->GetY(),0,0);
   graphNew->GetXaxis()->Set(newNbins,newBins);
   
   Char_t xName[50];
-  Double_t bin_unit = graphNew->GetXaxis()->GetNbins()/count;
   for(Int_t i=0;i<count;i++){
-    snprintf(xName,50,"%d",vrun.at(i));
+    snprintf(xName,50,"%d",Int_t(vrun[i]));
     graphNew->GetXaxis()->SetBinLabel(i+1,xName);
   }
   graphNew->GetHistogram()->SetTitle("");
-  
-  delete [] VV;
+  graphNew->SetMarkerStyle(mstyle); 
+  graphNew->SetMarkerColor(mcolor);
+  delete [] tempArray;
   delete [] index;
   delete [] newBins;
+  delete [] vrun;
   return graphNew;
 }