1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
17 ///////////////////////////////////////////////////////////////////////////
20 // Subset of matheamtical functions not included in the TMath
23 ///////////////////////////////////////////////////////////////////////////
25 #include "AliMathBase.h"
26 #include "Riostream.h"
29 #include "TLinearFitter.h"
32 // includes neccessary for test functions
37 #include "TStopwatch.h"
38 #include "TTreeStream.h"
40 ClassImp(AliMathBase) // Class implementation to enable ROOT I/O
42 AliMathBase::AliMathBase() : TObject()
45 // Default constructor
48 ///////////////////////////////////////////////////////////////////////////
49 AliMathBase::~AliMathBase()
57 //_____________________________________________________________________________
58 void AliMathBase::EvaluateUni(Int_t nvectors, Double_t *data, Double_t &mean
59 , Double_t &sigma, Int_t hh)
62 // Robust estimator in 1D case MI version - (faster than ROOT version)
64 // For the univariate case
65 // estimates of location and scatter are returned in mean and sigma parameters
66 // the algorithm works on the same principle as in multivariate case -
67 // it finds a subset of size hh with smallest sigma, and then returns mean and
68 // sigma of this subset
73 Double_t faclts[]={2.6477,2.5092,2.3826,2.2662,2.1587,2.0589,1.9660,1.879,1.7973,1.7203,1.6473};
74 Int_t *index=new Int_t[nvectors];
75 TMath::Sort(nvectors, data, index, kFALSE);
77 Int_t nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11);
78 Double_t factor = faclts[TMath::Max(0,nquant-1)];
83 Double_t bestmean = 0;
84 Double_t bestsigma = data[index[nvectors-1]]-data[index[0]]; // maximal possible sigma
85 for (Int_t i=0; i<hh; i++){
86 sumx += data[index[i]];
87 sumx2 += data[index[i]]*data[index[i]];
90 Double_t norm = 1./Double_t(hh);
91 Double_t norm2 = 1./Double_t(hh-1);
92 for (Int_t i=hh; i<nvectors; i++){
93 Double_t cmean = sumx*norm;
94 Double_t csigma = (sumx2 - hh*cmean*cmean)*norm2;
95 if (csigma<bestsigma){
101 sumx += data[index[i]]-data[index[i-hh]];
102 sumx2 += data[index[i]]*data[index[i]]-data[index[i-hh]]*data[index[i-hh]];
105 Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma));
114 void AliMathBase::EvaluateUniExternal(Int_t nvectors, Double_t *data, Double_t &mean, Double_t &sigma, Int_t hh, Float_t externalfactor)
116 // Modified version of ROOT robust EvaluateUni
117 // robust estimator in 1D case MI version
118 // added external factor to include precision of external measurement
123 Double_t faclts[]={2.6477,2.5092,2.3826,2.2662,2.1587,2.0589,1.9660,1.879,1.7973,1.7203,1.6473};
124 Int_t *index=new Int_t[nvectors];
125 TMath::Sort(nvectors, data, index, kFALSE);
127 Int_t nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11);
128 Double_t factor = faclts[0];
130 // fix proper normalization - Anja
131 factor = faclts[nquant-1];
138 Int_t bestindex = -1;
139 Double_t bestmean = 0;
140 Double_t bestsigma = -1;
141 for (Int_t i=0; i<hh; i++){
142 sumx += data[index[i]];
143 sumx2 += data[index[i]]*data[index[i]];
146 Double_t kfactor = 2.*externalfactor - externalfactor*externalfactor;
147 Double_t norm = 1./Double_t(hh);
148 for (Int_t i=hh; i<nvectors; i++){
149 Double_t cmean = sumx*norm;
150 Double_t csigma = (sumx2*norm - cmean*cmean*kfactor);
151 if (csigma<bestsigma || bestsigma<0){
158 sumx += data[index[i]]-data[index[i-hh]];
159 sumx2 += data[index[i]]*data[index[i]]-data[index[i-hh]]*data[index[i-hh]];
162 Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma));
169 //_____________________________________________________________________________
170 Int_t AliMathBase::Freq(Int_t n, const Int_t *inlist
171 , Int_t *outlist, Bool_t down)
174 // Sort eleements according occurancy
175 // The size of output array has is 2*n
178 Int_t * sindexS = new Int_t[n]; // temp array for sorting
179 Int_t * sindexF = new Int_t[2*n];
180 for (Int_t i=0;i<n;i++) sindexF[i]=0;
182 TMath::Sort(n,inlist, sindexS, down);
183 Int_t last = inlist[sindexS[0]];
190 for(Int_t i=1;i<n; i++){
191 val = inlist[sindexS[i]];
192 if (last == val) sindexF[countPos]++;
195 sindexF[countPos+n] = val;
200 if (last==val) countPos++;
201 // sort according frequency
202 TMath::Sort(countPos, sindexF, sindexS, kTRUE);
203 for (Int_t i=0;i<countPos;i++){
204 outlist[2*i ] = sindexF[sindexS[i]+n];
205 outlist[2*i+1] = sindexF[sindexS[i]];
214 //___AliMathBase__________________________________________________________________________
215 void AliMathBase::TruncatedMean(TH1F * his, TVectorD *param, Float_t down, Float_t up, Bool_t verbose){
219 Int_t nbins = his->GetNbinsX();
220 Float_t nentries = his->GetEntries();
225 for (Int_t ibin=1;ibin<nbins; ibin++){
226 ncumul+= his->GetBinContent(ibin);
227 Float_t fraction = Float_t(ncumul)/Float_t(nentries);
228 if (fraction>down && fraction<up){
229 sum+=his->GetBinContent(ibin);
230 mean+=his->GetBinCenter(ibin)*his->GetBinContent(ibin);
231 sigma2+=his->GetBinCenter(ibin)*his->GetBinCenter(ibin)*his->GetBinContent(ibin);
235 sigma2= TMath::Sqrt(TMath::Abs(sigma2/sum-mean*mean));
237 (*param)[0] = his->GetMaximum();
239 (*param)[2] = sigma2;
242 if (verbose) printf("Mean\t%f\t Sigma2\t%f\n", mean,sigma2);
245 void AliMathBase::LTM(TH1F * his, TVectorD *param , Float_t fraction, Bool_t verbose){
249 Int_t nbins = his->GetNbinsX();
250 Int_t nentries = (Int_t)his->GetEntries();
251 Double_t *data = new Double_t[nentries];
253 for (Int_t ibin=1;ibin<nbins; ibin++){
254 Float_t entriesI = his->GetBinContent(ibin);
255 Float_t xcenter= his->GetBinCenter(ibin);
256 for (Int_t ic=0; ic<entriesI; ic++){
257 if (npoints<nentries){
258 data[npoints]= xcenter;
263 Double_t mean, sigma;
264 Int_t npoints2=TMath::Min(Int_t(fraction*Float_t(npoints)),npoints-1);
265 npoints2=TMath::Max(Int_t(0.5*Float_t(npoints)),npoints2);
266 AliMathBase::EvaluateUni(npoints, data, mean,sigma,npoints2);
268 if (verbose) printf("Mean\t%f\t Sigma2\t%f\n", mean,sigma);if (param){
269 (*param)[0] = his->GetMaximum();
275 Double_t AliMathBase::FitGaus(TH1F* his, TVectorD *param, TMatrixD *matrix, Float_t xmin, Float_t xmax, Bool_t verbose){
277 // Fit histogram with gaussian function
280 // return value- chi2 - if negative ( not enough points)
281 // his - input histogram
282 // param - vector with parameters
283 // xmin, xmax - range to fit - if xmin=xmax=0 - the full histogram range used
285 // 1. Step - make logarithm
286 // 2. Linear fit (parabola) - more robust - always converge
287 // 3. In case of small statistic bins are averaged
289 static TLinearFitter fitter(3,"pol2");
293 if (his->GetMaximum()<4) return -1;
294 if (his->GetEntries()<12) return -1;
295 if (his->GetRMS()<mat.GetTol()) return -1;
296 Float_t maxEstimate = his->GetEntries()*his->GetBinWidth(1)/TMath::Sqrt((TMath::TwoPi()*his->GetRMS()));
297 Int_t dsmooth = TMath::Nint(6./TMath::Sqrt(maxEstimate));
299 if (maxEstimate<1) return -1;
300 Int_t nbins = his->GetNbinsX();
306 xmin = his->GetXaxis()->GetXmin();
307 xmax = his->GetXaxis()->GetXmax();
309 for (Int_t iter=0; iter<2; iter++){
310 fitter.ClearPoints();
312 for (Int_t ibin=1;ibin<nbins+1; ibin++){
314 Float_t entriesI = his->GetBinContent(ibin);
315 for (Int_t delta = -dsmooth; delta<=dsmooth; delta++){
316 if (ibin+delta>1 &&ibin+delta<nbins-1){
317 entriesI += his->GetBinContent(ibin+delta);
322 Double_t xcenter= his->GetBinCenter(ibin);
323 if (xcenter<xmin || xcenter>xmax) continue;
324 Double_t error=1./TMath::Sqrt(countB);
327 if (par[0]+par[1]*xcenter+par[2]*xcenter*xcenter>20) return 0;
328 cont = TMath::Exp(par[0]+par[1]*xcenter+par[2]*xcenter*xcenter);
329 if (cont>1.) error = 1./TMath::Sqrt(cont*Float_t(countB));
331 if (entriesI>1&&cont>1){
332 fitter.AddPoint(&xcenter,TMath::Log(Float_t(entriesI)),error);
338 fitter.GetParameters(par);
346 fitter.GetParameters(par);
347 fitter.GetCovarianceMatrix(mat);
348 if (TMath::Abs(par[1])<mat.GetTol()) return -1;
349 if (TMath::Abs(par[2])<mat.GetTol()) return -1;
350 Double_t chi2 = fitter.GetChisquare()/Float_t(npoints);
351 //fitter.GetParameters();
352 if (!param) param = new TVectorD(3);
353 if (!matrix) matrix = new TMatrixD(3,3);
354 (*param)[1] = par[1]/(-2.*par[2]);
355 (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
356 (*param)[0] = TMath::Exp(par[0]+ par[1]* (*param)[1] + par[2]*(*param)[1]*(*param)[1]);
361 printf("Chi2=%f\n",chi2);
362 TF1 * f1= new TF1("f1","[0]*exp(-(x-[1])^2/(2*[2]*[2]))",his->GetXaxis()->GetXmin(),his->GetXaxis()->GetXmax());
363 f1->SetParameter(0, (*param)[0]);
364 f1->SetParameter(1, (*param)[1]);
365 f1->SetParameter(2, (*param)[2]);
371 Double_t AliMathBase::FitGaus(Float_t *arr, Int_t nBins, Float_t xMin, Float_t xMax, TVectorD *param, TMatrixD *matrix, Bool_t verbose){
373 // Fit histogram with gaussian function
376 // nbins: size of the array and number of histogram bins
377 // xMin, xMax: histogram range
378 // param: paramters of the fit (0-Constant, 1-Mean, 2-Sigma)
379 // matrix: covariance matrix -- not implemented yet, pass dummy matrix!!!
382 // >0: the chi2 returned by TLinearFitter
383 // -3: only three points have been used for the calculation - no fitter was used
384 // -2: only two points have been used for the calculation - center of gravity was uesed for calculation
385 // -1: only one point has been used for the calculation - center of gravity was uesed for calculation
386 // -4: invalid result!!
389 // 1. Step - make logarithm
390 // 2. Linear fit (parabola) - more robust - always converge
392 static TLinearFitter fitter(3,"pol2");
393 static TMatrixD mat(3,3);
394 static Double_t kTol = mat.GetTol();
395 fitter.StoreData(kFALSE);
396 fitter.ClearPoints();
401 Float_t rms = TMath::RMS(nBins,arr);
402 Float_t max = TMath::MaxElement(nBins,arr);
403 Float_t binWidth = (xMax-xMin)/(Float_t)nBins;
412 for (Int_t i=0; i<nBins; i++){
414 if (arr[i]>0) nfilled++;
417 if (max<4) return -4;
418 if (entries<12) return -4;
419 if (rms<kTol) return -4;
425 for (Int_t ibin=0;ibin<nBins; ibin++){
426 Float_t entriesI = arr[ibin];
428 Double_t xcenter = xMin+(ibin+0.5)*binWidth;
430 Float_t error = 1./TMath::Sqrt(entriesI);
431 Float_t val = TMath::Log(Float_t(entriesI));
432 fitter.AddPoint(&xcenter,val,error);
435 A(npoints,1)=xcenter;
436 A(npoints,2)=xcenter*xcenter;
438 meanCOG+=xcenter*entriesI;
439 rms2COG +=xcenter*entriesI*xcenter;
450 //analytic calculation of the parameters for three points
459 // use fitter for more than three points
461 fitter.GetParameters(par);
462 fitter.GetCovarianceMatrix(mat);
463 chi2 = fitter.GetChisquare()/Float_t(npoints);
465 if (TMath::Abs(par[1])<kTol) return -4;
466 if (TMath::Abs(par[2])<kTol) return -4;
468 if (!param) param = new TVectorD(3);
469 if (!matrix) matrix = new TMatrixD(3,3); // !!!!might be a memory leek. use dummy matrix pointer to call this function!
471 (*param)[1] = par[1]/(-2.*par[2]);
472 (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
473 Double_t lnparam0 = par[0]+ par[1]* (*param)[1] + par[2]*(*param)[1]*(*param)[1];
474 if ( lnparam0>307 ) return -4;
475 (*param)[0] = TMath::Exp(lnparam0);
480 printf("Chi2=%f\n",chi2);
481 TF1 * f1= new TF1("f1","[0]*exp(-(x-[1])^2/(2*[2]*[2]))",xMin,xMax);
482 f1->SetParameter(0, (*param)[0]);
483 f1->SetParameter(1, (*param)[1]);
484 f1->SetParameter(2, (*param)[2]);
491 //use center of gravity for 2 points
495 (*param)[1] = meanCOG;
496 (*param)[2] = TMath::Sqrt(TMath::Abs(meanCOG*meanCOG-rms2COG));
502 (*param)[1] = meanCOG;
503 (*param)[2] = binWidth/TMath::Sqrt(12);
511 Float_t AliMathBase::GetCOG(Short_t *arr, Int_t nBins, Float_t xMin, Float_t xMax, Float_t *rms, Float_t *sum)
514 // calculate center of gravity rms and sum for array 'arr' with nBins an a x range xMin to xMax
515 // return COG; in case of failure return xMin
522 Float_t binWidth = (xMax-xMin)/(Float_t)nBins;
524 for (Int_t ibin=0; ibin<nBins; ibin++){
525 Float_t entriesI = (Float_t)arr[ibin];
526 Double_t xcenter = xMin+(ibin+0.5)*binWidth;
528 meanCOG += xcenter*entriesI;
529 rms2COG += xcenter*entriesI*xcenter;
534 if ( sumCOG == 0 ) return xMin;
539 (*rms) = TMath::Sqrt(TMath::Abs(meanCOG*meanCOG-rms2COG));
540 if ( npoints == 1 ) (*rms) = binWidth/TMath::Sqrt(12);
551 ///////////////////////////////////////////////////////////////
552 ////////////// TEST functions /////////////////////////
553 ///////////////////////////////////////////////////////////////
559 void AliMathBase::TestGausFit(Int_t nhistos){
561 // Test performance of the parabolic - gaussian fit - compare it with
563 // nhistos - number of histograms to be used for test
565 TTreeSRedirector *pcstream = new TTreeSRedirector("fitdebug.root");
567 Float_t *xTrue = new Float_t[nhistos];
568 Float_t *sTrue = new Float_t[nhistos];
569 TVectorD **par1 = new TVectorD*[nhistos];
570 TVectorD **par2 = new TVectorD*[nhistos];
574 TH1F **h1f = new TH1F*[nhistos];
575 TF1 *myg = new TF1("myg","gaus");
576 TF1 *fit = new TF1("fit","gaus");
580 for (Int_t i=0;i<nhistos; i++){
581 par1[i] = new TVectorD(3);
582 par2[i] = new TVectorD(3);
583 h1f[i] = new TH1F(Form("h1f%d",i),Form("h1f%d",i),20,-10,10);
584 xTrue[i]= gRandom->Rndm();
586 sTrue[i]= .75+gRandom->Rndm()*.5;
587 myg->SetParameters(1,xTrue[i],sTrue[i]);
588 h1f[i]->FillRandom("myg");
594 for (Int_t i=0; i<nhistos; i++){
595 h1f[i]->Fit(fit,"0q");
596 (*par1[i])(0) = fit->GetParameter(0);
597 (*par1[i])(1) = fit->GetParameter(1);
598 (*par1[i])(2) = fit->GetParameter(2);
601 printf("Gaussian fit\t");
605 //AliMathBase gaus fit
606 for (Int_t i=0; i<nhistos; i++){
607 AliMathBase::FitGaus(h1f[i]->GetArray()+1,h1f[i]->GetNbinsX(),h1f[i]->GetXaxis()->GetXmin(),h1f[i]->GetXaxis()->GetXmax(),par2[i],&dummy);
611 printf("Parabolic fit\t");
614 for (Int_t i=0;i<nhistos; i++){
615 Float_t xt = xTrue[i];
616 Float_t st = sTrue[i];
625 for (Int_t i=0;i<nhistos; i++){