#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]];
// 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:
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);
npoints++;
}
}
-
Double_t chi2 = 0;
if (npoints>=3){
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]);
}
+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 /////////////////////////
+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);
+}