#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
//
Double_t sumx2 =0;
Int_t bestindex = -1;
Double_t bestmean = 0;
- Double_t bestsigma = data[index[nvectors-1]]-data[index[0]]; // maximal possible sigma
+ Double_t bestsigma = (data[index[nvectors-1]]-data[index[0]]+1.); // maximal possible sigma
+ bestsigma *=bestsigma;
+
for (Int_t i=0; i<hh; i++){
sumx += data[index[i]];
sumx2 += data[index[i]]*data[index[i]];
return chi2;
}
-
-Double_t AliMathBase::FitGaus(Int_t size, Float_t *arr, Float_t firstBinX, Float_t binWidth, TVectorD *param, TMatrixD *matrix, Float_t xmin, Float_t xmax, Bool_t verbose){
+Double_t AliMathBase::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
//
// Prameters:
- // return value- chi2 - if negative ( not enough points)
- // his - input histogram
- // param - vector with parameters
- // xmin, xmax - range to fit - if xmin=xmax=0 - the full histogram range used
+ // 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, 3-Sum)
+ // matrix: covariance matrix -- not implemented yet, pass dummy matrix!!!
+ //
+ // Return values:
+ // >0: the chi2 returned by TLinearFitter
+ // -3: only three points have been used for the calculation - no fitter was used
+ // -2: only two points have been used for the calculation - center of gravity was uesed for calculation
+ // -1: only one point has been used for the calculation - center of gravity was uesed for calculation
+ // -4: invalid result!!
+ //
// Fitting:
// 1. Step - make logarithm
// 2. Linear fit (parabola) - more robust - always converge
- // 3. In case of small statistic bins are averaged
//
static TLinearFitter fitter(3,"pol2");
static TMatrixD mat(3,3);
TVectorD sigma(3);
TMatrixD A(3,3);
TMatrixD b(3,1);
- Float_t rms = TMath::RMS(size,arr);
- Float_t max = TMath::MaxElement(size,arr);
+ Float_t rms = TMath::RMS(nBins,arr);
+ Float_t max = TMath::MaxElement(nBins,arr);
+ Float_t binWidth = (xMax-xMin)/(Float_t)nBins;
Float_t meanCOG = 0;
Float_t rms2COG = 0;
Float_t entries = 0;
Int_t nfilled=0;
- for (Int_t i=0; i<size; i++){
+ for (Int_t i=0; i<nBins; i++){
entries+=arr[i];
if (arr[i]>0) nfilled++;
}
+ (*param)[0] = 0;
+ (*param)[1] = 0;
+ (*param)[2] = 0;
+ (*param)[3] = 0;
- if (max<4) return -1;
- if (entries<12) return -1;
- if (rms<kTol) return -1;
-
- Int_t npoints=0;
- //
+ if (max<4) return -4;
+ if (entries<12) return -4;
+ if (rms<kTol) return -4;
+ (*param)[3] = entries;
- if (xmin>=xmax){
- xmin = firstBinX;
- xmax = firstBinX+(size-1)*binWidth;
- }
- //
- for (Int_t ibin=0;ibin<size; ibin++){
- Float_t entriesI = arr[ibin];
+ Int_t npoints=0;
+ for (Int_t ibin=0;ibin<nBins; ibin++){
+ Float_t entriesI = arr[ibin];
if (entriesI>1){
- Double_t xcenter = firstBinX+ibin*binWidth;
- if (xcenter<xmin || xcenter>xmax) continue;
-
+ 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);
+ if (npoints<3){
+ A(npoints,0)=1;
+ A(npoints,1)=xcenter;
+ A(npoints,2)=xcenter*xcenter;
+ b(npoints,0)=val;
+ meanCOG+=xcenter*entriesI;
+ rms2COG +=xcenter*entriesI*xcenter;
+ sumCOG +=entriesI;
+ }
npoints++;
}
}
- if (npoints<=3){
- for (Int_t ibin=0;ibin<size; ibin++){
- Float_t entriesI = arr[ibin];
- if (entriesI>1){
- Double_t xcenter = firstBinX+ibin*binWidth;
- if (xcenter<xmin || xcenter>xmax) continue;
- Float_t val = TMath::Log(Float_t(entriesI));
- // if less than 3 point the fitter will crash!
- // for three points calculate the parameters analytically
- // for one and two points use center of gravity
- A(npoints,0)=1;
- A(npoints,1)=xcenter;
- A(npoints,2)=xcenter*xcenter;
- b(npoints,0)=val;
- meanCOG+=xcenter*val;
- rms2COG +=xcenter*val*xcenter*val;
- sumCOG +=val;
- npoints++;
- }
- }
- }
-
-
-
Double_t chi2 = 0;
if (npoints>=3){
if ( npoints == 3 ){
fitter.GetCovarianceMatrix(mat);
chi2 = fitter.GetChisquare()/Float_t(npoints);
}
- if (TMath::Abs(par[1])<kTol) return -1;
- if (TMath::Abs(par[2])<kTol) return -1;
+ 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]);
(*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]);
+ Double_t lnparam0 = par[0]+ par[1]* (*param)[1] + par[2]*(*param)[1]*(*param)[1];
+ if ( lnparam0>307 ) return -4;
+ (*param)[0] = TMath::Exp(lnparam0);
if (verbose){
par.Print();
mat.Print();
param->Print();
printf("Chi2=%f\n",chi2);
- TF1 * f1= new TF1("f1","[0]*exp(-(x-[1])^2/(2*[2]*[2]))",xmin,xmax);
+ TF1 * f1= new TF1("f1","[0]*exp(-(x-[1])^2/(2*[2]*[2]))",xMin,xMax);
f1->SetParameter(0, (*param)[0]);
f1->SetParameter(1, (*param)[1]);
f1->SetParameter(2, (*param)[2]);
chi2=-2.;
}
if ( npoints == 1 ){
+ meanCOG/=sumCOG;
(*param)[0] = max;
(*param)[1] = meanCOG;
(*param)[2] = binWidth/TMath::Sqrt(12);
}
+Float_t AliMathBase::GetCOG(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
+ // return COG; in case of failure return xMin
+ //
+ Float_t meanCOG = 0;
+ Float_t rms2COG = 0;
+ Float_t sumCOG = 0;
+ Int_t npoints = 0;
+
+ Float_t binWidth = (xMax-xMin)/(Float_t)nBins;
+
+ for (Int_t ibin=0; ibin<nBins; ibin++){
+ Float_t entriesI = (Float_t)arr[ibin];
+ Double_t xcenter = xMin+(ibin+0.5)*binWidth;
+ if ( entriesI>0 ){
+ meanCOG += xcenter*entriesI;
+ rms2COG += xcenter*entriesI*xcenter;
+ sumCOG += entriesI;
+ npoints++;
+ }
+ }
+ if ( sumCOG == 0 ) return xMin;
+ meanCOG/=sumCOG;
+
+ if ( rms ){
+ rms2COG /=sumCOG;
+ (*rms) = TMath::Sqrt(TMath::Abs(meanCOG*meanCOG-rms2COG));
+ if ( npoints == 1 ) (*rms) = binWidth/TMath::Sqrt(12);
+ }
+
+ if ( sum )
+ (*sum) = sumCOG;
+
+ return meanCOG;
+}
+
+
+Double_t AliMathBase::ErfcFast(Double_t x){
+ // Fast implementation of the complementary error function
+ // The error of the approximation is |eps(x)| < 5E-4
+ // See Abramowitz and Stegun, p.299, 7.1.27
+
+ Double_t z = TMath::Abs(x);
+ Double_t ans = 1+z*(0.278393+z*(0.230389+z*(0.000972+z*0.078108)));
+ ans = 1.0/ans;
+ ans *= ans;
+ ans *= ans;
+
+ return (x>=0.0 ? ans : 2.0 - ans);
+}
///////////////////////////////////////////////////////////////
////////////// TEST functions /////////////////////////
s.Start();
//AliMathBase gaus fit
for (Int_t i=0; i<nhistos; i++){
- AliMathBase::FitGaus(20,h1f[i]->GetArray()+1,-9.5,1,par2[i],&dummy);
+ AliMathBase::FitGaus(h1f[i]->GetArray()+1,h1f[i]->GetNbinsX(),h1f[i]->GetXaxis()->GetXmin(),h1f[i]->GetXaxis()->GetXmax(),par2[i],&dummy);
}
s.Stop();
+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);
+}