]>
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" |
28 | #include "TF1.h" | |
29 | #include "TLinearFitter.h" | |
284050f7 | 30 | |
31 | ClassImp(AliMathBase) // Class implementation to enable ROOT I/O | |
32 | ||
33 | AliMathBase::AliMathBase() : TObject() | |
34 | { | |
35 | // Default constructor | |
36 | } | |
37 | /////////////////////////////////////////////////////////////////////////// | |
38 | AliMathBase::~AliMathBase() | |
39 | { | |
40 | // Destructor | |
41 | } | |
42 | ||
43 | ||
44 | //_____________________________________________________________________________ | |
45 | void AliMathBase::EvaluateUni(Int_t nvectors, Double_t *data, Double_t &mean | |
46 | , Double_t &sigma, Int_t hh) | |
47 | { | |
48 | // | |
49 | // Robust estimator in 1D case MI version - (faster than ROOT version) | |
50 | // | |
51 | // For the univariate case | |
52 | // estimates of location and scatter are returned in mean and sigma parameters | |
53 | // the algorithm works on the same principle as in multivariate case - | |
54 | // it finds a subset of size hh with smallest sigma, and then returns mean and | |
55 | // sigma of this subset | |
56 | // | |
57 | ||
58 | if (hh==0) | |
59 | hh=(nvectors+2)/2; | |
60 | 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}; | |
61 | Int_t *index=new Int_t[nvectors]; | |
62 | TMath::Sort(nvectors, data, index, kFALSE); | |
63 | ||
64 | Int_t nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11); | |
d9e9045c | 65 | Double_t factor = faclts[TMath::Max(0,nquant-1)]; |
284050f7 | 66 | |
67 | Double_t sumx =0; | |
68 | Double_t sumx2 =0; | |
69 | Int_t bestindex = -1; | |
70 | Double_t bestmean = 0; | |
71 | Double_t bestsigma = data[index[nvectors-1]]-data[index[0]]; // maximal possible sigma | |
72 | for (Int_t i=0; i<hh; i++){ | |
73 | sumx += data[index[i]]; | |
74 | sumx2 += data[index[i]]*data[index[i]]; | |
75 | } | |
76 | ||
77 | Double_t norm = 1./Double_t(hh); | |
78 | Double_t norm2 = 1./Double_t(hh-1); | |
79 | for (Int_t i=hh; i<nvectors; i++){ | |
80 | Double_t cmean = sumx*norm; | |
81 | Double_t csigma = (sumx2 - hh*cmean*cmean)*norm2; | |
82 | if (csigma<bestsigma){ | |
83 | bestmean = cmean; | |
84 | bestsigma = csigma; | |
85 | bestindex = i-hh; | |
86 | } | |
87 | ||
88 | sumx += data[index[i]]-data[index[i-hh]]; | |
89 | sumx2 += data[index[i]]*data[index[i]]-data[index[i-hh]]*data[index[i-hh]]; | |
90 | } | |
91 | ||
92 | Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma)); | |
93 | mean = bestmean; | |
94 | sigma = bstd; | |
95 | delete [] index; | |
96 | ||
97 | } | |
98 | ||
99 | ||
100 | ||
101 | void AliMathBase::EvaluateUniExternal(Int_t nvectors, Double_t *data, Double_t &mean, Double_t &sigma, Int_t hh, Float_t externalfactor) | |
102 | { | |
103 | // Modified version of ROOT robust EvaluateUni | |
104 | // robust estimator in 1D case MI version | |
105 | // added external factor to include precision of external measurement | |
106 | // | |
107 | ||
108 | if (hh==0) | |
109 | hh=(nvectors+2)/2; | |
110 | 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}; | |
111 | Int_t *index=new Int_t[nvectors]; | |
112 | TMath::Sort(nvectors, data, index, kFALSE); | |
113 | // | |
114 | Int_t nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11); | |
115 | Double_t factor = faclts[0]; | |
116 | if (nquant>0){ | |
117 | // fix proper normalization - Anja | |
118 | factor = faclts[nquant-1]; | |
119 | } | |
120 | ||
121 | // | |
122 | // | |
123 | Double_t sumx =0; | |
124 | Double_t sumx2 =0; | |
125 | Int_t bestindex = -1; | |
126 | Double_t bestmean = 0; | |
127 | Double_t bestsigma = -1; | |
128 | for (Int_t i=0; i<hh; i++){ | |
129 | sumx += data[index[i]]; | |
130 | sumx2 += data[index[i]]*data[index[i]]; | |
131 | } | |
132 | // | |
133 | Double_t kfactor = 2.*externalfactor - externalfactor*externalfactor; | |
134 | Double_t norm = 1./Double_t(hh); | |
135 | for (Int_t i=hh; i<nvectors; i++){ | |
136 | Double_t cmean = sumx*norm; | |
137 | Double_t csigma = (sumx2*norm - cmean*cmean*kfactor); | |
138 | if (csigma<bestsigma || bestsigma<0){ | |
139 | bestmean = cmean; | |
140 | bestsigma = csigma; | |
141 | bestindex = i-hh; | |
142 | } | |
143 | // | |
144 | // | |
145 | sumx += data[index[i]]-data[index[i-hh]]; | |
146 | sumx2 += data[index[i]]*data[index[i]]-data[index[i-hh]]*data[index[i-hh]]; | |
147 | } | |
148 | ||
149 | Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma)); | |
150 | mean = bestmean; | |
151 | sigma = bstd; | |
152 | delete [] index; | |
153 | } | |
154 | ||
155 | ||
156 | //_____________________________________________________________________________ | |
157 | Int_t AliMathBase::Freq(Int_t n, const Int_t *inlist | |
158 | , Int_t *outlist, Bool_t down) | |
159 | { | |
160 | // | |
161 | // Sort eleements according occurancy | |
162 | // The size of output array has is 2*n | |
163 | // | |
164 | ||
165 | Int_t * sindexS = new Int_t[n]; // temp array for sorting | |
166 | Int_t * sindexF = new Int_t[2*n]; | |
167 | for (Int_t i=0;i<n;i++) sindexF[i]=0; | |
168 | // | |
169 | TMath::Sort(n,inlist, sindexS, down); | |
170 | Int_t last = inlist[sindexS[0]]; | |
171 | Int_t val = last; | |
172 | sindexF[0] = 1; | |
173 | sindexF[0+n] = last; | |
174 | Int_t countPos = 0; | |
175 | // | |
176 | // find frequency | |
177 | for(Int_t i=1;i<n; i++){ | |
178 | val = inlist[sindexS[i]]; | |
179 | if (last == val) sindexF[countPos]++; | |
180 | else{ | |
181 | countPos++; | |
182 | sindexF[countPos+n] = val; | |
183 | sindexF[countPos]++; | |
184 | last =val; | |
185 | } | |
186 | } | |
187 | if (last==val) countPos++; | |
188 | // sort according frequency | |
189 | TMath::Sort(countPos, sindexF, sindexS, kTRUE); | |
190 | for (Int_t i=0;i<countPos;i++){ | |
191 | outlist[2*i ] = sindexF[sindexS[i]+n]; | |
192 | outlist[2*i+1] = sindexF[sindexS[i]]; | |
193 | } | |
194 | delete [] sindexS; | |
195 | delete [] sindexF; | |
196 | ||
197 | return countPos; | |
198 | ||
199 | } | |
f6659a9d | 200 | |
201 | //___AliMathBase__________________________________________________________________________ | |
202 | void AliMathBase::TruncatedMean(TH1F * his, TVectorD *param, Float_t down, Float_t up, Bool_t verbose){ | |
203 | // | |
204 | // | |
205 | // | |
206 | Int_t nbins = his->GetNbinsX(); | |
207 | Float_t nentries = his->GetEntries(); | |
208 | Float_t sum =0; | |
209 | Float_t mean = 0; | |
210 | Float_t sigma2 = 0; | |
211 | Float_t ncumul=0; | |
212 | for (Int_t ibin=1;ibin<nbins; ibin++){ | |
213 | ncumul+= his->GetBinContent(ibin); | |
214 | Float_t fraction = Float_t(ncumul)/Float_t(nentries); | |
215 | if (fraction>down && fraction<up){ | |
216 | sum+=his->GetBinContent(ibin); | |
217 | mean+=his->GetBinCenter(ibin)*his->GetBinContent(ibin); | |
218 | sigma2+=his->GetBinCenter(ibin)*his->GetBinCenter(ibin)*his->GetBinContent(ibin); | |
219 | } | |
220 | } | |
221 | mean/=sum; | |
222 | sigma2= TMath::Sqrt(TMath::Abs(sigma2/sum-mean*mean)); | |
223 | if (param){ | |
224 | (*param)[0] = his->GetMaximum(); | |
225 | (*param)[1] = mean; | |
226 | (*param)[2] = sigma2; | |
227 | ||
228 | } | |
229 | if (verbose) printf("Mean\t%f\t Sigma2\t%f\n", mean,sigma2); | |
230 | } | |
231 | ||
232 | void AliMathBase::LTM(TH1F * his, TVectorD *param , Float_t fraction, Bool_t verbose){ | |
233 | // | |
234 | // LTM | |
235 | // | |
236 | Int_t nbins = his->GetNbinsX(); | |
237 | Int_t nentries = (Int_t)his->GetEntries(); | |
238 | Double_t *data = new Double_t[nentries]; | |
239 | Int_t npoints=0; | |
240 | for (Int_t ibin=1;ibin<nbins; ibin++){ | |
241 | Float_t entriesI = his->GetBinContent(ibin); | |
242 | Float_t xcenter= his->GetBinCenter(ibin); | |
243 | for (Int_t ic=0; ic<entriesI; ic++){ | |
244 | if (npoints<nentries){ | |
245 | data[npoints]= xcenter; | |
246 | npoints++; | |
247 | } | |
248 | } | |
249 | } | |
250 | Double_t mean, sigma; | |
251 | Int_t npoints2=TMath::Min(Int_t(fraction*Float_t(npoints)),npoints-1); | |
252 | npoints2=TMath::Max(Int_t(0.5*Float_t(npoints)),npoints2); | |
253 | AliMathBase::EvaluateUni(npoints, data, mean,sigma,npoints2); | |
254 | delete [] data; | |
255 | if (verbose) printf("Mean\t%f\t Sigma2\t%f\n", mean,sigma);if (param){ | |
256 | (*param)[0] = his->GetMaximum(); | |
257 | (*param)[1] = mean; | |
258 | (*param)[2] = sigma; | |
259 | } | |
260 | } | |
261 | ||
262 | Double_t AliMathBase::FitGaus(TH1F* his, TVectorD *param, TMatrixD *matrix, Float_t xmin, Float_t xmax, Bool_t verbose){ | |
263 | // | |
264 | // Fit histogram with gaussian function | |
265 | // | |
266 | // Prameters: | |
267 | // return value- chi2 - if negative ( not enough points) | |
268 | // his - input histogram | |
269 | // param - vector with parameters | |
270 | // xmin, xmax - range to fit - if xmin=xmax=0 - the full histogram range used | |
271 | // Fitting: | |
272 | // 1. Step - make logarithm | |
273 | // 2. Linear fit (parabola) - more robust - always converge | |
274 | // 3. In case of small statistic bins are averaged | |
275 | // | |
276 | static TLinearFitter fitter(3,"pol2"); | |
277 | TVectorD par(3); | |
278 | TVectorD sigma(3); | |
279 | TMatrixD mat(3,3); | |
280 | if (his->GetMaximum()<4) return -1; | |
281 | if (his->GetEntries()<12) return -1; | |
282 | if (his->GetRMS()<mat.GetTol()) return -1; | |
283 | Float_t maxEstimate = his->GetEntries()*his->GetBinWidth(1)/TMath::Sqrt((2.*TMath::Pi()*his->GetRMS())); | |
284 | Int_t dsmooth = TMath::Nint(6./TMath::Sqrt(maxEstimate)); | |
285 | ||
286 | if (maxEstimate<1) return -1; | |
287 | Int_t nbins = his->GetNbinsX(); | |
288 | Int_t npoints=0; | |
289 | // | |
290 | ||
291 | ||
292 | if (xmin>=xmax){ | |
293 | xmin = his->GetXaxis()->GetXmin(); | |
294 | xmax = his->GetXaxis()->GetXmax(); | |
295 | } | |
296 | for (Int_t iter=0; iter<2; iter++){ | |
297 | fitter.ClearPoints(); | |
298 | npoints=0; | |
299 | for (Int_t ibin=1;ibin<nbins-1; ibin++){ | |
300 | Int_t countB=1; | |
301 | Float_t entriesI = his->GetBinContent(ibin); | |
302 | for (Int_t delta = -dsmooth; delta<=dsmooth; delta++){ | |
303 | if (ibin+delta>1 &&ibin+delta<nbins-1){ | |
304 | entriesI += his->GetBinContent(ibin+delta); | |
305 | countB++; | |
306 | } | |
307 | } | |
308 | entriesI/=countB; | |
309 | Double_t xcenter= his->GetBinCenter(ibin); | |
310 | if (xcenter<xmin || xcenter>xmax) continue; | |
311 | Double_t error=1./TMath::Sqrt(countB); | |
312 | Float_t cont=2; | |
313 | if (iter>0){ | |
314 | if (par[0]+par[1]*xcenter+par[2]*xcenter*xcenter>20) return 0; | |
315 | cont = TMath::Exp(par[0]+par[1]*xcenter+par[2]*xcenter*xcenter); | |
316 | if (cont>1.) error = 1./TMath::Sqrt(cont*Float_t(countB)); | |
317 | } | |
318 | if (entriesI>1&&cont>1){ | |
319 | fitter.AddPoint(&xcenter,TMath::Log(Float_t(entriesI)),error); | |
320 | npoints++; | |
321 | } | |
322 | } | |
323 | if (npoints>3){ | |
324 | fitter.Eval(); | |
325 | fitter.GetParameters(par); | |
326 | }else{ | |
327 | break; | |
328 | } | |
329 | } | |
330 | if (npoints<=3){ | |
331 | return -1; | |
332 | } | |
333 | fitter.GetParameters(par); | |
334 | fitter.GetCovarianceMatrix(mat); | |
335 | if (TMath::Abs(par[1])<mat.GetTol()) return -1; | |
336 | if (TMath::Abs(par[2])<mat.GetTol()) return -1; | |
337 | Double_t chi2 = fitter.GetChisquare()/Float_t(npoints); | |
338 | //fitter.GetParameters(); | |
339 | if (!param) param = new TVectorD(3); | |
340 | if (!matrix) matrix = new TMatrixD(3,3); | |
341 | (*param)[1] = par[1]/(-2.*par[2]); | |
342 | (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2])); | |
343 | (*param)[0] = TMath::Exp(par[0]+ par[1]* (*param)[1] + par[2]*(*param)[1]*(*param)[1]); | |
344 | if (verbose){ | |
345 | par.Print(); | |
346 | mat.Print(); | |
347 | param->Print(); | |
348 | printf("Chi2=%f\n",chi2); | |
349 | TF1 * f1= new TF1("f1","[0]*exp(-(x-[1])^2/(2*[2]*[2]))",his->GetXaxis()->GetXmin(),his->GetXaxis()->GetXmax()); | |
350 | f1->SetParameter(0, (*param)[0]); | |
351 | f1->SetParameter(1, (*param)[1]); | |
352 | f1->SetParameter(2, (*param)[2]); | |
353 | f1->Draw("same"); | |
354 | } | |
355 | return chi2; | |
356 | } | |
357 | ||
358 |