]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - STEER/AliMathBase.cxx
Removing the circular dependency between the ESD and STEER libraries by moving the...
[u/mrichter/AliRoot.git] / STEER / AliMathBase.cxx
index 89d9a87c40a4a93f4eff770a3bd99e9b350c7790..21bee1a46d9548c1e36a771841154a1927ca7501 100644 (file)
 #include "AliMathBase.h"
 #include "Riostream.h"
 #include "TH1F.h"
+#include "TH3.h"
 #include "TF1.h"
 #include "TLinearFitter.h"
 
+#include "AliExternalTrackParam.h"
+
 //
 // includes neccessary for test functions
 //
@@ -377,7 +380,7 @@ Double_t  AliMathBase::FitGaus(Float_t *arr, Int_t nBins, Float_t xMin, Float_t
   //  Prameters:
   //     nbins: size of the array and number of histogram bins
   //     xMin, xMax: histogram range
-  //     param: paramters of the fit (0-Constant, 1-Mean, 2-Sigma)
+  //     param: paramters of the fit (0-Constant, 1-Mean, 2-Sigma, 3-Sum)
   //     matrix: covariance matrix -- not implemented yet, pass dummy matrix!!!
   //
   //  Return values:
@@ -415,20 +418,22 @@ Double_t  AliMathBase::FitGaus(Float_t *arr, Int_t nBins, Float_t xMin, Float_t
       entries+=arr[i];
       if (arr[i]>0) nfilled++;
   }
+  (*param)[0] = 0;
+  (*param)[1] = 0;
+  (*param)[2] = 0;
+  (*param)[3] = 0;
 
   if (max<4) return -4;
   if (entries<12) return -4;
   if (rms<kTol) return -4;
 
-  Int_t npoints=0;
-  //
+  (*param)[3] = entries;
 
-  //
+  Int_t npoints=0;
   for (Int_t ibin=0;ibin<nBins; ibin++){
       Float_t entriesI = arr[ibin];
     if (entriesI>1){
       Double_t xcenter = xMin+(ibin+0.5)*binWidth;
-      
       Float_t error    = 1./TMath::Sqrt(entriesI);
       Float_t val = TMath::Log(Float_t(entriesI));
       fitter.AddPoint(&xcenter,val,error);
@@ -444,7 +449,6 @@ Double_t  AliMathBase::FitGaus(Float_t *arr, Int_t nBins, Float_t xMin, Float_t
       npoints++;
     }
   }
-
   
   Double_t chi2 = 0;
   if (npoints>=3){
@@ -467,7 +471,8 @@ Double_t  AliMathBase::FitGaus(Float_t *arr, Int_t nBins, Float_t xMin, Float_t
       if (TMath::Abs(par[1])<kTol) return -4;
       if (TMath::Abs(par[2])<kTol) return -4;
 
-      if (!param)  param  = new TVectorD(3);
+      if (!param)  param  = new TVectorD(4);
+      if ( param->GetNrows()<4 ) param->ResizeTo(4);
       if (!matrix) matrix = new TMatrixD(3,3);  // !!!!might be a memory leek. use dummy matrix pointer to call this function!
 
       (*param)[1] = par[1]/(-2.*par[2]);
@@ -641,4 +646,129 @@ void AliMathBase::TestGausFit(Int_t nhistos){
 
 
 
+TGraph2D * AliMathBase::MakeStat2D(TH3 * his, Int_t delta0, Int_t delta1, Int_t type){
+  //
+  //
+  //
+  // delta - number of bins to integrate
+  // type - 0 - mean value
+
+  TAxis * xaxis  = his->GetXaxis();
+  TAxis * yaxis  = his->GetYaxis();
+  //  TAxis * zaxis  = his->GetZaxis();
+  Int_t   nbinx  = xaxis->GetNbins();
+  Int_t   nbiny  = yaxis->GetNbins();
+  char name[1000];
+  Int_t icount=0;
+  TGraph2D  *graph = new TGraph2D(nbinx*nbiny);
+  TF1 f1("f1","gaus");
+  for (Int_t ix=0; ix<nbinx;ix++)
+    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);
+      TH1 *projection = his->ProjectionZ(name,ix-delta0,ix+delta0,iy-delta1,iy+delta1);
+      Float_t stat= 0;
+      if (type==0) stat = projection->GetMean();
+      if (type==1) stat = projection->GetRMS();
+      if (type==2 || type==3){
+       TVectorD vec(3);
+       AliMathBase::LTM((TH1F*)projection,&vec,0.7);
+       if (type==2) stat= vec[1];
+       if (type==3) stat= vec[0];      
+      }
+      if (type==4|| type==5){
+       projection->Fit(&f1);
+       if (type==4) stat= f1.GetParameter(1);
+       if (type==5) stat= f1.GetParameter(2);
+      }
+      //printf("%d\t%f\t%f\t%f\n", icount,xcenter, ycenter, stat);
+      graph->SetPoint(icount,xcenter, ycenter, stat);
+      icount++;
+    }
+  return graph;
+}
 
+TGraph * AliMathBase::MakeStat1D(TH3 * his, Int_t delta1, Int_t type){
+  //
+  //
+  //
+  // delta - number of bins to integrate
+  // type - 0 - mean value
+
+  TAxis * xaxis  = his->GetXaxis();
+  TAxis * yaxis  = his->GetYaxis();
+  //  TAxis * zaxis  = his->GetZaxis();
+  Int_t   nbinx  = xaxis->GetNbins();
+  Int_t   nbiny  = yaxis->GetNbins();
+  char name[1000];
+  Int_t icount=0;
+  TGraph  *graph = new TGraph(nbinx);
+  TF1 f1("f1","gaus");
+  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);
+    TH1 *projection = his->ProjectionZ(name,ix-delta1,ix+delta1,0,nbiny);
+    Float_t stat= 0;
+    if (type==0) stat = projection->GetMean();
+    if (type==1) stat = projection->GetRMS();
+    if (type==2 || type==3){
+      TVectorD vec(3);
+       AliMathBase::LTM((TH1F*)projection,&vec,0.7);
+       if (type==2) stat= vec[1];
+       if (type==3) stat= vec[0];      
+    }
+    if (type==4|| type==5){
+      projection->Fit(&f1);
+      if (type==4) stat= f1.GetParameter(1);
+      if (type==5) stat= f1.GetParameter(2);
+    }
+      //printf("%d\t%f\t%f\t%f\n", icount,xcenter, ycenter, stat);
+    graph->SetPoint(icount,xcenter, stat);
+    icount++;
+  }
+  return graph;
+}
+
+Double_t AliMathBase::TruncatedGaus(Double_t mean, Double_t sigma, Double_t cutat)
+{
+  // return number generated according to a gaussian distribution N(mean,sigma) truncated at cutat
+  //
+  Double_t value;
+  do{
+    value=gRandom->Gaus(mean,sigma);
+  }while(TMath::Abs(value-mean)>cutat);
+  return value;
+}
+
+Double_t AliMathBase::TruncatedGaus(Double_t mean, Double_t sigma, Double_t leftCut, Double_t rightCut)
+{
+  // return number generated according to a gaussian distribution N(mean,sigma)
+  // truncated at leftCut and rightCut
+  //
+  Double_t value;
+  do{
+    value=gRandom->Gaus(mean,sigma);
+  }while((value-mean)<-leftCut || (value-mean)>rightCut);
+  return value;
+}
+
+Double_t AliMathBase::BetheBlochAleph(Double_t bg,
+         Double_t kp1,
+         Double_t kp2,
+         Double_t kp3,
+         Double_t kp4,
+         Double_t kp5) {
+  //
+  // This is the empirical ALEPH parameterization of the Bethe-Bloch formula.
+  // It is normalized to 1 at the minimum.
+  //
+  // bg - beta*gamma
+  // 
+  // The default values for the kp* parameters are for ALICE TPC.
+  // The returned value is in MIP units
+  //
+
+  return AliExternalTrackParam::BetheBlochAleph(bg,kp1,kp2,kp3,kp4,kp5);
+}