]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/AliMathBase.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Add fast gaussian fitting routine (Marian)
[u/mrichter/AliRoot.git] / STEER / AliMathBase.cxx
CommitLineData
284050f7 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
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 **************************************************************************/
15
16
17///////////////////////////////////////////////////////////////////////////
18// Class AliMathBase
19//
20// Subset of matheamtical functions not included in the TMath
21//
22
23///////////////////////////////////////////////////////////////////////////
24#include "TMath.h"
25#include "AliMathBase.h"
26#include "Riostream.h"
f6659a9d 27#include "TH1F.h"
28#include "TF1.h"
29#include "TLinearFitter.h"
5608e15a 30
31//
32// includes neccessary for test functions
33//
34
35#include "TSystem.h"
36#include "TRandom.h"
37#include "TStopwatch.h"
38#include "TTreeStream.h"
284050f7 39
40ClassImp(AliMathBase) // Class implementation to enable ROOT I/O
41
42AliMathBase::AliMathBase() : TObject()
43{
5608e15a 44 //
45 // Default constructor
46 //
284050f7 47}
48///////////////////////////////////////////////////////////////////////////
49AliMathBase::~AliMathBase()
50{
5608e15a 51 //
52 // Destructor
53 //
284050f7 54}
55
56
57//_____________________________________________________________________________
58void AliMathBase::EvaluateUni(Int_t nvectors, Double_t *data, Double_t &mean
59 , Double_t &sigma, Int_t hh)
60{
61 //
62 // Robust estimator in 1D case MI version - (faster than ROOT version)
63 //
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
69 //
70
71 if (hh==0)
72 hh=(nvectors+2)/2;
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);
76
77 Int_t nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11);
d9e9045c 78 Double_t factor = faclts[TMath::Max(0,nquant-1)];
284050f7 79
80 Double_t sumx =0;
81 Double_t sumx2 =0;
82 Int_t bestindex = -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]];
88 }
89
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){
96 bestmean = cmean;
97 bestsigma = csigma;
98 bestindex = i-hh;
99 }
100
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]];
103 }
104
105 Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma));
106 mean = bestmean;
107 sigma = bstd;
108 delete [] index;
109
110}
111
112
113
114void AliMathBase::EvaluateUniExternal(Int_t nvectors, Double_t *data, Double_t &mean, Double_t &sigma, Int_t hh, Float_t externalfactor)
115{
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
119 //
120
121 if (hh==0)
122 hh=(nvectors+2)/2;
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);
126 //
127 Int_t nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11);
128 Double_t factor = faclts[0];
129 if (nquant>0){
130 // fix proper normalization - Anja
131 factor = faclts[nquant-1];
132 }
133
134 //
135 //
136 Double_t sumx =0;
137 Double_t sumx2 =0;
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]];
144 }
145 //
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){
152 bestmean = cmean;
153 bestsigma = csigma;
154 bestindex = i-hh;
155 }
156 //
157 //
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]];
160 }
161
162 Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma));
163 mean = bestmean;
164 sigma = bstd;
165 delete [] index;
166}
167
168
169//_____________________________________________________________________________
170Int_t AliMathBase::Freq(Int_t n, const Int_t *inlist
171 , Int_t *outlist, Bool_t down)
172{
173 //
174 // Sort eleements according occurancy
175 // The size of output array has is 2*n
176 //
177
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;
181 //
182 TMath::Sort(n,inlist, sindexS, down);
183 Int_t last = inlist[sindexS[0]];
184 Int_t val = last;
185 sindexF[0] = 1;
186 sindexF[0+n] = last;
187 Int_t countPos = 0;
188 //
189 // find frequency
190 for(Int_t i=1;i<n; i++){
191 val = inlist[sindexS[i]];
192 if (last == val) sindexF[countPos]++;
193 else{
194 countPos++;
195 sindexF[countPos+n] = val;
196 sindexF[countPos]++;
197 last =val;
198 }
199 }
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]];
206 }
207 delete [] sindexS;
208 delete [] sindexF;
209
210 return countPos;
211
212}
f6659a9d 213
214//___AliMathBase__________________________________________________________________________
215void AliMathBase::TruncatedMean(TH1F * his, TVectorD *param, Float_t down, Float_t up, Bool_t verbose){
216 //
217 //
218 //
219 Int_t nbins = his->GetNbinsX();
220 Float_t nentries = his->GetEntries();
221 Float_t sum =0;
222 Float_t mean = 0;
223 Float_t sigma2 = 0;
224 Float_t ncumul=0;
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);
232 }
233 }
234 mean/=sum;
235 sigma2= TMath::Sqrt(TMath::Abs(sigma2/sum-mean*mean));
236 if (param){
237 (*param)[0] = his->GetMaximum();
238 (*param)[1] = mean;
239 (*param)[2] = sigma2;
240
241 }
242 if (verbose) printf("Mean\t%f\t Sigma2\t%f\n", mean,sigma2);
243}
244
245void AliMathBase::LTM(TH1F * his, TVectorD *param , Float_t fraction, Bool_t verbose){
246 //
247 // LTM
248 //
249 Int_t nbins = his->GetNbinsX();
250 Int_t nentries = (Int_t)his->GetEntries();
251 Double_t *data = new Double_t[nentries];
252 Int_t npoints=0;
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;
259 npoints++;
260 }
261 }
262 }
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);
267 delete [] data;
268 if (verbose) printf("Mean\t%f\t Sigma2\t%f\n", mean,sigma);if (param){
269 (*param)[0] = his->GetMaximum();
270 (*param)[1] = mean;
271 (*param)[2] = sigma;
272 }
273}
274
275Double_t AliMathBase::FitGaus(TH1F* his, TVectorD *param, TMatrixD *matrix, Float_t xmin, Float_t xmax, Bool_t verbose){
276 //
277 // Fit histogram with gaussian function
278 //
279 // Prameters:
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
284 // Fitting:
285 // 1. Step - make logarithm
286 // 2. Linear fit (parabola) - more robust - always converge
287 // 3. In case of small statistic bins are averaged
288 //
289 static TLinearFitter fitter(3,"pol2");
290 TVectorD par(3);
291 TVectorD sigma(3);
292 TMatrixD mat(3,3);
293 if (his->GetMaximum()<4) return -1;
294 if (his->GetEntries()<12) return -1;
295 if (his->GetRMS()<mat.GetTol()) return -1;
5608e15a 296 Float_t maxEstimate = his->GetEntries()*his->GetBinWidth(1)/TMath::Sqrt((TMath::TwoPi()*his->GetRMS()));
f6659a9d 297 Int_t dsmooth = TMath::Nint(6./TMath::Sqrt(maxEstimate));
298
299 if (maxEstimate<1) return -1;
300 Int_t nbins = his->GetNbinsX();
301 Int_t npoints=0;
302 //
303
304
305 if (xmin>=xmax){
306 xmin = his->GetXaxis()->GetXmin();
307 xmax = his->GetXaxis()->GetXmax();
308 }
309 for (Int_t iter=0; iter<2; iter++){
310 fitter.ClearPoints();
311 npoints=0;
5608e15a 312 for (Int_t ibin=1;ibin<nbins+1; ibin++){
f6659a9d 313 Int_t countB=1;
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);
318 countB++;
319 }
320 }
321 entriesI/=countB;
322 Double_t xcenter= his->GetBinCenter(ibin);
323 if (xcenter<xmin || xcenter>xmax) continue;
324 Double_t error=1./TMath::Sqrt(countB);
325 Float_t cont=2;
326 if (iter>0){
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));
330 }
331 if (entriesI>1&&cont>1){
332 fitter.AddPoint(&xcenter,TMath::Log(Float_t(entriesI)),error);
333 npoints++;
334 }
335 }
336 if (npoints>3){
337 fitter.Eval();
338 fitter.GetParameters(par);
339 }else{
340 break;
341 }
342 }
343 if (npoints<=3){
344 return -1;
345 }
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]);
357 if (verbose){
358 par.Print();
359 mat.Print();
360 param->Print();
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]);
366 f1->Draw("same");
367 }
368 return chi2;
369}
370
371
5608e15a 372Double_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){
373 //
374 // Fit histogram with gaussian function
375 //
376 // Prameters:
377 // return value- chi2 - if negative ( not enough points)
378 // his - input histogram
379 // param - vector with parameters
380 // xmin, xmax - range to fit - if xmin=xmax=0 - the full histogram range used
381 // Fitting:
382 // 1. Step - make logarithm
383 // 2. Linear fit (parabola) - more robust - always converge
384 // 3. In case of small statistic bins are averaged
385 //
386 static TLinearFitter fitter(3,"pol2");
387 static TMatrixD mat(3,3);
388 static Double_t kTol = mat.GetTol();
389 fitter.StoreData(kFALSE);
390 fitter.ClearPoints();
391 TVectorD par(3);
392 TVectorD sigma(3);
393 TMatrixD A(3,3);
394 TMatrixD b(3,1);
395 Float_t rms = TMath::RMS(size,arr);
396 Float_t max = TMath::MaxElement(size,arr);
397
398 Float_t meanCOG = 0;
399 Float_t rms2COG = 0;
400 Float_t sumCOG = 0;
401
402 Float_t entries = 0;
403 Int_t nfilled=0;
404
405 for (Int_t i=0; i<size; i++){
406 entries+=arr[i];
407 if (arr[i]>0) nfilled++;
408 }
409
410 if (max<4) return -1;
411 if (entries<12) return -1;
412 if (rms<kTol) return -1;
413
414 Int_t npoints=0;
415 //
416
417
418 if (xmin>=xmax){
419 xmin = firstBinX;
420 xmax = firstBinX+(size-1)*binWidth;
421 }
422 //
423 for (Int_t ibin=0;ibin<size; ibin++){
424 Float_t entriesI = arr[ibin];
425 if (entriesI>1){
426 Double_t xcenter = firstBinX+ibin*binWidth;
427 if (xcenter<xmin || xcenter>xmax) continue;
428
429 Float_t error = 1./TMath::Sqrt(entriesI);
430 Float_t val = TMath::Log(Float_t(entriesI));
431 fitter.AddPoint(&xcenter,val,error);
432 npoints++;
433 }
434 }
435
436 if (npoints<=3){
437 for (Int_t ibin=0;ibin<size; ibin++){
438 Float_t entriesI = arr[ibin];
439 if (entriesI>1){
440 Double_t xcenter = firstBinX+ibin*binWidth;
441 if (xcenter<xmin || xcenter>xmax) continue;
442 Float_t val = TMath::Log(Float_t(entriesI));
443 // if less than 3 point the fitter will crash!
444 // for three points calculate the parameters analytically
445 // for one and two points use center of gravity
446 A(npoints,0)=1;
447 A(npoints,1)=xcenter;
448 A(npoints,2)=xcenter*xcenter;
449 b(npoints,0)=val;
450 meanCOG+=xcenter*val;
451 rms2COG +=xcenter*val*xcenter*val;
452 sumCOG +=val;
453 npoints++;
454 }
455 }
456 }
457
458
459
460 Double_t chi2 = 0;
461 if (npoints>=3){
462 if ( npoints == 3 ){
463 //analytic calculation of the parameters for three points
464 A.Invert();
465 TMatrixD res(1,3);
466 res.Mult(A,b);
467 par[0]=res(0,0);
468 par[1]=res(0,1);
469 par[2]=res(0,2);
470 chi2 = -3.;
471 } else {
472 // use fitter for more than three points
473 fitter.Eval();
474 fitter.GetParameters(par);
475 fitter.GetCovarianceMatrix(mat);
476 chi2 = fitter.GetChisquare()/Float_t(npoints);
477 }
478 if (TMath::Abs(par[1])<kTol) return -1;
479 if (TMath::Abs(par[2])<kTol) return -1;
480
481 if (!param) param = new TVectorD(3);
482 if (!matrix) matrix = new TMatrixD(3,3); // !!!!might be a memory leek. use dummy matrix pointer to call this function!
483
484 (*param)[1] = par[1]/(-2.*par[2]);
485 (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
486 (*param)[0] = TMath::Exp(par[0]+ par[1]* (*param)[1] + par[2]*(*param)[1]*(*param)[1]);
487 if (verbose){
488 par.Print();
489 mat.Print();
490 param->Print();
491 printf("Chi2=%f\n",chi2);
492 TF1 * f1= new TF1("f1","[0]*exp(-(x-[1])^2/(2*[2]*[2]))",xmin,xmax);
493 f1->SetParameter(0, (*param)[0]);
494 f1->SetParameter(1, (*param)[1]);
495 f1->SetParameter(2, (*param)[2]);
496 f1->Draw("same");
497 }
498 return chi2;
499 }
500
501 if (npoints == 2){
502 //use center of gravity for 2 points
503 meanCOG/=sumCOG;
504 rms2COG /=sumCOG;
505 (*param)[0] = max;
506 (*param)[1] = meanCOG;
507 (*param)[2] = TMath::Sqrt(TMath::Abs(meanCOG*meanCOG-rms2COG));
508 chi2=-2.;
509 }
510 if ( npoints == 1 ){
511 (*param)[0] = max;
512 (*param)[1] = meanCOG;
513 (*param)[2] = binWidth/TMath::Sqrt(12);
514 chi2=-1.;
515 }
516 return chi2;
517
518}
519
520
521
522///////////////////////////////////////////////////////////////
523////////////// TEST functions /////////////////////////
524///////////////////////////////////////////////////////////////
525
526
527
528
529
530void AliMathBase::TestGausFit(Int_t nhistos){
531 //
532 // Test performance of the parabolic - gaussian fit - compare it with
533 // ROOT gauss fit
534 // nhistos - number of histograms to be used for test
535 //
536 TTreeSRedirector *pcstream = new TTreeSRedirector("fitdebug.root");
537
538 Float_t *xTrue = new Float_t[nhistos];
539 Float_t *sTrue = new Float_t[nhistos];
540 TVectorD **par1 = new TVectorD*[nhistos];
541 TVectorD **par2 = new TVectorD*[nhistos];
542 TMatrixD dummy(3,3);
543
544
545 TH1F **h1f = new TH1F*[nhistos];
546 TF1 *myg = new TF1("myg","gaus");
547 TF1 *fit = new TF1("fit","gaus");
548 gRandom->SetSeed(0);
549
550 //init
551 for (Int_t i=0;i<nhistos; i++){
552 par1[i] = new TVectorD(3);
553 par2[i] = new TVectorD(3);
554 h1f[i] = new TH1F(Form("h1f%d",i),Form("h1f%d",i),20,-10,10);
555 xTrue[i]= gRandom->Rndm();
556 gSystem->Sleep(2);
557 sTrue[i]= .75+gRandom->Rndm()*.5;
558 myg->SetParameters(1,xTrue[i],sTrue[i]);
559 h1f[i]->FillRandom("myg");
560 }
561
562 TStopwatch s;
563 s.Start();
564 //standard gaus fit
565 for (Int_t i=0; i<nhistos; i++){
566 h1f[i]->Fit(fit,"0q");
567 (*par1[i])(0) = fit->GetParameter(0);
568 (*par1[i])(1) = fit->GetParameter(1);
569 (*par1[i])(2) = fit->GetParameter(2);
570 }
571 s.Stop();
572 printf("Gaussian fit\t");
573 s.Print();
574
575 s.Start();
576 //AliMathBase gaus fit
577 for (Int_t i=0; i<nhistos; i++){
578 AliMathBase::FitGaus(20,h1f[i]->GetArray()+1,-9.5,1,par2[i],&dummy);
579 }
580
581 s.Stop();
582 printf("Parabolic fit\t");
583 s.Print();
584 //write stream
585 for (Int_t i=0;i<nhistos; i++){
586 Float_t xt = xTrue[i];
587 Float_t st = sTrue[i];
588 (*pcstream)<<"data"
589 <<"xTrue="<<xt
590 <<"sTrue="<<st
591 <<"pg.="<<(par1[i])
592 <<"pa.="<<(par2[i])
593 <<"\n";
594 }
595 //delete pointers
596 for (Int_t i=0;i<nhistos; i++){
597 delete par1[i];
598 delete par2[i];
599 delete h1f[i];
600 }
601 delete pcstream;
602 delete []h1f;
603 delete []xTrue;
604 delete []sTrue;
605 //
606 delete []par1;
607 delete []par2;
608
609}
610
611
612
613