Update master to aliroot
[u/mrichter/AliRoot.git] / STEER / STEERBase / AliMathBase.cxx
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 ClassImp(AliMathBase) // Class implementation to enable ROOT I/O
45  
46 AliMathBase::AliMathBase() : TObject()
47 {
48   //
49   // Default constructor
50   //
51 }
52 ///////////////////////////////////////////////////////////////////////////
53 AliMathBase::~AliMathBase()
54 {
55   //
56   // Destructor
57   //
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
75   if (nvectors<2) {
76     AliErrorClass(Form("nvectors = %d, should be > 1",nvectors));
77     return;
78   }
79   if (hh<2)
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);
86   Double_t factor = faclts[TMath::Max(0,nquant-1)];
87   
88   Double_t sumx  =0;
89   Double_t sumx2 =0;
90   Int_t    bestindex = -1;
91   Double_t bestmean  = 0; 
92   Double_t bestsigma = (data[index[nvectors-1]]-data[index[0]]+1.);   // maximal possible sigma
93   bestsigma *=bestsigma;
94
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 }
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
285 Double_t  AliMathBase::FitGaus(TH1F* his, TVectorD *param, TMatrixD */*matrix*/, Float_t xmin, Float_t xmax, Bool_t verbose){
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;
306   Float_t maxEstimate   = his->GetEntries()*his->GetBinWidth(1)/TMath::Sqrt((TMath::TwoPi()*his->GetRMS()));
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;
322     for (Int_t ibin=1;ibin<nbins+1; ibin++){
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);
363   //if (!matrix) matrix = new TMatrixD(3,3);
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
381 Double_t  AliMathBase::FitGaus(Float_t *arr, Int_t nBins, Float_t xMin, Float_t xMax, TVectorD *param, TMatrixD */*matrix*/, Bool_t verbose){
382   //
383   //  Fit histogram with gaussian function
384   //  
385   //  Prameters:
386   //     nbins: size of the array and number of histogram bins
387   //     xMin, xMax: histogram range
388   //     param: paramters of the fit (0-Constant, 1-Mean, 2-Sigma, 3-Sum)
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   //
398   //  Fitting:
399   //  1. Step - make logarithm
400   //  2. Linear  fit (parabola) - more robust - always converge
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);
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;
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
422   if (!param)  param  = new TVectorD(4);
423
424   for (Int_t i=0; i<nBins; i++){
425       entries+=arr[i];
426       if (arr[i]>0) nfilled++;
427   }
428   (*param)[0] = 0;
429   (*param)[1] = 0;
430   (*param)[2] = 0;
431   (*param)[3] = 0;
432
433   if (max<4) return -4;
434   if (entries<12) return -4;
435   if (rms<kTol) return -4;
436
437   (*param)[3] = entries;
438
439   Int_t npoints=0;
440   for (Int_t ibin=0;ibin<nBins; ibin++){
441       Float_t entriesI = arr[ibin];
442     if (entriesI>1){
443       Double_t xcenter = xMin+(ibin+0.5)*binWidth;
444       Float_t error    = 1./TMath::Sqrt(entriesI);
445       Float_t val = TMath::Log(Float_t(entriesI));
446       fitter.AddPoint(&xcenter,val,error);
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       }
456       npoints++;
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 -4;
479       if (TMath::Abs(par[2])<kTol) return -4;
480
481       //if (!param)  param  = new TVectorD(4);
482       if ( param->GetNrows()<4 ) param->ResizeTo(4);
483       //if (!matrix) matrix = new TMatrixD(3,3);  // !!!!might be a memory leek. use dummy matrix pointer to call this function!
484
485       (*param)[1] = par[1]/(-2.*par[2]);
486       (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
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);
490       if (verbose){
491           par.Print();
492           mat.Print();
493           param->Print();
494           printf("Chi2=%f\n",chi2);
495           TF1 * f1= new TF1("f1","[0]*exp(-(x-[1])^2/(2*[2]*[2]))",xMin,xMax);
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 ){
514       meanCOG/=sumCOG;
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
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
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 }
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++){
634     AliMathBase::FitGaus(h1f[i]->GetArray()+1,h1f[i]->GetNbinsX(),h1f[i]->GetXaxis()->GetXmin(),h1f[i]->GetXaxis()->GetXmax(),par2[i],&dummy);
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
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();
681   const Int_t nc=1000;
682   char name[nc];
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); 
690       snprintf(name,nc,"%s_%d_%d",his->GetName(), ix,iy);
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 }
712
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();
725   const Int_t nc=1000;
726   char name[nc];
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); 
733     snprintf(name,nc,"%s_%d",his->GetName(), ix);
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 }
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 }
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 }
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   //
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
793   //
794
795   return AliExternalTrackParam::BetheBlochAleph(bg,kp1,kp2,kp3,kp4,kp5);
796 }
797
798 Double_t AliMathBase::Gamma(Double_t k)
799 {
800   // from
801   // Hisashi Tanizaki
802   // http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.158.3866&rep=rep1&type=pdf
803   // A. Morsch 14/01/2014
804   static Double_t n=0;
805   static Double_t c1=0;
806   static Double_t c2=0;
807   static Double_t b1=0;
808   static Double_t b2=0;
809   if (k > 0) {
810     if (k < 0.4) 
811       n = 1./k;
812     else if (k >= 0.4 && k < 4) 
813       n = 1./k + (k - 0.4)/k/3.6;
814     else if (k >= 4.) 
815       n = 1./TMath::Sqrt(k);
816     b1 = k - 1./n;
817     b2 = k + 1./n;
818     c1 = (k < 0.4)? 0 : b1 * (TMath::Log(b1) - 1.)/2.;
819     c2 = b2 * (TMath::Log(b2) - 1.)/2.;
820   }
821   Double_t x;
822   Double_t y = -1.;
823   while (1) {
824     Double_t nu1 = gRandom->Rndm();
825     Double_t nu2 = gRandom->Rndm();
826     Double_t w1 = c1 + TMath::Log(nu1);
827     Double_t w2 = c2 + TMath::Log(nu2);
828     y = n * (b1 * w2 - b2 * w1);
829     if (y < 0) continue;
830     x = n * (w2 - w1);
831     if (TMath::Log(y) >= x) break;
832   }
833   return TMath::Exp(x);
834 }