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