#include "TMath.h"
#include "AliMathBase.h"
#include "Riostream.h"
+#include "TH1F.h"
+#include "TF1.h"
+#include "TLinearFitter.h"
ClassImp(AliMathBase) // Class implementation to enable ROOT I/O
return countPos;
}
+
+//___AliMathBase__________________________________________________________________________
+void AliMathBase::TruncatedMean(TH1F * his, TVectorD *param, Float_t down, Float_t up, Bool_t verbose){
+ //
+ //
+ //
+ Int_t nbins = his->GetNbinsX();
+ Float_t nentries = his->GetEntries();
+ Float_t sum =0;
+ Float_t mean = 0;
+ Float_t sigma2 = 0;
+ Float_t ncumul=0;
+ for (Int_t ibin=1;ibin<nbins; ibin++){
+ ncumul+= his->GetBinContent(ibin);
+ Float_t fraction = Float_t(ncumul)/Float_t(nentries);
+ if (fraction>down && fraction<up){
+ sum+=his->GetBinContent(ibin);
+ mean+=his->GetBinCenter(ibin)*his->GetBinContent(ibin);
+ sigma2+=his->GetBinCenter(ibin)*his->GetBinCenter(ibin)*his->GetBinContent(ibin);
+ }
+ }
+ mean/=sum;
+ sigma2= TMath::Sqrt(TMath::Abs(sigma2/sum-mean*mean));
+ if (param){
+ (*param)[0] = his->GetMaximum();
+ (*param)[1] = mean;
+ (*param)[2] = sigma2;
+
+ }
+ if (verbose) printf("Mean\t%f\t Sigma2\t%f\n", mean,sigma2);
+}
+
+void AliMathBase::LTM(TH1F * his, TVectorD *param , Float_t fraction, Bool_t verbose){
+ //
+ // LTM
+ //
+ Int_t nbins = his->GetNbinsX();
+ Int_t nentries = (Int_t)his->GetEntries();
+ Double_t *data = new Double_t[nentries];
+ Int_t npoints=0;
+ for (Int_t ibin=1;ibin<nbins; ibin++){
+ Float_t entriesI = his->GetBinContent(ibin);
+ Float_t xcenter= his->GetBinCenter(ibin);
+ for (Int_t ic=0; ic<entriesI; ic++){
+ if (npoints<nentries){
+ data[npoints]= xcenter;
+ npoints++;
+ }
+ }
+ }
+ Double_t mean, sigma;
+ Int_t npoints2=TMath::Min(Int_t(fraction*Float_t(npoints)),npoints-1);
+ npoints2=TMath::Max(Int_t(0.5*Float_t(npoints)),npoints2);
+ AliMathBase::EvaluateUni(npoints, data, mean,sigma,npoints2);
+ delete [] data;
+ if (verbose) printf("Mean\t%f\t Sigma2\t%f\n", mean,sigma);if (param){
+ (*param)[0] = his->GetMaximum();
+ (*param)[1] = mean;
+ (*param)[2] = sigma;
+ }
+}
+
+Double_t AliMathBase::FitGaus(TH1F* his, TVectorD *param, TMatrixD *matrix, Float_t xmin, Float_t xmax, 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
+ // 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");
+ TVectorD par(3);
+ TVectorD sigma(3);
+ TMatrixD mat(3,3);
+ if (his->GetMaximum()<4) return -1;
+ if (his->GetEntries()<12) return -1;
+ if (his->GetRMS()<mat.GetTol()) return -1;
+ Float_t maxEstimate = his->GetEntries()*his->GetBinWidth(1)/TMath::Sqrt((2.*TMath::Pi()*his->GetRMS()));
+ Int_t dsmooth = TMath::Nint(6./TMath::Sqrt(maxEstimate));
+
+ if (maxEstimate<1) return -1;
+ Int_t nbins = his->GetNbinsX();
+ Int_t npoints=0;
+ //
+
+
+ if (xmin>=xmax){
+ xmin = his->GetXaxis()->GetXmin();
+ xmax = his->GetXaxis()->GetXmax();
+ }
+ for (Int_t iter=0; iter<2; iter++){
+ fitter.ClearPoints();
+ npoints=0;
+ for (Int_t ibin=1;ibin<nbins-1; ibin++){
+ Int_t countB=1;
+ Float_t entriesI = his->GetBinContent(ibin);
+ for (Int_t delta = -dsmooth; delta<=dsmooth; delta++){
+ if (ibin+delta>1 &&ibin+delta<nbins-1){
+ entriesI += his->GetBinContent(ibin+delta);
+ countB++;
+ }
+ }
+ entriesI/=countB;
+ Double_t xcenter= his->GetBinCenter(ibin);
+ if (xcenter<xmin || xcenter>xmax) continue;
+ Double_t error=1./TMath::Sqrt(countB);
+ Float_t cont=2;
+ if (iter>0){
+ if (par[0]+par[1]*xcenter+par[2]*xcenter*xcenter>20) return 0;
+ cont = TMath::Exp(par[0]+par[1]*xcenter+par[2]*xcenter*xcenter);
+ if (cont>1.) error = 1./TMath::Sqrt(cont*Float_t(countB));
+ }
+ if (entriesI>1&&cont>1){
+ fitter.AddPoint(&xcenter,TMath::Log(Float_t(entriesI)),error);
+ npoints++;
+ }
+ }
+ if (npoints>3){
+ fitter.Eval();
+ fitter.GetParameters(par);
+ }else{
+ break;
+ }
+ }
+ if (npoints<=3){
+ return -1;
+ }
+ fitter.GetParameters(par);
+ fitter.GetCovarianceMatrix(mat);
+ if (TMath::Abs(par[1])<mat.GetTol()) return -1;
+ if (TMath::Abs(par[2])<mat.GetTol()) return -1;
+ Double_t chi2 = fitter.GetChisquare()/Float_t(npoints);
+ //fitter.GetParameters();
+ if (!param) param = new TVectorD(3);
+ if (!matrix) matrix = new TMatrixD(3,3);
+ (*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]);
+ 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]))",his->GetXaxis()->GetXmin(),his->GetXaxis()->GetXmax());
+ f1->SetParameter(0, (*param)[0]);
+ f1->SetParameter(1, (*param)[1]);
+ f1->SetParameter(2, (*param)[2]);
+ f1->Draw("same");
+ }
+ return chi2;
+}
+
+
#include "TObject.h"
+#include "TVectorD.h"
+#include "TMatrixD.h"
+class TH1F;
class AliMathBase : public TObject
{
static void EvaluateUni(Int_t nvectors, Double_t *data, Double_t &mean, Double_t &sigma, Int_t hh);
static void EvaluateUniExternal(Int_t nvectors, Double_t *data, Double_t &mean, Double_t &sigma, Int_t hh, Float_t externalfactor=1);
static Int_t Freq(Int_t n, const Int_t *inlist, Int_t *outlist, Bool_t down);
-
+ static void TruncatedMean(TH1F * his, TVectorD *param, Float_t down=0, Float_t up=1.0, Bool_t verbose=kFALSE);
+ static void LTM(TH1F * his, TVectorD *param=0 , Float_t fraction=1, Bool_t verbose=kFALSE);
+ static Double_t FitGaus(TH1F* his, TVectorD *param=0, TMatrixD *matrix=0, Float_t xmin=0, Float_t xmax=0, Bool_t verbose=kFALSE);
+
ClassDef(AliMathBase,0) // Various mathematical tools for physics analysis - which are not included in ROOT TMath
};