]>
Commit | Line | Data |
---|---|---|
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" |
3392b4c9 | 28 | #include "TH3.h" |
f6659a9d | 29 | #include "TF1.h" |
30 | #include "TLinearFitter.h" | |
a8752eef | 31 | #include "AliLog.h" |
5608e15a | 32 | |
9c56b409 | 33 | #include "AliExternalTrackParam.h" |
34 | ||
5608e15a | 35 | // |
36 | // includes neccessary for test functions | |
37 | // | |
38 | ||
39 | #include "TSystem.h" | |
40 | #include "TRandom.h" | |
41 | #include "TStopwatch.h" | |
42 | #include "TTreeStream.h" | |
284050f7 | 43 | |
44 | ClassImp(AliMathBase) // Class implementation to enable ROOT I/O | |
45 | ||
46 | AliMathBase::AliMathBase() : TObject() | |
47 | { | |
5608e15a | 48 | // |
49 | // Default constructor | |
50 | // | |
284050f7 | 51 | } |
52 | /////////////////////////////////////////////////////////////////////////// | |
53 | AliMathBase::~AliMathBase() | |
54 | { | |
5608e15a | 55 | // |
56 | // Destructor | |
57 | // | |
284050f7 | 58 | } |
59 | ||
60 | ||
61 | //_____________________________________________________________________________ | |
62 | void AliMathBase::EvaluateUni(Int_t nvectors, Double_t *data, Double_t &mean | |
63 | , Double_t &sigma, Int_t hh) | |
64 | { | |
65 | // | |
66 | // Robust estimator in 1D case MI version - (faster than ROOT version) | |
67 | // | |
68 | // For the univariate case | |
69 | // estimates of location and scatter are returned in mean and sigma parameters | |
70 | // the algorithm works on the same principle as in multivariate case - | |
71 | // it finds a subset of size hh with smallest sigma, and then returns mean and | |
72 | // sigma of this subset | |
73 | // | |
74 | ||
a8752eef | 75 | if (nvectors<2) { |
76 | AliErrorClass(Form("nvectors = %d, should be > 1",nvectors)); | |
77 | return; | |
78 | } | |
79 | if (hh<2) | |
284050f7 | 80 | hh=(nvectors+2)/2; |
81 | 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}; | |
82 | Int_t *index=new Int_t[nvectors]; | |
83 | TMath::Sort(nvectors, data, index, kFALSE); | |
84 | ||
85 | Int_t nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11); | |
d9e9045c | 86 | Double_t factor = faclts[TMath::Max(0,nquant-1)]; |
284050f7 | 87 | |
88 | Double_t sumx =0; | |
89 | Double_t sumx2 =0; | |
90 | Int_t bestindex = -1; | |
91 | Double_t bestmean = 0; | |
07d955de | 92 | Double_t bestsigma = (data[index[nvectors-1]]-data[index[0]]+1.); // maximal possible sigma |
93 | bestsigma *=bestsigma; | |
94 | ||
284050f7 | 95 | for (Int_t i=0; i<hh; i++){ |
96 | sumx += data[index[i]]; | |
97 | sumx2 += data[index[i]]*data[index[i]]; | |
98 | } | |
99 | ||
100 | Double_t norm = 1./Double_t(hh); | |
101 | Double_t norm2 = 1./Double_t(hh-1); | |
102 | for (Int_t i=hh; i<nvectors; i++){ | |
103 | Double_t cmean = sumx*norm; | |
104 | Double_t csigma = (sumx2 - hh*cmean*cmean)*norm2; | |
105 | if (csigma<bestsigma){ | |
106 | bestmean = cmean; | |
107 | bestsigma = csigma; | |
108 | bestindex = i-hh; | |
109 | } | |
110 | ||
111 | sumx += data[index[i]]-data[index[i-hh]]; | |
112 | sumx2 += data[index[i]]*data[index[i]]-data[index[i-hh]]*data[index[i-hh]]; | |
113 | } | |
114 | ||
115 | Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma)); | |
116 | mean = bestmean; | |
117 | sigma = bstd; | |
118 | delete [] index; | |
119 | ||
120 | } | |
121 | ||
122 | ||
123 | ||
124 | void AliMathBase::EvaluateUniExternal(Int_t nvectors, Double_t *data, Double_t &mean, Double_t &sigma, Int_t hh, Float_t externalfactor) | |
125 | { | |
126 | // Modified version of ROOT robust EvaluateUni | |
127 | // robust estimator in 1D case MI version | |
128 | // added external factor to include precision of external measurement | |
129 | // | |
130 | ||
131 | if (hh==0) | |
132 | hh=(nvectors+2)/2; | |
133 | 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}; | |
134 | Int_t *index=new Int_t[nvectors]; | |
135 | TMath::Sort(nvectors, data, index, kFALSE); | |
136 | // | |
137 | Int_t nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11); | |
138 | Double_t factor = faclts[0]; | |
139 | if (nquant>0){ | |
140 | // fix proper normalization - Anja | |
141 | factor = faclts[nquant-1]; | |
142 | } | |
143 | ||
144 | // | |
145 | // | |
146 | Double_t sumx =0; | |
147 | Double_t sumx2 =0; | |
148 | Int_t bestindex = -1; | |
149 | Double_t bestmean = 0; | |
150 | Double_t bestsigma = -1; | |
151 | for (Int_t i=0; i<hh; i++){ | |
152 | sumx += data[index[i]]; | |
153 | sumx2 += data[index[i]]*data[index[i]]; | |
154 | } | |
155 | // | |
156 | Double_t kfactor = 2.*externalfactor - externalfactor*externalfactor; | |
157 | Double_t norm = 1./Double_t(hh); | |
158 | for (Int_t i=hh; i<nvectors; i++){ | |
159 | Double_t cmean = sumx*norm; | |
160 | Double_t csigma = (sumx2*norm - cmean*cmean*kfactor); | |
161 | if (csigma<bestsigma || bestsigma<0){ | |
162 | bestmean = cmean; | |
163 | bestsigma = csigma; | |
164 | bestindex = i-hh; | |
165 | } | |
166 | // | |
167 | // | |
168 | sumx += data[index[i]]-data[index[i-hh]]; | |
169 | sumx2 += data[index[i]]*data[index[i]]-data[index[i-hh]]*data[index[i-hh]]; | |
170 | } | |
171 | ||
172 | Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma)); | |
173 | mean = bestmean; | |
174 | sigma = bstd; | |
175 | delete [] index; | |
176 | } | |
177 | ||
178 | ||
179 | //_____________________________________________________________________________ | |
180 | Int_t AliMathBase::Freq(Int_t n, const Int_t *inlist | |
181 | , Int_t *outlist, Bool_t down) | |
182 | { | |
183 | // | |
184 | // Sort eleements according occurancy | |
185 | // The size of output array has is 2*n | |
186 | // | |
187 | ||
188 | Int_t * sindexS = new Int_t[n]; // temp array for sorting | |
189 | Int_t * sindexF = new Int_t[2*n]; | |
190 | for (Int_t i=0;i<n;i++) sindexF[i]=0; | |
191 | // | |
192 | TMath::Sort(n,inlist, sindexS, down); | |
193 | Int_t last = inlist[sindexS[0]]; | |
194 | Int_t val = last; | |
195 | sindexF[0] = 1; | |
196 | sindexF[0+n] = last; | |
197 | Int_t countPos = 0; | |
198 | // | |
199 | // find frequency | |
200 | for(Int_t i=1;i<n; i++){ | |
201 | val = inlist[sindexS[i]]; | |
202 | if (last == val) sindexF[countPos]++; | |
203 | else{ | |
204 | countPos++; | |
205 | sindexF[countPos+n] = val; | |
206 | sindexF[countPos]++; | |
207 | last =val; | |
208 | } | |
209 | } | |
210 | if (last==val) countPos++; | |
211 | // sort according frequency | |
212 | TMath::Sort(countPos, sindexF, sindexS, kTRUE); | |
213 | for (Int_t i=0;i<countPos;i++){ | |
214 | outlist[2*i ] = sindexF[sindexS[i]+n]; | |
215 | outlist[2*i+1] = sindexF[sindexS[i]]; | |
216 | } | |
217 | delete [] sindexS; | |
218 | delete [] sindexF; | |
219 | ||
220 | return countPos; | |
221 | ||
222 | } | |
f6659a9d | 223 | |
224 | //___AliMathBase__________________________________________________________________________ | |
225 | void AliMathBase::TruncatedMean(TH1F * his, TVectorD *param, Float_t down, Float_t up, Bool_t verbose){ | |
226 | // | |
227 | // | |
228 | // | |
229 | Int_t nbins = his->GetNbinsX(); | |
230 | Float_t nentries = his->GetEntries(); | |
231 | Float_t sum =0; | |
232 | Float_t mean = 0; | |
233 | Float_t sigma2 = 0; | |
234 | Float_t ncumul=0; | |
235 | for (Int_t ibin=1;ibin<nbins; ibin++){ | |
236 | ncumul+= his->GetBinContent(ibin); | |
237 | Float_t fraction = Float_t(ncumul)/Float_t(nentries); | |
238 | if (fraction>down && fraction<up){ | |
239 | sum+=his->GetBinContent(ibin); | |
240 | mean+=his->GetBinCenter(ibin)*his->GetBinContent(ibin); | |
241 | sigma2+=his->GetBinCenter(ibin)*his->GetBinCenter(ibin)*his->GetBinContent(ibin); | |
242 | } | |
243 | } | |
244 | mean/=sum; | |
245 | sigma2= TMath::Sqrt(TMath::Abs(sigma2/sum-mean*mean)); | |
246 | if (param){ | |
247 | (*param)[0] = his->GetMaximum(); | |
248 | (*param)[1] = mean; | |
249 | (*param)[2] = sigma2; | |
250 | ||
251 | } | |
252 | if (verbose) printf("Mean\t%f\t Sigma2\t%f\n", mean,sigma2); | |
253 | } | |
254 | ||
255 | void AliMathBase::LTM(TH1F * his, TVectorD *param , Float_t fraction, Bool_t verbose){ | |
256 | // | |
257 | // LTM | |
258 | // | |
259 | Int_t nbins = his->GetNbinsX(); | |
260 | Int_t nentries = (Int_t)his->GetEntries(); | |
261 | Double_t *data = new Double_t[nentries]; | |
262 | Int_t npoints=0; | |
263 | for (Int_t ibin=1;ibin<nbins; ibin++){ | |
264 | Float_t entriesI = his->GetBinContent(ibin); | |
265 | Float_t xcenter= his->GetBinCenter(ibin); | |
266 | for (Int_t ic=0; ic<entriesI; ic++){ | |
267 | if (npoints<nentries){ | |
268 | data[npoints]= xcenter; | |
269 | npoints++; | |
270 | } | |
271 | } | |
272 | } | |
273 | Double_t mean, sigma; | |
274 | Int_t npoints2=TMath::Min(Int_t(fraction*Float_t(npoints)),npoints-1); | |
275 | npoints2=TMath::Max(Int_t(0.5*Float_t(npoints)),npoints2); | |
276 | AliMathBase::EvaluateUni(npoints, data, mean,sigma,npoints2); | |
277 | delete [] data; | |
278 | if (verbose) printf("Mean\t%f\t Sigma2\t%f\n", mean,sigma);if (param){ | |
279 | (*param)[0] = his->GetMaximum(); | |
280 | (*param)[1] = mean; | |
281 | (*param)[2] = sigma; | |
282 | } | |
283 | } | |
284 | ||
d9a129b4 | 285 | Double_t AliMathBase::FitGaus(TH1F* his, TVectorD *param, TMatrixD */*matrix*/, Float_t xmin, Float_t xmax, Bool_t verbose){ |
f6659a9d | 286 | // |
287 | // Fit histogram with gaussian function | |
288 | // | |
289 | // Prameters: | |
290 | // return value- chi2 - if negative ( not enough points) | |
291 | // his - input histogram | |
292 | // param - vector with parameters | |
293 | // xmin, xmax - range to fit - if xmin=xmax=0 - the full histogram range used | |
294 | // Fitting: | |
295 | // 1. Step - make logarithm | |
296 | // 2. Linear fit (parabola) - more robust - always converge | |
297 | // 3. In case of small statistic bins are averaged | |
298 | // | |
299 | static TLinearFitter fitter(3,"pol2"); | |
300 | TVectorD par(3); | |
301 | TVectorD sigma(3); | |
302 | TMatrixD mat(3,3); | |
303 | if (his->GetMaximum()<4) return -1; | |
304 | if (his->GetEntries()<12) return -1; | |
305 | if (his->GetRMS()<mat.GetTol()) return -1; | |
5608e15a | 306 | Float_t maxEstimate = his->GetEntries()*his->GetBinWidth(1)/TMath::Sqrt((TMath::TwoPi()*his->GetRMS())); |
f6659a9d | 307 | Int_t dsmooth = TMath::Nint(6./TMath::Sqrt(maxEstimate)); |
308 | ||
309 | if (maxEstimate<1) return -1; | |
310 | Int_t nbins = his->GetNbinsX(); | |
311 | Int_t npoints=0; | |
312 | // | |
313 | ||
314 | ||
315 | if (xmin>=xmax){ | |
316 | xmin = his->GetXaxis()->GetXmin(); | |
317 | xmax = his->GetXaxis()->GetXmax(); | |
318 | } | |
319 | for (Int_t iter=0; iter<2; iter++){ | |
320 | fitter.ClearPoints(); | |
321 | npoints=0; | |
5608e15a | 322 | for (Int_t ibin=1;ibin<nbins+1; ibin++){ |
f6659a9d | 323 | Int_t countB=1; |
324 | Float_t entriesI = his->GetBinContent(ibin); | |
325 | for (Int_t delta = -dsmooth; delta<=dsmooth; delta++){ | |
326 | if (ibin+delta>1 &&ibin+delta<nbins-1){ | |
327 | entriesI += his->GetBinContent(ibin+delta); | |
328 | countB++; | |
329 | } | |
330 | } | |
331 | entriesI/=countB; | |
332 | Double_t xcenter= his->GetBinCenter(ibin); | |
333 | if (xcenter<xmin || xcenter>xmax) continue; | |
334 | Double_t error=1./TMath::Sqrt(countB); | |
335 | Float_t cont=2; | |
336 | if (iter>0){ | |
337 | if (par[0]+par[1]*xcenter+par[2]*xcenter*xcenter>20) return 0; | |
338 | cont = TMath::Exp(par[0]+par[1]*xcenter+par[2]*xcenter*xcenter); | |
339 | if (cont>1.) error = 1./TMath::Sqrt(cont*Float_t(countB)); | |
340 | } | |
341 | if (entriesI>1&&cont>1){ | |
342 | fitter.AddPoint(&xcenter,TMath::Log(Float_t(entriesI)),error); | |
343 | npoints++; | |
344 | } | |
345 | } | |
346 | if (npoints>3){ | |
347 | fitter.Eval(); | |
348 | fitter.GetParameters(par); | |
349 | }else{ | |
350 | break; | |
351 | } | |
352 | } | |
353 | if (npoints<=3){ | |
354 | return -1; | |
355 | } | |
356 | fitter.GetParameters(par); | |
357 | fitter.GetCovarianceMatrix(mat); | |
358 | if (TMath::Abs(par[1])<mat.GetTol()) return -1; | |
359 | if (TMath::Abs(par[2])<mat.GetTol()) return -1; | |
360 | Double_t chi2 = fitter.GetChisquare()/Float_t(npoints); | |
361 | //fitter.GetParameters(); | |
362 | if (!param) param = new TVectorD(3); | |
d9a129b4 | 363 | //if (!matrix) matrix = new TMatrixD(3,3); |
f6659a9d | 364 | (*param)[1] = par[1]/(-2.*par[2]); |
365 | (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2])); | |
366 | (*param)[0] = TMath::Exp(par[0]+ par[1]* (*param)[1] + par[2]*(*param)[1]*(*param)[1]); | |
367 | if (verbose){ | |
368 | par.Print(); | |
369 | mat.Print(); | |
370 | param->Print(); | |
371 | printf("Chi2=%f\n",chi2); | |
372 | TF1 * f1= new TF1("f1","[0]*exp(-(x-[1])^2/(2*[2]*[2]))",his->GetXaxis()->GetXmin(),his->GetXaxis()->GetXmax()); | |
373 | f1->SetParameter(0, (*param)[0]); | |
374 | f1->SetParameter(1, (*param)[1]); | |
375 | f1->SetParameter(2, (*param)[2]); | |
376 | f1->Draw("same"); | |
377 | } | |
378 | return chi2; | |
379 | } | |
380 | ||
d9a129b4 | 381 | Double_t AliMathBase::FitGaus(Float_t *arr, Int_t nBins, Float_t xMin, Float_t xMax, TVectorD *param, TMatrixD */*matrix*/, Bool_t verbose){ |
5608e15a | 382 | // |
383 | // Fit histogram with gaussian function | |
384 | // | |
385 | // Prameters: | |
5f645a6e | 386 | // nbins: size of the array and number of histogram bins |
387 | // xMin, xMax: histogram range | |
00bb7de0 | 388 | // param: paramters of the fit (0-Constant, 1-Mean, 2-Sigma, 3-Sum) |
5f645a6e | 389 | // matrix: covariance matrix -- not implemented yet, pass dummy matrix!!! |
390 | // | |
391 | // Return values: | |
392 | // >0: the chi2 returned by TLinearFitter | |
393 | // -3: only three points have been used for the calculation - no fitter was used | |
394 | // -2: only two points have been used for the calculation - center of gravity was uesed for calculation | |
395 | // -1: only one point has been used for the calculation - center of gravity was uesed for calculation | |
396 | // -4: invalid result!! | |
397 | // | |
5608e15a | 398 | // Fitting: |
399 | // 1. Step - make logarithm | |
400 | // 2. Linear fit (parabola) - more robust - always converge | |
5608e15a | 401 | // |
402 | static TLinearFitter fitter(3,"pol2"); | |
403 | static TMatrixD mat(3,3); | |
404 | static Double_t kTol = mat.GetTol(); | |
405 | fitter.StoreData(kFALSE); | |
406 | fitter.ClearPoints(); | |
407 | TVectorD par(3); | |
408 | TVectorD sigma(3); | |
409 | TMatrixD A(3,3); | |
410 | TMatrixD b(3,1); | |
5f645a6e | 411 | Float_t rms = TMath::RMS(nBins,arr); |
412 | Float_t max = TMath::MaxElement(nBins,arr); | |
413 | Float_t binWidth = (xMax-xMin)/(Float_t)nBins; | |
5608e15a | 414 | |
415 | Float_t meanCOG = 0; | |
416 | Float_t rms2COG = 0; | |
417 | Float_t sumCOG = 0; | |
418 | ||
419 | Float_t entries = 0; | |
420 | Int_t nfilled=0; | |
421 | ||
d9a129b4 | 422 | if (!param) param = new TVectorD(4); |
61d38705 | 423 | |
5f645a6e | 424 | for (Int_t i=0; i<nBins; i++){ |
5608e15a | 425 | entries+=arr[i]; |
426 | if (arr[i]>0) nfilled++; | |
427 | } | |
00bb7de0 | 428 | (*param)[0] = 0; |
429 | (*param)[1] = 0; | |
430 | (*param)[2] = 0; | |
431 | (*param)[3] = 0; | |
5608e15a | 432 | |
5f645a6e | 433 | if (max<4) return -4; |
434 | if (entries<12) return -4; | |
435 | if (rms<kTol) return -4; | |
5608e15a | 436 | |
00bb7de0 | 437 | (*param)[3] = entries; |
5608e15a | 438 | |
00bb7de0 | 439 | Int_t npoints=0; |
5f645a6e | 440 | for (Int_t ibin=0;ibin<nBins; ibin++){ |
441 | Float_t entriesI = arr[ibin]; | |
5608e15a | 442 | if (entriesI>1){ |
5f645a6e | 443 | Double_t xcenter = xMin+(ibin+0.5)*binWidth; |
5608e15a | 444 | Float_t error = 1./TMath::Sqrt(entriesI); |
445 | Float_t val = TMath::Log(Float_t(entriesI)); | |
446 | fitter.AddPoint(&xcenter,val,error); | |
5f645a6e | 447 | if (npoints<3){ |
448 | A(npoints,0)=1; | |
449 | A(npoints,1)=xcenter; | |
450 | A(npoints,2)=xcenter*xcenter; | |
451 | b(npoints,0)=val; | |
452 | meanCOG+=xcenter*entriesI; | |
453 | rms2COG +=xcenter*entriesI*xcenter; | |
454 | sumCOG +=entriesI; | |
455 | } | |
5608e15a | 456 | npoints++; |
457 | } | |
458 | } | |
5608e15a | 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 | } | |
5f645a6e | 478 | if (TMath::Abs(par[1])<kTol) return -4; |
479 | if (TMath::Abs(par[2])<kTol) return -4; | |
5608e15a | 480 | |
d9a129b4 | 481 | //if (!param) param = new TVectorD(4); |
00bb7de0 | 482 | if ( param->GetNrows()<4 ) param->ResizeTo(4); |
d9a129b4 | 483 | //if (!matrix) matrix = new TMatrixD(3,3); // !!!!might be a memory leek. use dummy matrix pointer to call this function! |
5608e15a | 484 | |
485 | (*param)[1] = par[1]/(-2.*par[2]); | |
486 | (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2])); | |
5f645a6e | 487 | Double_t lnparam0 = par[0]+ par[1]* (*param)[1] + par[2]*(*param)[1]*(*param)[1]; |
488 | if ( lnparam0>307 ) return -4; | |
489 | (*param)[0] = TMath::Exp(lnparam0); | |
5608e15a | 490 | if (verbose){ |
491 | par.Print(); | |
492 | mat.Print(); | |
493 | param->Print(); | |
494 | printf("Chi2=%f\n",chi2); | |
5f645a6e | 495 | TF1 * f1= new TF1("f1","[0]*exp(-(x-[1])^2/(2*[2]*[2]))",xMin,xMax); |
5608e15a | 496 | f1->SetParameter(0, (*param)[0]); |
497 | f1->SetParameter(1, (*param)[1]); | |
498 | f1->SetParameter(2, (*param)[2]); | |
499 | f1->Draw("same"); | |
500 | } | |
501 | return chi2; | |
502 | } | |
503 | ||
504 | if (npoints == 2){ | |
505 | //use center of gravity for 2 points | |
506 | meanCOG/=sumCOG; | |
507 | rms2COG /=sumCOG; | |
508 | (*param)[0] = max; | |
509 | (*param)[1] = meanCOG; | |
510 | (*param)[2] = TMath::Sqrt(TMath::Abs(meanCOG*meanCOG-rms2COG)); | |
511 | chi2=-2.; | |
512 | } | |
513 | if ( npoints == 1 ){ | |
5f645a6e | 514 | meanCOG/=sumCOG; |
5608e15a | 515 | (*param)[0] = max; |
516 | (*param)[1] = meanCOG; | |
517 | (*param)[2] = binWidth/TMath::Sqrt(12); | |
518 | chi2=-1.; | |
519 | } | |
520 | return chi2; | |
521 | ||
522 | } | |
523 | ||
524 | ||
5f645a6e | 525 | Float_t AliMathBase::GetCOG(Short_t *arr, Int_t nBins, Float_t xMin, Float_t xMax, Float_t *rms, Float_t *sum) |
526 | { | |
527 | // | |
528 | // calculate center of gravity rms and sum for array 'arr' with nBins an a x range xMin to xMax | |
529 | // return COG; in case of failure return xMin | |
530 | // | |
531 | Float_t meanCOG = 0; | |
532 | Float_t rms2COG = 0; | |
533 | Float_t sumCOG = 0; | |
534 | Int_t npoints = 0; | |
535 | ||
536 | Float_t binWidth = (xMax-xMin)/(Float_t)nBins; | |
537 | ||
538 | for (Int_t ibin=0; ibin<nBins; ibin++){ | |
539 | Float_t entriesI = (Float_t)arr[ibin]; | |
540 | Double_t xcenter = xMin+(ibin+0.5)*binWidth; | |
541 | if ( entriesI>0 ){ | |
542 | meanCOG += xcenter*entriesI; | |
543 | rms2COG += xcenter*entriesI*xcenter; | |
544 | sumCOG += entriesI; | |
545 | npoints++; | |
546 | } | |
547 | } | |
548 | if ( sumCOG == 0 ) return xMin; | |
549 | meanCOG/=sumCOG; | |
550 | ||
551 | if ( rms ){ | |
552 | rms2COG /=sumCOG; | |
553 | (*rms) = TMath::Sqrt(TMath::Abs(meanCOG*meanCOG-rms2COG)); | |
554 | if ( npoints == 1 ) (*rms) = binWidth/TMath::Sqrt(12); | |
555 | } | |
556 | ||
557 | if ( sum ) | |
558 | (*sum) = sumCOG; | |
559 | ||
560 | return meanCOG; | |
561 | } | |
562 | ||
563 | ||
bb7e41dd | 564 | Double_t AliMathBase::ErfcFast(Double_t x){ |
565 | // Fast implementation of the complementary error function | |
566 | // The error of the approximation is |eps(x)| < 5E-4 | |
567 | // See Abramowitz and Stegun, p.299, 7.1.27 | |
568 | ||
569 | Double_t z = TMath::Abs(x); | |
570 | Double_t ans = 1+z*(0.278393+z*(0.230389+z*(0.000972+z*0.078108))); | |
571 | ans = 1.0/ans; | |
572 | ans *= ans; | |
573 | ans *= ans; | |
574 | ||
575 | return (x>=0.0 ? ans : 2.0 - ans); | |
576 | } | |
5608e15a | 577 | |
578 | /////////////////////////////////////////////////////////////// | |
579 | ////////////// TEST functions ///////////////////////// | |
580 | /////////////////////////////////////////////////////////////// | |
581 | ||
582 | ||
583 | ||
584 | ||
585 | ||
586 | void AliMathBase::TestGausFit(Int_t nhistos){ | |
587 | // | |
588 | // Test performance of the parabolic - gaussian fit - compare it with | |
589 | // ROOT gauss fit | |
590 | // nhistos - number of histograms to be used for test | |
591 | // | |
592 | TTreeSRedirector *pcstream = new TTreeSRedirector("fitdebug.root"); | |
593 | ||
594 | Float_t *xTrue = new Float_t[nhistos]; | |
595 | Float_t *sTrue = new Float_t[nhistos]; | |
596 | TVectorD **par1 = new TVectorD*[nhistos]; | |
597 | TVectorD **par2 = new TVectorD*[nhistos]; | |
598 | TMatrixD dummy(3,3); | |
599 | ||
600 | ||
601 | TH1F **h1f = new TH1F*[nhistos]; | |
602 | TF1 *myg = new TF1("myg","gaus"); | |
603 | TF1 *fit = new TF1("fit","gaus"); | |
604 | gRandom->SetSeed(0); | |
605 | ||
606 | //init | |
607 | for (Int_t i=0;i<nhistos; i++){ | |
608 | par1[i] = new TVectorD(3); | |
609 | par2[i] = new TVectorD(3); | |
610 | h1f[i] = new TH1F(Form("h1f%d",i),Form("h1f%d",i),20,-10,10); | |
611 | xTrue[i]= gRandom->Rndm(); | |
612 | gSystem->Sleep(2); | |
613 | sTrue[i]= .75+gRandom->Rndm()*.5; | |
614 | myg->SetParameters(1,xTrue[i],sTrue[i]); | |
615 | h1f[i]->FillRandom("myg"); | |
616 | } | |
617 | ||
618 | TStopwatch s; | |
619 | s.Start(); | |
620 | //standard gaus fit | |
621 | for (Int_t i=0; i<nhistos; i++){ | |
622 | h1f[i]->Fit(fit,"0q"); | |
623 | (*par1[i])(0) = fit->GetParameter(0); | |
624 | (*par1[i])(1) = fit->GetParameter(1); | |
625 | (*par1[i])(2) = fit->GetParameter(2); | |
626 | } | |
627 | s.Stop(); | |
628 | printf("Gaussian fit\t"); | |
629 | s.Print(); | |
630 | ||
631 | s.Start(); | |
632 | //AliMathBase gaus fit | |
633 | for (Int_t i=0; i<nhistos; i++){ | |
5f645a6e | 634 | AliMathBase::FitGaus(h1f[i]->GetArray()+1,h1f[i]->GetNbinsX(),h1f[i]->GetXaxis()->GetXmin(),h1f[i]->GetXaxis()->GetXmax(),par2[i],&dummy); |
5608e15a | 635 | } |
636 | ||
637 | s.Stop(); | |
638 | printf("Parabolic fit\t"); | |
639 | s.Print(); | |
640 | //write stream | |
641 | for (Int_t i=0;i<nhistos; i++){ | |
642 | Float_t xt = xTrue[i]; | |
643 | Float_t st = sTrue[i]; | |
644 | (*pcstream)<<"data" | |
645 | <<"xTrue="<<xt | |
646 | <<"sTrue="<<st | |
647 | <<"pg.="<<(par1[i]) | |
648 | <<"pa.="<<(par2[i]) | |
649 | <<"\n"; | |
650 | } | |
651 | //delete pointers | |
652 | for (Int_t i=0;i<nhistos; i++){ | |
653 | delete par1[i]; | |
654 | delete par2[i]; | |
655 | delete h1f[i]; | |
656 | } | |
657 | delete pcstream; | |
658 | delete []h1f; | |
659 | delete []xTrue; | |
660 | delete []sTrue; | |
661 | // | |
662 | delete []par1; | |
663 | delete []par2; | |
664 | ||
665 | } | |
666 | ||
667 | ||
668 | ||
3392b4c9 | 669 | TGraph2D * AliMathBase::MakeStat2D(TH3 * his, Int_t delta0, Int_t delta1, Int_t type){ |
670 | // | |
671 | // | |
672 | // | |
673 | // delta - number of bins to integrate | |
674 | // type - 0 - mean value | |
675 | ||
676 | TAxis * xaxis = his->GetXaxis(); | |
677 | TAxis * yaxis = his->GetYaxis(); | |
678 | // TAxis * zaxis = his->GetZaxis(); | |
679 | Int_t nbinx = xaxis->GetNbins(); | |
680 | Int_t nbiny = yaxis->GetNbins(); | |
61d38705 | 681 | const Int_t nc=1000; |
682 | char name[nc]; | |
3392b4c9 | 683 | Int_t icount=0; |
684 | TGraph2D *graph = new TGraph2D(nbinx*nbiny); | |
685 | TF1 f1("f1","gaus"); | |
686 | for (Int_t ix=0; ix<nbinx;ix++) | |
687 | for (Int_t iy=0; iy<nbiny;iy++){ | |
688 | Float_t xcenter = xaxis->GetBinCenter(ix); | |
689 | Float_t ycenter = yaxis->GetBinCenter(iy); | |
61d38705 | 690 | snprintf(name,nc,"%s_%d_%d",his->GetName(), ix,iy); |
3392b4c9 | 691 | TH1 *projection = his->ProjectionZ(name,ix-delta0,ix+delta0,iy-delta1,iy+delta1); |
692 | Float_t stat= 0; | |
693 | if (type==0) stat = projection->GetMean(); | |
694 | if (type==1) stat = projection->GetRMS(); | |
695 | if (type==2 || type==3){ | |
696 | TVectorD vec(3); | |
697 | AliMathBase::LTM((TH1F*)projection,&vec,0.7); | |
698 | if (type==2) stat= vec[1]; | |
699 | if (type==3) stat= vec[0]; | |
700 | } | |
701 | if (type==4|| type==5){ | |
702 | projection->Fit(&f1); | |
703 | if (type==4) stat= f1.GetParameter(1); | |
704 | if (type==5) stat= f1.GetParameter(2); | |
705 | } | |
706 | //printf("%d\t%f\t%f\t%f\n", icount,xcenter, ycenter, stat); | |
707 | graph->SetPoint(icount,xcenter, ycenter, stat); | |
708 | icount++; | |
709 | } | |
710 | return graph; | |
711 | } | |
5608e15a | 712 | |
3392b4c9 | 713 | TGraph * AliMathBase::MakeStat1D(TH3 * his, Int_t delta1, Int_t type){ |
714 | // | |
715 | // | |
716 | // | |
717 | // delta - number of bins to integrate | |
718 | // type - 0 - mean value | |
719 | ||
720 | TAxis * xaxis = his->GetXaxis(); | |
721 | TAxis * yaxis = his->GetYaxis(); | |
722 | // TAxis * zaxis = his->GetZaxis(); | |
723 | Int_t nbinx = xaxis->GetNbins(); | |
724 | Int_t nbiny = yaxis->GetNbins(); | |
61d38705 | 725 | const Int_t nc=1000; |
726 | char name[nc]; | |
3392b4c9 | 727 | Int_t icount=0; |
728 | TGraph *graph = new TGraph(nbinx); | |
729 | TF1 f1("f1","gaus"); | |
730 | for (Int_t ix=0; ix<nbinx;ix++){ | |
731 | Float_t xcenter = xaxis->GetBinCenter(ix); | |
732 | // Float_t ycenter = yaxis->GetBinCenter(iy); | |
61d38705 | 733 | snprintf(name,nc,"%s_%d",his->GetName(), ix); |
3392b4c9 | 734 | TH1 *projection = his->ProjectionZ(name,ix-delta1,ix+delta1,0,nbiny); |
735 | Float_t stat= 0; | |
736 | if (type==0) stat = projection->GetMean(); | |
737 | if (type==1) stat = projection->GetRMS(); | |
738 | if (type==2 || type==3){ | |
739 | TVectorD vec(3); | |
740 | AliMathBase::LTM((TH1F*)projection,&vec,0.7); | |
741 | if (type==2) stat= vec[1]; | |
742 | if (type==3) stat= vec[0]; | |
743 | } | |
744 | if (type==4|| type==5){ | |
745 | projection->Fit(&f1); | |
746 | if (type==4) stat= f1.GetParameter(1); | |
747 | if (type==5) stat= f1.GetParameter(2); | |
748 | } | |
749 | //printf("%d\t%f\t%f\t%f\n", icount,xcenter, ycenter, stat); | |
750 | graph->SetPoint(icount,xcenter, stat); | |
751 | icount++; | |
752 | } | |
753 | return graph; | |
754 | } | |
09a29b67 | 755 | |
756 | Double_t AliMathBase::TruncatedGaus(Double_t mean, Double_t sigma, Double_t cutat) | |
757 | { | |
758 | // return number generated according to a gaussian distribution N(mean,sigma) truncated at cutat | |
759 | // | |
760 | Double_t value; | |
761 | do{ | |
762 | value=gRandom->Gaus(mean,sigma); | |
763 | }while(TMath::Abs(value-mean)>cutat); | |
764 | return value; | |
765 | } | |
4c5881a0 | 766 | |
767 | Double_t AliMathBase::TruncatedGaus(Double_t mean, Double_t sigma, Double_t leftCut, Double_t rightCut) | |
768 | { | |
769 | // return number generated according to a gaussian distribution N(mean,sigma) | |
770 | // truncated at leftCut and rightCut | |
771 | // | |
772 | Double_t value; | |
773 | do{ | |
774 | value=gRandom->Gaus(mean,sigma); | |
775 | }while((value-mean)<-leftCut || (value-mean)>rightCut); | |
776 | return value; | |
777 | } | |
f6d87824 | 778 | |
779 | Double_t AliMathBase::BetheBlochAleph(Double_t bg, | |
780 | Double_t kp1, | |
781 | Double_t kp2, | |
782 | Double_t kp3, | |
783 | Double_t kp4, | |
784 | Double_t kp5) { | |
785 | // | |
d46683db | 786 | // This is the empirical ALEPH parameterization of the Bethe-Bloch formula. |
787 | // It is normalized to 1 at the minimum. | |
788 | // | |
789 | // bg - beta*gamma | |
790 | // | |
791 | // The default values for the kp* parameters are for ALICE TPC. | |
792 | // The returned value is in MIP units | |
f6d87824 | 793 | // |
794 | ||
9c56b409 | 795 | return AliExternalTrackParam::BetheBlochAleph(bg,kp1,kp2,kp3,kp4,kp5); |
d46683db | 796 | } |