]>
Commit | Line | Data |
---|---|---|
21f3a443 | 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 TStatToolkit | |
19 | // | |
20 | // Subset of matheamtical functions not included in the TMath | |
21 | // | |
22 | ||
23 | /////////////////////////////////////////////////////////////////////////// | |
24 | #include "TMath.h" | |
25 | #include "Riostream.h" | |
26 | #include "TH1F.h" | |
27 | #include "TH3.h" | |
28 | #include "TF1.h" | |
29 | #include "TTree.h" | |
30 | #include "TChain.h" | |
31 | #include "TObjString.h" | |
32 | #include "TLinearFitter.h" | |
3d7cc0b4 | 33 | #include "TGraph2D.h" |
34 | #include "TGraph.h" | |
b3453fe7 | 35 | #include "TGraphErrors.h" |
21f3a443 | 36 | |
37 | // | |
38 | // includes neccessary for test functions | |
39 | // | |
40 | #include "TSystem.h" | |
41 | #include "TRandom.h" | |
42 | #include "TStopwatch.h" | |
43 | #include "TTreeStream.h" | |
44 | ||
45 | #include "TStatToolkit.h" | |
46 | ||
47 | ||
48 | ClassImp(TStatToolkit) // Class implementation to enable ROOT I/O | |
49 | ||
50 | TStatToolkit::TStatToolkit() : TObject() | |
51 | { | |
52 | // | |
53 | // Default constructor | |
54 | // | |
55 | } | |
56 | /////////////////////////////////////////////////////////////////////////// | |
57 | TStatToolkit::~TStatToolkit() | |
58 | { | |
59 | // | |
60 | // Destructor | |
61 | // | |
62 | } | |
63 | ||
64 | ||
65 | //_____________________________________________________________________________ | |
66 | void TStatToolkit::EvaluateUni(Int_t nvectors, Double_t *data, Double_t &mean | |
67 | , Double_t &sigma, Int_t hh) | |
68 | { | |
69 | // | |
70 | // Robust estimator in 1D case MI version - (faster than ROOT version) | |
71 | // | |
72 | // For the univariate case | |
73 | // estimates of location and scatter are returned in mean and sigma parameters | |
74 | // the algorithm works on the same principle as in multivariate case - | |
75 | // it finds a subset of size hh with smallest sigma, and then returns mean and | |
76 | // sigma of this subset | |
77 | // | |
78 | ||
79 | if (hh==0) | |
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); | |
bd7b4d18 | 101 | Double_t norm2 = (hh-1)>0 ? 1./Double_t(hh-1):1; |
21f3a443 | 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 TStatToolkit::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 TStatToolkit::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]; | |
b8072cce | 190 | for (Int_t i=0;i<n;i++) sindexS[i]=0; |
191 | for (Int_t i=0;i<2*n;i++) sindexF[i]=0; | |
21f3a443 | 192 | // |
193 | TMath::Sort(n,inlist, sindexS, down); | |
194 | Int_t last = inlist[sindexS[0]]; | |
195 | Int_t val = last; | |
196 | sindexF[0] = 1; | |
197 | sindexF[0+n] = last; | |
198 | Int_t countPos = 0; | |
199 | // | |
200 | // find frequency | |
201 | for(Int_t i=1;i<n; i++){ | |
202 | val = inlist[sindexS[i]]; | |
203 | if (last == val) sindexF[countPos]++; | |
204 | else{ | |
205 | countPos++; | |
206 | sindexF[countPos+n] = val; | |
207 | sindexF[countPos]++; | |
208 | last =val; | |
209 | } | |
210 | } | |
211 | if (last==val) countPos++; | |
212 | // sort according frequency | |
213 | TMath::Sort(countPos, sindexF, sindexS, kTRUE); | |
214 | for (Int_t i=0;i<countPos;i++){ | |
215 | outlist[2*i ] = sindexF[sindexS[i]+n]; | |
216 | outlist[2*i+1] = sindexF[sindexS[i]]; | |
217 | } | |
218 | delete [] sindexS; | |
219 | delete [] sindexF; | |
220 | ||
221 | return countPos; | |
222 | ||
223 | } | |
224 | ||
225 | //___TStatToolkit__________________________________________________________________________ | |
3d7cc0b4 | 226 | void TStatToolkit::TruncatedMean(const TH1 * his, TVectorD *param, Float_t down, Float_t up, Bool_t verbose){ |
21f3a443 | 227 | // |
228 | // | |
229 | // | |
230 | Int_t nbins = his->GetNbinsX(); | |
231 | Float_t nentries = his->GetEntries(); | |
232 | Float_t sum =0; | |
233 | Float_t mean = 0; | |
234 | Float_t sigma2 = 0; | |
235 | Float_t ncumul=0; | |
236 | for (Int_t ibin=1;ibin<nbins; ibin++){ | |
237 | ncumul+= his->GetBinContent(ibin); | |
238 | Float_t fraction = Float_t(ncumul)/Float_t(nentries); | |
239 | if (fraction>down && fraction<up){ | |
240 | sum+=his->GetBinContent(ibin); | |
241 | mean+=his->GetBinCenter(ibin)*his->GetBinContent(ibin); | |
242 | sigma2+=his->GetBinCenter(ibin)*his->GetBinCenter(ibin)*his->GetBinContent(ibin); | |
243 | } | |
244 | } | |
245 | mean/=sum; | |
246 | sigma2= TMath::Sqrt(TMath::Abs(sigma2/sum-mean*mean)); | |
247 | if (param){ | |
248 | (*param)[0] = his->GetMaximum(); | |
249 | (*param)[1] = mean; | |
250 | (*param)[2] = sigma2; | |
251 | ||
252 | } | |
253 | if (verbose) printf("Mean\t%f\t Sigma2\t%f\n", mean,sigma2); | |
254 | } | |
255 | ||
256 | void TStatToolkit::LTM(TH1F * his, TVectorD *param , Float_t fraction, Bool_t verbose){ | |
257 | // | |
258 | // LTM | |
259 | // | |
260 | Int_t nbins = his->GetNbinsX(); | |
261 | Int_t nentries = (Int_t)his->GetEntries(); | |
262 | Double_t *data = new Double_t[nentries]; | |
263 | Int_t npoints=0; | |
264 | for (Int_t ibin=1;ibin<nbins; ibin++){ | |
265 | Float_t entriesI = his->GetBinContent(ibin); | |
266 | Float_t xcenter= his->GetBinCenter(ibin); | |
267 | for (Int_t ic=0; ic<entriesI; ic++){ | |
268 | if (npoints<nentries){ | |
269 | data[npoints]= xcenter; | |
270 | npoints++; | |
271 | } | |
272 | } | |
273 | } | |
274 | Double_t mean, sigma; | |
275 | Int_t npoints2=TMath::Min(Int_t(fraction*Float_t(npoints)),npoints-1); | |
276 | npoints2=TMath::Max(Int_t(0.5*Float_t(npoints)),npoints2); | |
277 | TStatToolkit::EvaluateUni(npoints, data, mean,sigma,npoints2); | |
278 | delete [] data; | |
279 | if (verbose) printf("Mean\t%f\t Sigma2\t%f\n", mean,sigma);if (param){ | |
280 | (*param)[0] = his->GetMaximum(); | |
281 | (*param)[1] = mean; | |
282 | (*param)[2] = sigma; | |
283 | } | |
284 | } | |
285 | ||
94a43b22 | 286 | Double_t TStatToolkit::FitGaus(TH1* his, TVectorD *param, TMatrixD */*matrix*/, Float_t xmin, Float_t xmax, Bool_t verbose){ |
21f3a443 | 287 | // |
288 | // Fit histogram with gaussian function | |
289 | // | |
290 | // Prameters: | |
291 | // return value- chi2 - if negative ( not enough points) | |
292 | // his - input histogram | |
293 | // param - vector with parameters | |
294 | // xmin, xmax - range to fit - if xmin=xmax=0 - the full histogram range used | |
295 | // Fitting: | |
296 | // 1. Step - make logarithm | |
297 | // 2. Linear fit (parabola) - more robust - always converge | |
298 | // 3. In case of small statistic bins are averaged | |
299 | // | |
300 | static TLinearFitter fitter(3,"pol2"); | |
301 | TVectorD par(3); | |
302 | TVectorD sigma(3); | |
303 | TMatrixD mat(3,3); | |
304 | if (his->GetMaximum()<4) return -1; | |
305 | if (his->GetEntries()<12) return -1; | |
306 | if (his->GetRMS()<mat.GetTol()) return -1; | |
307 | Float_t maxEstimate = his->GetEntries()*his->GetBinWidth(1)/TMath::Sqrt((TMath::TwoPi()*his->GetRMS())); | |
308 | Int_t dsmooth = TMath::Nint(6./TMath::Sqrt(maxEstimate)); | |
309 | ||
310 | if (maxEstimate<1) return -1; | |
311 | Int_t nbins = his->GetNbinsX(); | |
312 | Int_t npoints=0; | |
313 | // | |
314 | ||
315 | ||
316 | if (xmin>=xmax){ | |
317 | xmin = his->GetXaxis()->GetXmin(); | |
318 | xmax = his->GetXaxis()->GetXmax(); | |
319 | } | |
320 | for (Int_t iter=0; iter<2; iter++){ | |
321 | fitter.ClearPoints(); | |
322 | npoints=0; | |
323 | for (Int_t ibin=1;ibin<nbins+1; ibin++){ | |
324 | Int_t countB=1; | |
325 | Float_t entriesI = his->GetBinContent(ibin); | |
326 | for (Int_t delta = -dsmooth; delta<=dsmooth; delta++){ | |
327 | if (ibin+delta>1 &&ibin+delta<nbins-1){ | |
328 | entriesI += his->GetBinContent(ibin+delta); | |
329 | countB++; | |
330 | } | |
331 | } | |
332 | entriesI/=countB; | |
333 | Double_t xcenter= his->GetBinCenter(ibin); | |
334 | if (xcenter<xmin || xcenter>xmax) continue; | |
335 | Double_t error=1./TMath::Sqrt(countB); | |
336 | Float_t cont=2; | |
337 | if (iter>0){ | |
338 | if (par[0]+par[1]*xcenter+par[2]*xcenter*xcenter>20) return 0; | |
339 | cont = TMath::Exp(par[0]+par[1]*xcenter+par[2]*xcenter*xcenter); | |
340 | if (cont>1.) error = 1./TMath::Sqrt(cont*Float_t(countB)); | |
341 | } | |
342 | if (entriesI>1&&cont>1){ | |
343 | fitter.AddPoint(&xcenter,TMath::Log(Float_t(entriesI)),error); | |
344 | npoints++; | |
345 | } | |
346 | } | |
347 | if (npoints>3){ | |
348 | fitter.Eval(); | |
349 | fitter.GetParameters(par); | |
350 | }else{ | |
351 | break; | |
352 | } | |
353 | } | |
354 | if (npoints<=3){ | |
355 | return -1; | |
356 | } | |
357 | fitter.GetParameters(par); | |
358 | fitter.GetCovarianceMatrix(mat); | |
359 | if (TMath::Abs(par[1])<mat.GetTol()) return -1; | |
360 | if (TMath::Abs(par[2])<mat.GetTol()) return -1; | |
361 | Double_t chi2 = fitter.GetChisquare()/Float_t(npoints); | |
362 | //fitter.GetParameters(); | |
363 | if (!param) param = new TVectorD(3); | |
cb1d20de | 364 | // if (!matrix) matrix = new TMatrixD(3,3); // Covariance matrix to be implemented |
21f3a443 | 365 | (*param)[1] = par[1]/(-2.*par[2]); |
366 | (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2])); | |
367 | (*param)[0] = TMath::Exp(par[0]+ par[1]* (*param)[1] + par[2]*(*param)[1]*(*param)[1]); | |
368 | if (verbose){ | |
369 | par.Print(); | |
370 | mat.Print(); | |
371 | param->Print(); | |
372 | printf("Chi2=%f\n",chi2); | |
373 | TF1 * f1= new TF1("f1","[0]*exp(-(x-[1])^2/(2*[2]*[2]))",his->GetXaxis()->GetXmin(),his->GetXaxis()->GetXmax()); | |
374 | f1->SetParameter(0, (*param)[0]); | |
375 | f1->SetParameter(1, (*param)[1]); | |
376 | f1->SetParameter(2, (*param)[2]); | |
377 | f1->Draw("same"); | |
378 | } | |
379 | return chi2; | |
380 | } | |
381 | ||
cb1d20de | 382 | Double_t TStatToolkit::FitGaus(Float_t *arr, Int_t nBins, Float_t xMin, Float_t xMax, TVectorD *param, TMatrixD */*matrix*/, Bool_t verbose){ |
21f3a443 | 383 | // |
384 | // Fit histogram with gaussian function | |
385 | // | |
386 | // Prameters: | |
387 | // nbins: size of the array and number of histogram bins | |
388 | // xMin, xMax: histogram range | |
389 | // param: paramters of the fit (0-Constant, 1-Mean, 2-Sigma) | |
390 | // matrix: covariance matrix -- not implemented yet, pass dummy matrix!!! | |
391 | // | |
392 | // Return values: | |
393 | // >0: the chi2 returned by TLinearFitter | |
394 | // -3: only three points have been used for the calculation - no fitter was used | |
395 | // -2: only two points have been used for the calculation - center of gravity was uesed for calculation | |
396 | // -1: only one point has been used for the calculation - center of gravity was uesed for calculation | |
397 | // -4: invalid result!! | |
398 | // | |
399 | // Fitting: | |
400 | // 1. Step - make logarithm | |
401 | // 2. Linear fit (parabola) - more robust - always converge | |
402 | // | |
403 | static TLinearFitter fitter(3,"pol2"); | |
404 | static TMatrixD mat(3,3); | |
405 | static Double_t kTol = mat.GetTol(); | |
406 | fitter.StoreData(kFALSE); | |
407 | fitter.ClearPoints(); | |
408 | TVectorD par(3); | |
409 | TVectorD sigma(3); | |
3d7cc0b4 | 410 | TMatrixD matA(3,3); |
21f3a443 | 411 | TMatrixD b(3,1); |
412 | Float_t rms = TMath::RMS(nBins,arr); | |
413 | Float_t max = TMath::MaxElement(nBins,arr); | |
414 | Float_t binWidth = (xMax-xMin)/(Float_t)nBins; | |
415 | ||
416 | Float_t meanCOG = 0; | |
417 | Float_t rms2COG = 0; | |
418 | Float_t sumCOG = 0; | |
419 | ||
420 | Float_t entries = 0; | |
421 | Int_t nfilled=0; | |
422 | ||
423 | for (Int_t i=0; i<nBins; i++){ | |
424 | entries+=arr[i]; | |
425 | if (arr[i]>0) nfilled++; | |
426 | } | |
427 | ||
428 | if (max<4) return -4; | |
429 | if (entries<12) return -4; | |
430 | if (rms<kTol) return -4; | |
431 | ||
432 | Int_t npoints=0; | |
433 | // | |
434 | ||
435 | // | |
436 | for (Int_t ibin=0;ibin<nBins; ibin++){ | |
437 | Float_t entriesI = arr[ibin]; | |
438 | if (entriesI>1){ | |
439 | Double_t xcenter = xMin+(ibin+0.5)*binWidth; | |
440 | ||
441 | Float_t error = 1./TMath::Sqrt(entriesI); | |
442 | Float_t val = TMath::Log(Float_t(entriesI)); | |
443 | fitter.AddPoint(&xcenter,val,error); | |
444 | if (npoints<3){ | |
3d7cc0b4 | 445 | matA(npoints,0)=1; |
446 | matA(npoints,1)=xcenter; | |
447 | matA(npoints,2)=xcenter*xcenter; | |
21f3a443 | 448 | b(npoints,0)=val; |
449 | meanCOG+=xcenter*entriesI; | |
450 | rms2COG +=xcenter*entriesI*xcenter; | |
451 | sumCOG +=entriesI; | |
452 | } | |
453 | npoints++; | |
454 | } | |
455 | } | |
456 | ||
457 | ||
458 | Double_t chi2 = 0; | |
459 | if (npoints>=3){ | |
460 | if ( npoints == 3 ){ | |
461 | //analytic calculation of the parameters for three points | |
3d7cc0b4 | 462 | matA.Invert(); |
21f3a443 | 463 | TMatrixD res(1,3); |
3d7cc0b4 | 464 | res.Mult(matA,b); |
21f3a443 | 465 | par[0]=res(0,0); |
466 | par[1]=res(0,1); | |
467 | par[2]=res(0,2); | |
468 | chi2 = -3.; | |
469 | } else { | |
470 | // use fitter for more than three points | |
471 | fitter.Eval(); | |
472 | fitter.GetParameters(par); | |
473 | fitter.GetCovarianceMatrix(mat); | |
474 | chi2 = fitter.GetChisquare()/Float_t(npoints); | |
475 | } | |
476 | if (TMath::Abs(par[1])<kTol) return -4; | |
477 | if (TMath::Abs(par[2])<kTol) return -4; | |
478 | ||
479 | if (!param) param = new TVectorD(3); | |
cb1d20de | 480 | //if (!matrix) matrix = new TMatrixD(3,3); // !!!!might be a memory leek. use dummy matrix pointer to call this function! // Covariance matrix to be implemented |
21f3a443 | 481 | |
482 | (*param)[1] = par[1]/(-2.*par[2]); | |
483 | (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2])); | |
484 | Double_t lnparam0 = par[0]+ par[1]* (*param)[1] + par[2]*(*param)[1]*(*param)[1]; | |
485 | if ( lnparam0>307 ) return -4; | |
486 | (*param)[0] = TMath::Exp(lnparam0); | |
487 | if (verbose){ | |
488 | par.Print(); | |
489 | mat.Print(); | |
490 | param->Print(); | |
491 | printf("Chi2=%f\n",chi2); | |
492 | TF1 * f1= new TF1("f1","[0]*exp(-(x-[1])^2/(2*[2]*[2]))",xMin,xMax); | |
493 | f1->SetParameter(0, (*param)[0]); | |
494 | f1->SetParameter(1, (*param)[1]); | |
495 | f1->SetParameter(2, (*param)[2]); | |
496 | f1->Draw("same"); | |
497 | } | |
498 | return chi2; | |
499 | } | |
500 | ||
501 | if (npoints == 2){ | |
502 | //use center of gravity for 2 points | |
503 | meanCOG/=sumCOG; | |
504 | rms2COG /=sumCOG; | |
505 | (*param)[0] = max; | |
506 | (*param)[1] = meanCOG; | |
507 | (*param)[2] = TMath::Sqrt(TMath::Abs(meanCOG*meanCOG-rms2COG)); | |
508 | chi2=-2.; | |
509 | } | |
510 | if ( npoints == 1 ){ | |
511 | meanCOG/=sumCOG; | |
512 | (*param)[0] = max; | |
513 | (*param)[1] = meanCOG; | |
514 | (*param)[2] = binWidth/TMath::Sqrt(12); | |
515 | chi2=-1.; | |
516 | } | |
517 | return chi2; | |
518 | ||
519 | } | |
520 | ||
521 | ||
3d7cc0b4 | 522 | Float_t TStatToolkit::GetCOG(const Short_t *arr, Int_t nBins, Float_t xMin, Float_t xMax, Float_t *rms, Float_t *sum) |
21f3a443 | 523 | { |
524 | // | |
525 | // calculate center of gravity rms and sum for array 'arr' with nBins an a x range xMin to xMax | |
526 | // return COG; in case of failure return xMin | |
527 | // | |
528 | Float_t meanCOG = 0; | |
529 | Float_t rms2COG = 0; | |
530 | Float_t sumCOG = 0; | |
531 | Int_t npoints = 0; | |
532 | ||
533 | Float_t binWidth = (xMax-xMin)/(Float_t)nBins; | |
534 | ||
535 | for (Int_t ibin=0; ibin<nBins; ibin++){ | |
536 | Float_t entriesI = (Float_t)arr[ibin]; | |
537 | Double_t xcenter = xMin+(ibin+0.5)*binWidth; | |
538 | if ( entriesI>0 ){ | |
539 | meanCOG += xcenter*entriesI; | |
540 | rms2COG += xcenter*entriesI*xcenter; | |
541 | sumCOG += entriesI; | |
542 | npoints++; | |
543 | } | |
544 | } | |
545 | if ( sumCOG == 0 ) return xMin; | |
546 | meanCOG/=sumCOG; | |
547 | ||
548 | if ( rms ){ | |
549 | rms2COG /=sumCOG; | |
550 | (*rms) = TMath::Sqrt(TMath::Abs(meanCOG*meanCOG-rms2COG)); | |
551 | if ( npoints == 1 ) (*rms) = binWidth/TMath::Sqrt(12); | |
552 | } | |
553 | ||
554 | if ( sum ) | |
555 | (*sum) = sumCOG; | |
556 | ||
557 | return meanCOG; | |
558 | } | |
559 | ||
560 | ||
561 | ||
562 | /////////////////////////////////////////////////////////////// | |
563 | ////////////// TEST functions ///////////////////////// | |
564 | /////////////////////////////////////////////////////////////// | |
565 | ||
566 | ||
567 | ||
568 | ||
569 | ||
570 | void TStatToolkit::TestGausFit(Int_t nhistos){ | |
571 | // | |
572 | // Test performance of the parabolic - gaussian fit - compare it with | |
573 | // ROOT gauss fit | |
574 | // nhistos - number of histograms to be used for test | |
575 | // | |
576 | TTreeSRedirector *pcstream = new TTreeSRedirector("fitdebug.root"); | |
577 | ||
578 | Float_t *xTrue = new Float_t[nhistos]; | |
579 | Float_t *sTrue = new Float_t[nhistos]; | |
580 | TVectorD **par1 = new TVectorD*[nhistos]; | |
581 | TVectorD **par2 = new TVectorD*[nhistos]; | |
582 | TMatrixD dummy(3,3); | |
583 | ||
584 | ||
585 | TH1F **h1f = new TH1F*[nhistos]; | |
586 | TF1 *myg = new TF1("myg","gaus"); | |
587 | TF1 *fit = new TF1("fit","gaus"); | |
588 | gRandom->SetSeed(0); | |
589 | ||
590 | //init | |
591 | for (Int_t i=0;i<nhistos; i++){ | |
592 | par1[i] = new TVectorD(3); | |
593 | par2[i] = new TVectorD(3); | |
594 | h1f[i] = new TH1F(Form("h1f%d",i),Form("h1f%d",i),20,-10,10); | |
595 | xTrue[i]= gRandom->Rndm(); | |
596 | gSystem->Sleep(2); | |
597 | sTrue[i]= .75+gRandom->Rndm()*.5; | |
598 | myg->SetParameters(1,xTrue[i],sTrue[i]); | |
599 | h1f[i]->FillRandom("myg"); | |
600 | } | |
601 | ||
602 | TStopwatch s; | |
603 | s.Start(); | |
604 | //standard gaus fit | |
605 | for (Int_t i=0; i<nhistos; i++){ | |
606 | h1f[i]->Fit(fit,"0q"); | |
607 | (*par1[i])(0) = fit->GetParameter(0); | |
608 | (*par1[i])(1) = fit->GetParameter(1); | |
609 | (*par1[i])(2) = fit->GetParameter(2); | |
610 | } | |
611 | s.Stop(); | |
612 | printf("Gaussian fit\t"); | |
613 | s.Print(); | |
614 | ||
615 | s.Start(); | |
616 | //TStatToolkit gaus fit | |
617 | for (Int_t i=0; i<nhistos; i++){ | |
618 | TStatToolkit::FitGaus(h1f[i]->GetArray()+1,h1f[i]->GetNbinsX(),h1f[i]->GetXaxis()->GetXmin(),h1f[i]->GetXaxis()->GetXmax(),par2[i],&dummy); | |
619 | } | |
620 | ||
621 | s.Stop(); | |
622 | printf("Parabolic fit\t"); | |
623 | s.Print(); | |
624 | //write stream | |
625 | for (Int_t i=0;i<nhistos; i++){ | |
626 | Float_t xt = xTrue[i]; | |
627 | Float_t st = sTrue[i]; | |
628 | (*pcstream)<<"data" | |
629 | <<"xTrue="<<xt | |
630 | <<"sTrue="<<st | |
631 | <<"pg.="<<(par1[i]) | |
632 | <<"pa.="<<(par2[i]) | |
633 | <<"\n"; | |
634 | } | |
635 | //delete pointers | |
636 | for (Int_t i=0;i<nhistos; i++){ | |
637 | delete par1[i]; | |
638 | delete par2[i]; | |
639 | delete h1f[i]; | |
640 | } | |
641 | delete pcstream; | |
642 | delete []h1f; | |
643 | delete []xTrue; | |
644 | delete []sTrue; | |
645 | // | |
646 | delete []par1; | |
647 | delete []par2; | |
648 | ||
649 | } | |
650 | ||
651 | ||
652 | ||
653 | TGraph2D * TStatToolkit::MakeStat2D(TH3 * his, Int_t delta0, Int_t delta1, Int_t type){ | |
654 | // | |
655 | // | |
656 | // | |
657 | // delta - number of bins to integrate | |
658 | // type - 0 - mean value | |
659 | ||
660 | TAxis * xaxis = his->GetXaxis(); | |
661 | TAxis * yaxis = his->GetYaxis(); | |
662 | // TAxis * zaxis = his->GetZaxis(); | |
663 | Int_t nbinx = xaxis->GetNbins(); | |
664 | Int_t nbiny = yaxis->GetNbins(); | |
665 | char name[1000]; | |
666 | Int_t icount=0; | |
667 | TGraph2D *graph = new TGraph2D(nbinx*nbiny); | |
668 | TF1 f1("f1","gaus"); | |
669 | for (Int_t ix=0; ix<nbinx;ix++) | |
670 | for (Int_t iy=0; iy<nbiny;iy++){ | |
671 | Float_t xcenter = xaxis->GetBinCenter(ix); | |
672 | Float_t ycenter = yaxis->GetBinCenter(iy); | |
cb1d20de | 673 | snprintf(name,1000,"%s_%d_%d",his->GetName(), ix,iy); |
21f3a443 | 674 | TH1 *projection = his->ProjectionZ(name,ix-delta0,ix+delta0,iy-delta1,iy+delta1); |
675 | Float_t stat= 0; | |
676 | if (type==0) stat = projection->GetMean(); | |
677 | if (type==1) stat = projection->GetRMS(); | |
678 | if (type==2 || type==3){ | |
679 | TVectorD vec(3); | |
680 | TStatToolkit::LTM((TH1F*)projection,&vec,0.7); | |
681 | if (type==2) stat= vec[1]; | |
682 | if (type==3) stat= vec[0]; | |
683 | } | |
684 | if (type==4|| type==5){ | |
685 | projection->Fit(&f1); | |
686 | if (type==4) stat= f1.GetParameter(1); | |
687 | if (type==5) stat= f1.GetParameter(2); | |
688 | } | |
689 | //printf("%d\t%f\t%f\t%f\n", icount,xcenter, ycenter, stat); | |
690 | graph->SetPoint(icount,xcenter, ycenter, stat); | |
691 | icount++; | |
692 | } | |
693 | return graph; | |
694 | } | |
695 | ||
696 | TGraph * TStatToolkit::MakeStat1D(TH3 * his, Int_t delta1, Int_t type){ | |
697 | // | |
698 | // | |
699 | // | |
700 | // delta - number of bins to integrate | |
701 | // type - 0 - mean value | |
702 | ||
703 | TAxis * xaxis = his->GetXaxis(); | |
704 | TAxis * yaxis = his->GetYaxis(); | |
705 | // TAxis * zaxis = his->GetZaxis(); | |
706 | Int_t nbinx = xaxis->GetNbins(); | |
707 | Int_t nbiny = yaxis->GetNbins(); | |
708 | char name[1000]; | |
709 | Int_t icount=0; | |
710 | TGraph *graph = new TGraph(nbinx); | |
711 | TF1 f1("f1","gaus"); | |
712 | for (Int_t ix=0; ix<nbinx;ix++){ | |
713 | Float_t xcenter = xaxis->GetBinCenter(ix); | |
714 | // Float_t ycenter = yaxis->GetBinCenter(iy); | |
cb1d20de | 715 | snprintf(name,1000,"%s_%d",his->GetName(), ix); |
21f3a443 | 716 | TH1 *projection = his->ProjectionZ(name,ix-delta1,ix+delta1,0,nbiny); |
717 | Float_t stat= 0; | |
718 | if (type==0) stat = projection->GetMean(); | |
719 | if (type==1) stat = projection->GetRMS(); | |
720 | if (type==2 || type==3){ | |
721 | TVectorD vec(3); | |
722 | TStatToolkit::LTM((TH1F*)projection,&vec,0.7); | |
723 | if (type==2) stat= vec[1]; | |
724 | if (type==3) stat= vec[0]; | |
725 | } | |
726 | if (type==4|| type==5){ | |
727 | projection->Fit(&f1); | |
728 | if (type==4) stat= f1.GetParameter(1); | |
729 | if (type==5) stat= f1.GetParameter(2); | |
730 | } | |
731 | //printf("%d\t%f\t%f\t%f\n", icount,xcenter, ycenter, stat); | |
732 | graph->SetPoint(icount,xcenter, stat); | |
733 | icount++; | |
734 | } | |
735 | return graph; | |
736 | } | |
737 | ||
738 | ||
739 | ||
740 | ||
741 | ||
88b1c775 | 742 | TString* TStatToolkit::FitPlane(TTree *tree, const char* drawCommand, const char* formula, const char* cuts, Double_t & chi2, Int_t &npoints, TVectorD &fitParam, TMatrixD &covMatrix, Float_t frac, Int_t start, Int_t stop,Bool_t fix0){ |
21f3a443 | 743 | // |
744 | // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts | |
745 | // returns chi2, fitParam and covMatrix | |
746 | // returns TString with fitted formula | |
747 | // | |
dd46129c | 748 | |
21f3a443 | 749 | TString formulaStr(formula); |
750 | TString drawStr(drawCommand); | |
751 | TString cutStr(cuts); | |
dd46129c | 752 | TString ferr("1"); |
753 | ||
754 | TString strVal(drawCommand); | |
755 | if (strVal.Contains(":")){ | |
756 | TObjArray* valTokens = strVal.Tokenize(":"); | |
757 | drawStr = valTokens->At(0)->GetName(); | |
758 | ferr = valTokens->At(1)->GetName(); | |
759 | } | |
760 | ||
21f3a443 | 761 | |
762 | formulaStr.ReplaceAll("++", "~"); | |
763 | TObjArray* formulaTokens = formulaStr.Tokenize("~"); | |
764 | Int_t dim = formulaTokens->GetEntriesFast(); | |
765 | ||
766 | fitParam.ResizeTo(dim); | |
767 | covMatrix.ResizeTo(dim,dim); | |
768 | ||
769 | TLinearFitter* fitter = new TLinearFitter(dim+1, Form("hyp%d",dim)); | |
770 | fitter->StoreData(kTRUE); | |
771 | fitter->ClearPoints(); | |
772 | ||
773 | Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start, start); | |
774 | if (entries == -1) return new TString("An ERROR has occured during fitting!"); | |
bd7b4d18 | 775 | Double_t **values = new Double_t*[dim+1] ; |
776 | for (Int_t i=0; i<dim+1; i++) values[i]=NULL; | |
dd46129c | 777 | // |
778 | entries = tree->Draw(ferr.Data(), cutStr.Data(), "goff", stop-start, start); | |
b8072cce | 779 | if (entries == -1) { |
780 | delete []values; | |
781 | return new TString("An ERROR has occured during fitting!"); | |
782 | } | |
dd46129c | 783 | Double_t *errors = new Double_t[entries]; |
784 | memcpy(errors, tree->GetV1(), entries*sizeof(Double_t)); | |
21f3a443 | 785 | |
786 | for (Int_t i = 0; i < dim + 1; i++){ | |
787 | Int_t centries = 0; | |
788 | if (i < dim) centries = tree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start); | |
789 | else centries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start); | |
790 | ||
b8072cce | 791 | if (entries != centries) { |
792 | delete []errors; | |
793 | delete []values; | |
794 | return new TString("An ERROR has occured during fitting!"); | |
795 | } | |
21f3a443 | 796 | values[i] = new Double_t[entries]; |
797 | memcpy(values[i], tree->GetV1(), entries*sizeof(Double_t)); | |
798 | } | |
799 | ||
800 | // add points to the fitter | |
801 | for (Int_t i = 0; i < entries; i++){ | |
802 | Double_t x[1000]; | |
803 | for (Int_t j=0; j<dim;j++) x[j]=values[j][i]; | |
dd46129c | 804 | fitter->AddPoint(x, values[dim][i], errors[i]); |
21f3a443 | 805 | } |
806 | ||
807 | fitter->Eval(); | |
2c629c56 | 808 | if (frac>0.5 && frac<1){ |
809 | fitter->EvalRobust(frac); | |
88b1c775 | 810 | }else{ |
811 | if (fix0) { | |
812 | fitter->FixParameter(0,0); | |
813 | fitter->Eval(); | |
814 | } | |
2c629c56 | 815 | } |
21f3a443 | 816 | fitter->GetParameters(fitParam); |
817 | fitter->GetCovarianceMatrix(covMatrix); | |
818 | chi2 = fitter->GetChisquare(); | |
b8072cce | 819 | npoints = entries; |
21f3a443 | 820 | TString *preturnFormula = new TString(Form("( %f+",fitParam[0])), &returnFormula = *preturnFormula; |
821 | ||
822 | for (Int_t iparam = 0; iparam < dim; iparam++) { | |
823 | returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam+1])); | |
824 | if (iparam < dim-1) returnFormula.Append("+"); | |
825 | } | |
826 | returnFormula.Append(" )"); | |
4d61c301 | 827 | |
828 | ||
b8072cce | 829 | for (Int_t j=0; j<dim+1;j++) delete [] values[j]; |
4d61c301 | 830 | |
831 | ||
cb1d20de | 832 | delete formulaTokens; |
833 | delete fitter; | |
834 | delete[] values; | |
b8072cce | 835 | delete[] errors; |
cb1d20de | 836 | return preturnFormula; |
837 | } | |
838 | ||
839 | TString* TStatToolkit::FitPlaneConstrain(TTree *tree, const char* drawCommand, const char* formula, const char* cuts, Double_t & chi2, Int_t &npoints, TVectorD &fitParam, TMatrixD &covMatrix, Float_t frac, Int_t start, Int_t stop,Double_t constrain){ | |
840 | // | |
841 | // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts | |
842 | // returns chi2, fitParam and covMatrix | |
843 | // returns TString with fitted formula | |
844 | // | |
845 | ||
846 | TString formulaStr(formula); | |
847 | TString drawStr(drawCommand); | |
848 | TString cutStr(cuts); | |
849 | TString ferr("1"); | |
850 | ||
851 | TString strVal(drawCommand); | |
852 | if (strVal.Contains(":")){ | |
853 | TObjArray* valTokens = strVal.Tokenize(":"); | |
854 | drawStr = valTokens->At(0)->GetName(); | |
855 | ferr = valTokens->At(1)->GetName(); | |
856 | } | |
857 | ||
858 | ||
859 | formulaStr.ReplaceAll("++", "~"); | |
860 | TObjArray* formulaTokens = formulaStr.Tokenize("~"); | |
861 | Int_t dim = formulaTokens->GetEntriesFast(); | |
862 | ||
863 | fitParam.ResizeTo(dim); | |
864 | covMatrix.ResizeTo(dim,dim); | |
865 | ||
866 | TLinearFitter* fitter = new TLinearFitter(dim+1, Form("hyp%d",dim)); | |
867 | fitter->StoreData(kTRUE); | |
868 | fitter->ClearPoints(); | |
869 | ||
870 | Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start, start); | |
871 | if (entries == -1) return new TString("An ERROR has occured during fitting!"); | |
872 | Double_t **values = new Double_t*[dim+1] ; | |
bd7b4d18 | 873 | for (Int_t i=0; i<dim+1; i++) values[i]=NULL; |
cb1d20de | 874 | // |
875 | entries = tree->Draw(ferr.Data(), cutStr.Data(), "goff", stop-start, start); | |
b8072cce | 876 | if (entries == -1) { |
877 | delete [] values; | |
878 | return new TString("An ERROR has occured during fitting!"); | |
879 | } | |
cb1d20de | 880 | Double_t *errors = new Double_t[entries]; |
881 | memcpy(errors, tree->GetV1(), entries*sizeof(Double_t)); | |
882 | ||
883 | for (Int_t i = 0; i < dim + 1; i++){ | |
884 | Int_t centries = 0; | |
885 | if (i < dim) centries = tree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start); | |
886 | else centries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start); | |
887 | ||
b8072cce | 888 | if (entries != centries) { |
889 | delete []errors; | |
890 | delete []values; | |
891 | return new TString("An ERROR has occured during fitting!"); | |
892 | } | |
cb1d20de | 893 | values[i] = new Double_t[entries]; |
894 | memcpy(values[i], tree->GetV1(), entries*sizeof(Double_t)); | |
895 | } | |
896 | ||
897 | // add points to the fitter | |
898 | for (Int_t i = 0; i < entries; i++){ | |
899 | Double_t x[1000]; | |
900 | for (Int_t j=0; j<dim;j++) x[j]=values[j][i]; | |
901 | fitter->AddPoint(x, values[dim][i], errors[i]); | |
902 | } | |
903 | if (constrain>0){ | |
904 | for (Int_t i = 0; i < dim; i++){ | |
905 | Double_t x[1000]; | |
906 | for (Int_t j=0; j<dim;j++) if (i!=j) x[j]=0; | |
907 | x[i]=1.; | |
908 | fitter->AddPoint(x, 0, constrain); | |
909 | } | |
910 | } | |
911 | ||
912 | ||
913 | fitter->Eval(); | |
914 | if (frac>0.5 && frac<1){ | |
915 | fitter->EvalRobust(frac); | |
916 | } | |
917 | fitter->GetParameters(fitParam); | |
918 | fitter->GetCovarianceMatrix(covMatrix); | |
919 | chi2 = fitter->GetChisquare(); | |
920 | npoints = entries; | |
cb1d20de | 921 | |
922 | TString *preturnFormula = new TString(Form("( %f+",fitParam[0])), &returnFormula = *preturnFormula; | |
923 | ||
924 | for (Int_t iparam = 0; iparam < dim; iparam++) { | |
925 | returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam+1])); | |
926 | if (iparam < dim-1) returnFormula.Append("+"); | |
927 | } | |
928 | returnFormula.Append(" )"); | |
929 | ||
b8072cce | 930 | for (Int_t j=0; j<dim+1;j++) delete [] values[j]; |
cb1d20de | 931 | |
932 | ||
933 | ||
934 | delete formulaTokens; | |
935 | delete fitter; | |
936 | delete[] values; | |
b8072cce | 937 | delete[] errors; |
cb1d20de | 938 | return preturnFormula; |
939 | } | |
940 | ||
941 | ||
942 | ||
943 | TString* TStatToolkit::FitPlaneFixed(TTree *tree, const char* drawCommand, const char* formula, const char* cuts, Double_t & chi2, Int_t &npoints, TVectorD &fitParam, TMatrixD &covMatrix, Float_t frac, Int_t start, Int_t stop){ | |
944 | // | |
945 | // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts | |
946 | // returns chi2, fitParam and covMatrix | |
947 | // returns TString with fitted formula | |
948 | // | |
949 | ||
950 | TString formulaStr(formula); | |
951 | TString drawStr(drawCommand); | |
952 | TString cutStr(cuts); | |
953 | TString ferr("1"); | |
954 | ||
955 | TString strVal(drawCommand); | |
956 | if (strVal.Contains(":")){ | |
957 | TObjArray* valTokens = strVal.Tokenize(":"); | |
958 | drawStr = valTokens->At(0)->GetName(); | |
959 | ferr = valTokens->At(1)->GetName(); | |
960 | } | |
961 | ||
962 | ||
963 | formulaStr.ReplaceAll("++", "~"); | |
964 | TObjArray* formulaTokens = formulaStr.Tokenize("~"); | |
965 | Int_t dim = formulaTokens->GetEntriesFast(); | |
966 | ||
967 | fitParam.ResizeTo(dim); | |
968 | covMatrix.ResizeTo(dim,dim); | |
969 | TString fitString="x0"; | |
970 | for (Int_t i=1; i<dim; i++) fitString+=Form("++x%d",i); | |
971 | TLinearFitter* fitter = new TLinearFitter(dim, fitString.Data()); | |
972 | fitter->StoreData(kTRUE); | |
973 | fitter->ClearPoints(); | |
974 | ||
975 | Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start, start); | |
976 | if (entries == -1) return new TString("An ERROR has occured during fitting!"); | |
977 | Double_t **values = new Double_t*[dim+1] ; | |
bd7b4d18 | 978 | for (Int_t i=0; i<dim+1; i++) values[i]=NULL; |
cb1d20de | 979 | // |
980 | entries = tree->Draw(ferr.Data(), cutStr.Data(), "goff", stop-start, start); | |
b8072cce | 981 | if (entries == -1) { |
982 | delete []values; | |
983 | return new TString("An ERROR has occured during fitting!"); | |
984 | } | |
cb1d20de | 985 | Double_t *errors = new Double_t[entries]; |
986 | memcpy(errors, tree->GetV1(), entries*sizeof(Double_t)); | |
987 | ||
988 | for (Int_t i = 0; i < dim + 1; i++){ | |
989 | Int_t centries = 0; | |
990 | if (i < dim) centries = tree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start); | |
991 | else centries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start); | |
992 | ||
b8072cce | 993 | if (entries != centries) { |
994 | delete []errors; | |
995 | delete []values; | |
996 | return new TString("An ERROR has occured during fitting!"); | |
997 | } | |
cb1d20de | 998 | values[i] = new Double_t[entries]; |
999 | memcpy(values[i], tree->GetV1(), entries*sizeof(Double_t)); | |
1000 | } | |
1001 | ||
1002 | // add points to the fitter | |
1003 | for (Int_t i = 0; i < entries; i++){ | |
1004 | Double_t x[1000]; | |
1005 | for (Int_t j=0; j<dim;j++) x[j]=values[j][i]; | |
1006 | fitter->AddPoint(x, values[dim][i], errors[i]); | |
1007 | } | |
1008 | ||
1009 | fitter->Eval(); | |
1010 | if (frac>0.5 && frac<1){ | |
1011 | fitter->EvalRobust(frac); | |
1012 | } | |
1013 | fitter->GetParameters(fitParam); | |
1014 | fitter->GetCovarianceMatrix(covMatrix); | |
1015 | chi2 = fitter->GetChisquare(); | |
1016 | npoints = entries; | |
cb1d20de | 1017 | |
1018 | TString *preturnFormula = new TString("("), &returnFormula = *preturnFormula; | |
1019 | ||
1020 | for (Int_t iparam = 0; iparam < dim; iparam++) { | |
1021 | returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam])); | |
1022 | if (iparam < dim-1) returnFormula.Append("+"); | |
1023 | } | |
1024 | returnFormula.Append(" )"); | |
1025 | ||
1026 | ||
b8072cce | 1027 | for (Int_t j=0; j<dim+1;j++) delete [] values[j]; |
cb1d20de | 1028 | |
21f3a443 | 1029 | delete formulaTokens; |
1030 | delete fitter; | |
1031 | delete[] values; | |
b8072cce | 1032 | delete[] errors; |
21f3a443 | 1033 | return preturnFormula; |
1034 | } | |
7c9cf6e4 | 1035 | |
1036 | ||
1037 | ||
1038 | ||
1039 | ||
3d7cc0b4 | 1040 | Int_t TStatToolkit::GetFitIndex(const TString fString, const TString subString){ |
7c9cf6e4 | 1041 | // |
1042 | // fitString - ++ separated list of fits | |
1043 | // substring - ++ separated list of the requiered substrings | |
1044 | // | |
1045 | // return the last occurance of substring in fit string | |
1046 | // | |
1047 | TObjArray *arrFit = fString.Tokenize("++"); | |
1048 | TObjArray *arrSub = subString.Tokenize("++"); | |
1049 | Int_t index=-1; | |
1050 | for (Int_t i=0; i<arrFit->GetEntries(); i++){ | |
1051 | Bool_t isOK=kTRUE; | |
1052 | TString str =arrFit->At(i)->GetName(); | |
1053 | for (Int_t isub=0; isub<arrSub->GetEntries(); isub++){ | |
1054 | if (str.Contains(arrSub->At(isub)->GetName())==0) isOK=kFALSE; | |
1055 | } | |
1056 | if (isOK) index=i; | |
1057 | } | |
1058 | return index; | |
1059 | } | |
1060 | ||
1061 | ||
3d7cc0b4 | 1062 | TString TStatToolkit::FilterFit(const TString &input, const TString filter, TVectorD ¶m, TMatrixD & covar){ |
7c9cf6e4 | 1063 | // |
1064 | // Filter fit expression make sub-fit | |
1065 | // | |
1066 | TObjArray *array0= input.Tokenize("++"); | |
1067 | TObjArray *array1= filter.Tokenize("++"); | |
1068 | //TString *presult=new TString("(0"); | |
1069 | TString result="(0.0"; | |
1070 | for (Int_t i=0; i<array0->GetEntries(); i++){ | |
1071 | Bool_t isOK=kTRUE; | |
1072 | TString str(array0->At(i)->GetName()); | |
1073 | for (Int_t j=0; j<array1->GetEntries(); j++){ | |
1074 | if (str.Contains(array1->At(j)->GetName())==0) isOK=kFALSE; | |
1075 | } | |
1076 | if (isOK) { | |
1077 | result+="+"+str; | |
1078 | result+=Form("*(%f)",param[i+1]); | |
1079 | printf("%f\t%f\t%s\n",param[i+1], TMath::Sqrt(covar(i+1,i+1)),str.Data()); | |
1080 | } | |
1081 | } | |
1082 | result+="-0.)"; | |
1083 | return result; | |
1084 | } | |
1085 | ||
1086 | void TStatToolkit::Update1D(Double_t delta, Double_t sigma, Int_t s1, TMatrixD &vecXk, TMatrixD &covXk){ | |
1087 | // | |
1088 | // Update parameters and covariance - with one measurement | |
1089 | // Input: | |
1090 | // vecXk - input vector - Updated in function | |
1091 | // covXk - covariance matrix - Updated in function | |
1092 | // delta, sigma, s1 - new measurement, rms of new measurement and the index of measurement | |
1093 | const Int_t knMeas=1; | |
1094 | Int_t knElem=vecXk.GetNrows(); | |
1095 | ||
1096 | TMatrixD mat1(knElem,knElem); // update covariance matrix | |
1097 | TMatrixD matHk(1,knElem); // vector to mesurement | |
1098 | TMatrixD vecYk(knMeas,1); // Innovation or measurement residual | |
1099 | TMatrixD matHkT(knElem,knMeas); // helper matrix Hk transpose | |
1100 | TMatrixD matSk(knMeas,knMeas); // Innovation (or residual) covariance | |
1101 | TMatrixD matKk(knElem,knMeas); // Optimal Kalman gain | |
1102 | TMatrixD covXk2(knElem,knElem); // helper matrix | |
1103 | TMatrixD covXk3(knElem,knElem); // helper matrix | |
1104 | TMatrixD vecZk(1,1); | |
1105 | TMatrixD measR(1,1); | |
1106 | vecZk(0,0)=delta; | |
1107 | measR(0,0)=sigma*sigma; | |
1108 | // | |
1109 | // reset matHk | |
1110 | for (Int_t iel=0;iel<knElem;iel++) | |
1111 | for (Int_t ip=0;ip<knMeas;ip++) matHk(ip,iel)=0; | |
1112 | //mat1 | |
1113 | for (Int_t iel=0;iel<knElem;iel++) { | |
1114 | for (Int_t jel=0;jel<knElem;jel++) mat1(iel,jel)=0; | |
1115 | mat1(iel,iel)=1; | |
1116 | } | |
1117 | // | |
1118 | matHk(0, s1)=1; | |
1119 | vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual | |
1120 | matHkT=matHk.T(); matHk.T(); | |
1121 | matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance | |
1122 | matSk.Invert(); | |
1123 | matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain | |
1124 | vecXk += matKk*vecYk; // updated vector | |
1125 | covXk2= (mat1-(matKk*matHk)); | |
1126 | covXk3 = covXk2*covXk; | |
1127 | covXk = covXk3; | |
1128 | Int_t nrows=covXk3.GetNrows(); | |
1129 | ||
1130 | for (Int_t irow=0; irow<nrows; irow++) | |
1131 | for (Int_t icol=0; icol<nrows; icol++){ | |
1132 | // rounding problems - make matrix again symteric | |
1133 | covXk(irow,icol)=(covXk3(irow,icol)+covXk3(icol,irow))*0.5; | |
1134 | } | |
1135 | } | |
1136 | ||
1137 | ||
1138 | ||
3d7cc0b4 | 1139 | void TStatToolkit::Constrain1D(const TString &input, const TString filter, TVectorD ¶m, TMatrixD & covar, Double_t mean, Double_t sigma){ |
7c9cf6e4 | 1140 | // |
1141 | // constrain linear fit | |
1142 | // input - string description of fit function | |
1143 | // filter - string filter to select sub fits | |
1144 | // param,covar - parameters and covariance matrix of the fit | |
1145 | // mean,sigma - new measurement uning which the fit is updated | |
1146 | // | |
ae45c94d | 1147 | |
7c9cf6e4 | 1148 | TObjArray *array0= input.Tokenize("++"); |
1149 | TObjArray *array1= filter.Tokenize("++"); | |
1150 | TMatrixD paramM(param.GetNrows(),1); | |
1151 | for (Int_t i=0; i<=array0->GetEntries(); i++){paramM(i,0)=param(i);} | |
1152 | ||
ae45c94d | 1153 | if (filter.Length()==0){ |
1154 | TStatToolkit::Update1D(mean, sigma, 0, paramM, covar);// | |
1155 | }else{ | |
1156 | for (Int_t i=0; i<array0->GetEntries(); i++){ | |
1157 | Bool_t isOK=kTRUE; | |
1158 | TString str(array0->At(i)->GetName()); | |
1159 | for (Int_t j=0; j<array1->GetEntries(); j++){ | |
1160 | if (str.Contains(array1->At(j)->GetName())==0) isOK=kFALSE; | |
1161 | } | |
1162 | if (isOK) { | |
1163 | TStatToolkit::Update1D(mean, sigma, i+1, paramM, covar);// | |
1164 | } | |
7c9cf6e4 | 1165 | } |
1166 | } | |
1167 | for (Int_t i=0; i<=array0->GetEntries(); i++){ | |
1168 | param(i)=paramM(i,0); | |
1169 | } | |
1170 | } | |
1171 | ||
ae45c94d | 1172 | TString TStatToolkit::MakeFitString(const TString &input, const TVectorD ¶m, const TMatrixD & covar, Bool_t verbose){ |
7c9cf6e4 | 1173 | // |
1174 | // | |
1175 | // | |
1176 | TObjArray *array0= input.Tokenize("++"); | |
ae45c94d | 1177 | TString result=Form("(%f",param[0]); |
1178 | printf("%f\t%f\t\n", param[0], TMath::Sqrt(covar(0,0))); | |
7c9cf6e4 | 1179 | for (Int_t i=0; i<array0->GetEntries(); i++){ |
1180 | TString str(array0->At(i)->GetName()); | |
1181 | result+="+"+str; | |
1182 | result+=Form("*(%f)",param[i+1]); | |
ae45c94d | 1183 | if (verbose) printf("%f\t%f\t%s\n", param[i+1], TMath::Sqrt(covar(i+1,i+1)),str.Data()); |
7c9cf6e4 | 1184 | } |
1185 | result+="-0.)"; | |
1186 | return result; | |
1187 | } | |
df0a2a0a | 1188 | |
1189 | ||
b3453fe7 | 1190 | TGraph * TStatToolkit::MakeGraphSparse(TTree * tree, const char * expr, const char * cut, Int_t mstyle, Int_t mcolor){ |
df0a2a0a | 1191 | // |
1192 | // Make a sparse draw of the variables | |
b3453fe7 | 1193 | // Writen by Weilin.Yu |
df0a2a0a | 1194 | const Int_t entries = tree->Draw(expr,cut,"goff"); |
1195 | // TGraph * graph = (TGraph*)gPad->GetPrimitive("Graph"); // 2D | |
b3453fe7 | 1196 | TGraph * graph = 0; |
1197 | if (tree->GetV3()) graph = new TGraphErrors (entries, tree->GetV2(),tree->GetV1(),0,tree->GetV3()); | |
1198 | graph = new TGraphErrors (entries, tree->GetV2(),tree->GetV1(),0,0); | |
1199 | graph->SetMarkerStyle(mstyle); | |
1200 | graph->SetMarkerColor(mcolor); | |
df0a2a0a | 1201 | // |
1202 | Int_t *index = new Int_t[entries]; | |
1203 | TMath::Sort(entries,graph->GetX(),index,kFALSE); | |
1204 | ||
3d7cc0b4 | 1205 | Double_t *tempArray = new Double_t[entries]; |
df0a2a0a | 1206 | |
1207 | Double_t count = 0.5; | |
baa0041d | 1208 | Double_t *vrun = new Double_t[entries]; |
1209 | Int_t icount=0; | |
1210 | // | |
3d7cc0b4 | 1211 | tempArray[index[0]] = count; |
baa0041d | 1212 | vrun[0] = graph->GetX()[index[0]]; |
df0a2a0a | 1213 | for(Int_t i=1;i<entries;i++){ |
1214 | if(graph->GetX()[index[i]]==graph->GetX()[index[i-1]]) | |
3d7cc0b4 | 1215 | tempArray[index[i]] = count; |
df0a2a0a | 1216 | else if(graph->GetX()[index[i]]!=graph->GetX()[index[i-1]]){ |
1217 | count++; | |
baa0041d | 1218 | icount++; |
3d7cc0b4 | 1219 | tempArray[index[i]] = count; |
baa0041d | 1220 | vrun[icount]=graph->GetX()[index[i]]; |
df0a2a0a | 1221 | } |
1222 | } | |
1223 | ||
1224 | const Int_t newNbins = int(count+0.5); | |
1225 | Double_t *newBins = new Double_t[newNbins+1]; | |
1226 | for(Int_t i=0; i<=count+1;i++){ | |
1227 | newBins[i] = i; | |
1228 | } | |
1229 | ||
b3453fe7 | 1230 | TGraph *graphNew = 0; |
1231 | if (tree->GetV3()) graphNew = new TGraphErrors(entries,tempArray,graph->GetY(),0,tree->GetV3()); | |
1232 | else | |
1233 | graphNew = new TGraphErrors(entries,tempArray,graph->GetY(),0,0); | |
df0a2a0a | 1234 | graphNew->GetXaxis()->Set(newNbins,newBins); |
1235 | ||
1236 | Char_t xName[50]; | |
df0a2a0a | 1237 | for(Int_t i=0;i<count;i++){ |
baa0041d | 1238 | snprintf(xName,50,"%d",Int_t(vrun[i])); |
df0a2a0a | 1239 | graphNew->GetXaxis()->SetBinLabel(i+1,xName); |
1240 | } | |
1241 | graphNew->GetHistogram()->SetTitle(""); | |
b3453fe7 | 1242 | graphNew->SetMarkerStyle(mstyle); |
1243 | graphNew->SetMarkerColor(mcolor); | |
3d7cc0b4 | 1244 | delete [] tempArray; |
df0a2a0a | 1245 | delete [] index; |
1246 | delete [] newBins; | |
baa0041d | 1247 | delete [] vrun; |
df0a2a0a | 1248 | return graphNew; |
1249 | } | |
1250 |