]>
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" | |
3240a856 | 27 | #include "TH2F.h" |
21f3a443 | 28 | #include "TH3.h" |
29 | #include "TF1.h" | |
30 | #include "TTree.h" | |
31 | #include "TChain.h" | |
32 | #include "TObjString.h" | |
33 | #include "TLinearFitter.h" | |
3d7cc0b4 | 34 | #include "TGraph2D.h" |
35 | #include "TGraph.h" | |
b3453fe7 | 36 | #include "TGraphErrors.h" |
377a7d60 | 37 | #include "TMultiGraph.h" |
38 | #include "TCanvas.h" | |
39 | #include "TLatex.h" | |
40 | #include "TCut.h" | |
21f3a443 | 41 | |
42 | // | |
43 | // includes neccessary for test functions | |
44 | // | |
45 | #include "TSystem.h" | |
46 | #include "TRandom.h" | |
47 | #include "TStopwatch.h" | |
48 | #include "TTreeStream.h" | |
49 | ||
50 | #include "TStatToolkit.h" | |
51 | ||
52 | ||
53 | ClassImp(TStatToolkit) // Class implementation to enable ROOT I/O | |
54 | ||
55 | TStatToolkit::TStatToolkit() : TObject() | |
56 | { | |
57 | // | |
58 | // Default constructor | |
59 | // | |
60 | } | |
61 | /////////////////////////////////////////////////////////////////////////// | |
62 | TStatToolkit::~TStatToolkit() | |
63 | { | |
64 | // | |
65 | // Destructor | |
66 | // | |
67 | } | |
68 | ||
69 | ||
70 | //_____________________________________________________________________________ | |
71 | void TStatToolkit::EvaluateUni(Int_t nvectors, Double_t *data, Double_t &mean | |
72 | , Double_t &sigma, Int_t hh) | |
73 | { | |
74 | // | |
75 | // Robust estimator in 1D case MI version - (faster than ROOT version) | |
76 | // | |
77 | // For the univariate case | |
78 | // estimates of location and scatter are returned in mean and sigma parameters | |
79 | // the algorithm works on the same principle as in multivariate case - | |
80 | // it finds a subset of size hh with smallest sigma, and then returns mean and | |
81 | // sigma of this subset | |
82 | // | |
83 | ||
84 | if (hh==0) | |
85 | hh=(nvectors+2)/2; | |
86 | 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}; | |
87 | Int_t *index=new Int_t[nvectors]; | |
88 | TMath::Sort(nvectors, data, index, kFALSE); | |
89 | ||
90 | Int_t nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11); | |
91 | Double_t factor = faclts[TMath::Max(0,nquant-1)]; | |
92 | ||
93 | Double_t sumx =0; | |
94 | Double_t sumx2 =0; | |
95 | Int_t bestindex = -1; | |
96 | Double_t bestmean = 0; | |
97 | Double_t bestsigma = (data[index[nvectors-1]]-data[index[0]]+1.); // maximal possible sigma | |
98 | bestsigma *=bestsigma; | |
99 | ||
100 | for (Int_t i=0; i<hh; i++){ | |
101 | sumx += data[index[i]]; | |
102 | sumx2 += data[index[i]]*data[index[i]]; | |
103 | } | |
104 | ||
105 | Double_t norm = 1./Double_t(hh); | |
bd7b4d18 | 106 | Double_t norm2 = (hh-1)>0 ? 1./Double_t(hh-1):1; |
21f3a443 | 107 | for (Int_t i=hh; i<nvectors; i++){ |
108 | Double_t cmean = sumx*norm; | |
109 | Double_t csigma = (sumx2 - hh*cmean*cmean)*norm2; | |
110 | if (csigma<bestsigma){ | |
111 | bestmean = cmean; | |
112 | bestsigma = csigma; | |
113 | bestindex = i-hh; | |
114 | } | |
115 | ||
116 | sumx += data[index[i]]-data[index[i-hh]]; | |
117 | sumx2 += data[index[i]]*data[index[i]]-data[index[i-hh]]*data[index[i-hh]]; | |
118 | } | |
119 | ||
120 | Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma)); | |
121 | mean = bestmean; | |
122 | sigma = bstd; | |
123 | delete [] index; | |
124 | ||
125 | } | |
126 | ||
127 | ||
128 | ||
129 | void TStatToolkit::EvaluateUniExternal(Int_t nvectors, Double_t *data, Double_t &mean, Double_t &sigma, Int_t hh, Float_t externalfactor) | |
130 | { | |
131 | // Modified version of ROOT robust EvaluateUni | |
132 | // robust estimator in 1D case MI version | |
133 | // added external factor to include precision of external measurement | |
134 | // | |
135 | ||
136 | if (hh==0) | |
137 | hh=(nvectors+2)/2; | |
138 | 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}; | |
139 | Int_t *index=new Int_t[nvectors]; | |
140 | TMath::Sort(nvectors, data, index, kFALSE); | |
141 | // | |
142 | Int_t nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11); | |
143 | Double_t factor = faclts[0]; | |
144 | if (nquant>0){ | |
145 | // fix proper normalization - Anja | |
146 | factor = faclts[nquant-1]; | |
147 | } | |
148 | ||
149 | // | |
150 | // | |
151 | Double_t sumx =0; | |
152 | Double_t sumx2 =0; | |
153 | Int_t bestindex = -1; | |
154 | Double_t bestmean = 0; | |
155 | Double_t bestsigma = -1; | |
156 | for (Int_t i=0; i<hh; i++){ | |
157 | sumx += data[index[i]]; | |
158 | sumx2 += data[index[i]]*data[index[i]]; | |
159 | } | |
160 | // | |
161 | Double_t kfactor = 2.*externalfactor - externalfactor*externalfactor; | |
162 | Double_t norm = 1./Double_t(hh); | |
163 | for (Int_t i=hh; i<nvectors; i++){ | |
164 | Double_t cmean = sumx*norm; | |
165 | Double_t csigma = (sumx2*norm - cmean*cmean*kfactor); | |
166 | if (csigma<bestsigma || bestsigma<0){ | |
167 | bestmean = cmean; | |
168 | bestsigma = csigma; | |
169 | bestindex = i-hh; | |
170 | } | |
171 | // | |
172 | // | |
173 | sumx += data[index[i]]-data[index[i-hh]]; | |
174 | sumx2 += data[index[i]]*data[index[i]]-data[index[i-hh]]*data[index[i-hh]]; | |
175 | } | |
176 | ||
177 | Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma)); | |
178 | mean = bestmean; | |
179 | sigma = bstd; | |
180 | delete [] index; | |
181 | } | |
182 | ||
183 | ||
184 | //_____________________________________________________________________________ | |
185 | Int_t TStatToolkit::Freq(Int_t n, const Int_t *inlist | |
186 | , Int_t *outlist, Bool_t down) | |
187 | { | |
188 | // | |
189 | // Sort eleements according occurancy | |
190 | // The size of output array has is 2*n | |
191 | // | |
192 | ||
193 | Int_t * sindexS = new Int_t[n]; // temp array for sorting | |
194 | Int_t * sindexF = new Int_t[2*n]; | |
b8072cce | 195 | for (Int_t i=0;i<n;i++) sindexS[i]=0; |
196 | for (Int_t i=0;i<2*n;i++) sindexF[i]=0; | |
21f3a443 | 197 | // |
198 | TMath::Sort(n,inlist, sindexS, down); | |
199 | Int_t last = inlist[sindexS[0]]; | |
200 | Int_t val = last; | |
201 | sindexF[0] = 1; | |
202 | sindexF[0+n] = last; | |
203 | Int_t countPos = 0; | |
204 | // | |
205 | // find frequency | |
206 | for(Int_t i=1;i<n; i++){ | |
207 | val = inlist[sindexS[i]]; | |
208 | if (last == val) sindexF[countPos]++; | |
209 | else{ | |
210 | countPos++; | |
211 | sindexF[countPos+n] = val; | |
212 | sindexF[countPos]++; | |
213 | last =val; | |
214 | } | |
215 | } | |
216 | if (last==val) countPos++; | |
217 | // sort according frequency | |
218 | TMath::Sort(countPos, sindexF, sindexS, kTRUE); | |
219 | for (Int_t i=0;i<countPos;i++){ | |
220 | outlist[2*i ] = sindexF[sindexS[i]+n]; | |
221 | outlist[2*i+1] = sindexF[sindexS[i]]; | |
222 | } | |
223 | delete [] sindexS; | |
224 | delete [] sindexF; | |
225 | ||
226 | return countPos; | |
227 | ||
228 | } | |
229 | ||
230 | //___TStatToolkit__________________________________________________________________________ | |
3d7cc0b4 | 231 | void TStatToolkit::TruncatedMean(const TH1 * his, TVectorD *param, Float_t down, Float_t up, Bool_t verbose){ |
21f3a443 | 232 | // |
233 | // | |
234 | // | |
235 | Int_t nbins = his->GetNbinsX(); | |
236 | 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 | for (Int_t ibin=1;ibin<nbins; ibin++){ | |
242 | ncumul+= his->GetBinContent(ibin); | |
243 | Float_t fraction = Float_t(ncumul)/Float_t(nentries); | |
244 | if (fraction>down && fraction<up){ | |
245 | sum+=his->GetBinContent(ibin); | |
246 | mean+=his->GetBinCenter(ibin)*his->GetBinContent(ibin); | |
247 | sigma2+=his->GetBinCenter(ibin)*his->GetBinCenter(ibin)*his->GetBinContent(ibin); | |
248 | } | |
249 | } | |
250 | mean/=sum; | |
251 | sigma2= TMath::Sqrt(TMath::Abs(sigma2/sum-mean*mean)); | |
252 | if (param){ | |
253 | (*param)[0] = his->GetMaximum(); | |
254 | (*param)[1] = mean; | |
255 | (*param)[2] = sigma2; | |
256 | ||
257 | } | |
258 | if (verbose) printf("Mean\t%f\t Sigma2\t%f\n", mean,sigma2); | |
259 | } | |
260 | ||
261 | void TStatToolkit::LTM(TH1F * his, TVectorD *param , Float_t fraction, Bool_t verbose){ | |
262 | // | |
263 | // LTM | |
264 | // | |
265 | Int_t nbins = his->GetNbinsX(); | |
266 | Int_t nentries = (Int_t)his->GetEntries(); | |
267 | Double_t *data = new Double_t[nentries]; | |
268 | Int_t npoints=0; | |
269 | for (Int_t ibin=1;ibin<nbins; ibin++){ | |
270 | Float_t entriesI = his->GetBinContent(ibin); | |
271 | Float_t xcenter= his->GetBinCenter(ibin); | |
272 | for (Int_t ic=0; ic<entriesI; ic++){ | |
273 | if (npoints<nentries){ | |
274 | data[npoints]= xcenter; | |
275 | npoints++; | |
276 | } | |
277 | } | |
278 | } | |
279 | Double_t mean, sigma; | |
280 | Int_t npoints2=TMath::Min(Int_t(fraction*Float_t(npoints)),npoints-1); | |
281 | npoints2=TMath::Max(Int_t(0.5*Float_t(npoints)),npoints2); | |
282 | TStatToolkit::EvaluateUni(npoints, data, mean,sigma,npoints2); | |
283 | delete [] data; | |
284 | if (verbose) printf("Mean\t%f\t Sigma2\t%f\n", mean,sigma);if (param){ | |
285 | (*param)[0] = his->GetMaximum(); | |
286 | (*param)[1] = mean; | |
287 | (*param)[2] = sigma; | |
288 | } | |
289 | } | |
290 | ||
94a43b22 | 291 | Double_t TStatToolkit::FitGaus(TH1* his, TVectorD *param, TMatrixD */*matrix*/, Float_t xmin, Float_t xmax, Bool_t verbose){ |
21f3a443 | 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 | static TLinearFitter fitter(3,"pol2"); | |
306 | TVectorD par(3); | |
307 | TVectorD sigma(3); | |
308 | TMatrixD mat(3,3); | |
309 | if (his->GetMaximum()<4) return -1; | |
310 | if (his->GetEntries()<12) return -1; | |
311 | if (his->GetRMS()<mat.GetTol()) return -1; | |
312 | Float_t maxEstimate = his->GetEntries()*his->GetBinWidth(1)/TMath::Sqrt((TMath::TwoPi()*his->GetRMS())); | |
313 | Int_t dsmooth = TMath::Nint(6./TMath::Sqrt(maxEstimate)); | |
314 | ||
315 | if (maxEstimate<1) return -1; | |
316 | Int_t nbins = his->GetNbinsX(); | |
317 | Int_t npoints=0; | |
318 | // | |
319 | ||
320 | ||
321 | if (xmin>=xmax){ | |
322 | xmin = his->GetXaxis()->GetXmin(); | |
323 | xmax = his->GetXaxis()->GetXmax(); | |
324 | } | |
325 | for (Int_t iter=0; iter<2; iter++){ | |
326 | fitter.ClearPoints(); | |
327 | npoints=0; | |
328 | for (Int_t ibin=1;ibin<nbins+1; ibin++){ | |
329 | Int_t countB=1; | |
330 | Float_t entriesI = his->GetBinContent(ibin); | |
331 | for (Int_t delta = -dsmooth; delta<=dsmooth; delta++){ | |
332 | if (ibin+delta>1 &&ibin+delta<nbins-1){ | |
333 | entriesI += his->GetBinContent(ibin+delta); | |
334 | countB++; | |
335 | } | |
336 | } | |
337 | entriesI/=countB; | |
338 | Double_t xcenter= his->GetBinCenter(ibin); | |
339 | if (xcenter<xmin || xcenter>xmax) continue; | |
340 | Double_t error=1./TMath::Sqrt(countB); | |
341 | Float_t cont=2; | |
342 | if (iter>0){ | |
343 | if (par[0]+par[1]*xcenter+par[2]*xcenter*xcenter>20) return 0; | |
344 | cont = TMath::Exp(par[0]+par[1]*xcenter+par[2]*xcenter*xcenter); | |
345 | if (cont>1.) error = 1./TMath::Sqrt(cont*Float_t(countB)); | |
346 | } | |
347 | if (entriesI>1&&cont>1){ | |
348 | fitter.AddPoint(&xcenter,TMath::Log(Float_t(entriesI)),error); | |
349 | npoints++; | |
350 | } | |
351 | } | |
352 | if (npoints>3){ | |
353 | fitter.Eval(); | |
354 | fitter.GetParameters(par); | |
355 | }else{ | |
356 | break; | |
357 | } | |
358 | } | |
359 | if (npoints<=3){ | |
360 | return -1; | |
361 | } | |
362 | fitter.GetParameters(par); | |
363 | fitter.GetCovarianceMatrix(mat); | |
364 | if (TMath::Abs(par[1])<mat.GetTol()) return -1; | |
365 | if (TMath::Abs(par[2])<mat.GetTol()) return -1; | |
366 | Double_t chi2 = fitter.GetChisquare()/Float_t(npoints); | |
367 | //fitter.GetParameters(); | |
368 | if (!param) param = new TVectorD(3); | |
cb1d20de | 369 | // if (!matrix) matrix = new TMatrixD(3,3); // Covariance matrix to be implemented |
21f3a443 | 370 | (*param)[1] = par[1]/(-2.*par[2]); |
371 | (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2])); | |
372 | (*param)[0] = TMath::Exp(par[0]+ par[1]* (*param)[1] + par[2]*(*param)[1]*(*param)[1]); | |
373 | if (verbose){ | |
374 | par.Print(); | |
375 | mat.Print(); | |
376 | param->Print(); | |
377 | printf("Chi2=%f\n",chi2); | |
378 | TF1 * f1= new TF1("f1","[0]*exp(-(x-[1])^2/(2*[2]*[2]))",his->GetXaxis()->GetXmin(),his->GetXaxis()->GetXmax()); | |
379 | f1->SetParameter(0, (*param)[0]); | |
380 | f1->SetParameter(1, (*param)[1]); | |
381 | f1->SetParameter(2, (*param)[2]); | |
382 | f1->Draw("same"); | |
383 | } | |
384 | return chi2; | |
385 | } | |
386 | ||
cb1d20de | 387 | Double_t TStatToolkit::FitGaus(Float_t *arr, Int_t nBins, Float_t xMin, Float_t xMax, TVectorD *param, TMatrixD */*matrix*/, Bool_t verbose){ |
21f3a443 | 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) | |
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 | static TLinearFitter fitter(3,"pol2"); | |
409 | static TMatrixD mat(3,3); | |
410 | static Double_t kTol = mat.GetTol(); | |
411 | fitter.StoreData(kFALSE); | |
412 | fitter.ClearPoints(); | |
413 | TVectorD par(3); | |
414 | TVectorD sigma(3); | |
3d7cc0b4 | 415 | TMatrixD matA(3,3); |
21f3a443 | 416 | TMatrixD b(3,1); |
417 | Float_t rms = TMath::RMS(nBins,arr); | |
418 | Float_t max = TMath::MaxElement(nBins,arr); | |
419 | 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 | for (Int_t i=0; i<nBins; i++){ | |
429 | entries+=arr[i]; | |
430 | if (arr[i]>0) nfilled++; | |
431 | } | |
432 | ||
433 | if (max<4) return -4; | |
434 | if (entries<12) return -4; | |
435 | if (rms<kTol) return -4; | |
436 | ||
437 | Int_t npoints=0; | |
438 | // | |
439 | ||
440 | // | |
441 | for (Int_t ibin=0;ibin<nBins; ibin++){ | |
442 | Float_t entriesI = arr[ibin]; | |
443 | if (entriesI>1){ | |
444 | Double_t xcenter = xMin+(ibin+0.5)*binWidth; | |
445 | ||
446 | Float_t error = 1./TMath::Sqrt(entriesI); | |
447 | Float_t val = TMath::Log(Float_t(entriesI)); | |
448 | fitter.AddPoint(&xcenter,val,error); | |
449 | if (npoints<3){ | |
3d7cc0b4 | 450 | matA(npoints,0)=1; |
451 | matA(npoints,1)=xcenter; | |
452 | matA(npoints,2)=xcenter*xcenter; | |
21f3a443 | 453 | b(npoints,0)=val; |
454 | meanCOG+=xcenter*entriesI; | |
455 | rms2COG +=xcenter*entriesI*xcenter; | |
456 | sumCOG +=entriesI; | |
457 | } | |
458 | npoints++; | |
459 | } | |
460 | } | |
461 | ||
462 | ||
463 | Double_t chi2 = 0; | |
464 | if (npoints>=3){ | |
465 | if ( npoints == 3 ){ | |
466 | //analytic calculation of the parameters for three points | |
3d7cc0b4 | 467 | matA.Invert(); |
21f3a443 | 468 | TMatrixD res(1,3); |
3d7cc0b4 | 469 | res.Mult(matA,b); |
21f3a443 | 470 | par[0]=res(0,0); |
471 | par[1]=res(0,1); | |
472 | par[2]=res(0,2); | |
473 | chi2 = -3.; | |
474 | } else { | |
475 | // use fitter for more than three points | |
476 | fitter.Eval(); | |
477 | fitter.GetParameters(par); | |
478 | fitter.GetCovarianceMatrix(mat); | |
479 | chi2 = fitter.GetChisquare()/Float_t(npoints); | |
480 | } | |
481 | if (TMath::Abs(par[1])<kTol) return -4; | |
482 | if (TMath::Abs(par[2])<kTol) return -4; | |
483 | ||
484 | if (!param) param = new TVectorD(3); | |
cb1d20de | 485 | //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 | 486 | |
487 | (*param)[1] = par[1]/(-2.*par[2]); | |
488 | (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2])); | |
489 | Double_t lnparam0 = par[0]+ par[1]* (*param)[1] + par[2]*(*param)[1]*(*param)[1]; | |
490 | if ( lnparam0>307 ) return -4; | |
491 | (*param)[0] = TMath::Exp(lnparam0); | |
492 | if (verbose){ | |
493 | par.Print(); | |
494 | mat.Print(); | |
495 | param->Print(); | |
496 | printf("Chi2=%f\n",chi2); | |
497 | TF1 * f1= new TF1("f1","[0]*exp(-(x-[1])^2/(2*[2]*[2]))",xMin,xMax); | |
498 | f1->SetParameter(0, (*param)[0]); | |
499 | f1->SetParameter(1, (*param)[1]); | |
500 | f1->SetParameter(2, (*param)[2]); | |
501 | f1->Draw("same"); | |
502 | } | |
503 | return chi2; | |
504 | } | |
505 | ||
506 | if (npoints == 2){ | |
507 | //use center of gravity for 2 points | |
508 | meanCOG/=sumCOG; | |
509 | rms2COG /=sumCOG; | |
510 | (*param)[0] = max; | |
511 | (*param)[1] = meanCOG; | |
512 | (*param)[2] = TMath::Sqrt(TMath::Abs(meanCOG*meanCOG-rms2COG)); | |
513 | chi2=-2.; | |
514 | } | |
515 | if ( npoints == 1 ){ | |
516 | meanCOG/=sumCOG; | |
517 | (*param)[0] = max; | |
518 | (*param)[1] = meanCOG; | |
519 | (*param)[2] = binWidth/TMath::Sqrt(12); | |
520 | chi2=-1.; | |
521 | } | |
522 | return chi2; | |
523 | ||
524 | } | |
525 | ||
526 | ||
3d7cc0b4 | 527 | Float_t TStatToolkit::GetCOG(const Short_t *arr, Int_t nBins, Float_t xMin, Float_t xMax, Float_t *rms, Float_t *sum) |
21f3a443 | 528 | { |
529 | // | |
530 | // calculate center of gravity rms and sum for array 'arr' with nBins an a x range xMin to xMax | |
531 | // return COG; in case of failure return xMin | |
532 | // | |
533 | Float_t meanCOG = 0; | |
534 | Float_t rms2COG = 0; | |
535 | Float_t sumCOG = 0; | |
536 | Int_t npoints = 0; | |
537 | ||
538 | Float_t binWidth = (xMax-xMin)/(Float_t)nBins; | |
539 | ||
540 | for (Int_t ibin=0; ibin<nBins; ibin++){ | |
541 | Float_t entriesI = (Float_t)arr[ibin]; | |
542 | Double_t xcenter = xMin+(ibin+0.5)*binWidth; | |
543 | if ( entriesI>0 ){ | |
544 | meanCOG += xcenter*entriesI; | |
545 | rms2COG += xcenter*entriesI*xcenter; | |
546 | sumCOG += entriesI; | |
547 | npoints++; | |
548 | } | |
549 | } | |
550 | if ( sumCOG == 0 ) return xMin; | |
551 | meanCOG/=sumCOG; | |
552 | ||
553 | if ( rms ){ | |
554 | rms2COG /=sumCOG; | |
555 | (*rms) = TMath::Sqrt(TMath::Abs(meanCOG*meanCOG-rms2COG)); | |
556 | if ( npoints == 1 ) (*rms) = binWidth/TMath::Sqrt(12); | |
557 | } | |
558 | ||
559 | if ( sum ) | |
560 | (*sum) = sumCOG; | |
561 | ||
562 | return meanCOG; | |
563 | } | |
564 | ||
565 | ||
566 | ||
567 | /////////////////////////////////////////////////////////////// | |
568 | ////////////// TEST functions ///////////////////////// | |
569 | /////////////////////////////////////////////////////////////// | |
570 | ||
571 | ||
572 | ||
573 | ||
574 | ||
575 | void TStatToolkit::TestGausFit(Int_t nhistos){ | |
576 | // | |
577 | // Test performance of the parabolic - gaussian fit - compare it with | |
578 | // ROOT gauss fit | |
579 | // nhistos - number of histograms to be used for test | |
580 | // | |
581 | TTreeSRedirector *pcstream = new TTreeSRedirector("fitdebug.root"); | |
582 | ||
583 | Float_t *xTrue = new Float_t[nhistos]; | |
584 | Float_t *sTrue = new Float_t[nhistos]; | |
585 | TVectorD **par1 = new TVectorD*[nhistos]; | |
586 | TVectorD **par2 = new TVectorD*[nhistos]; | |
587 | TMatrixD dummy(3,3); | |
588 | ||
589 | ||
590 | TH1F **h1f = new TH1F*[nhistos]; | |
591 | TF1 *myg = new TF1("myg","gaus"); | |
592 | TF1 *fit = new TF1("fit","gaus"); | |
593 | gRandom->SetSeed(0); | |
594 | ||
595 | //init | |
596 | for (Int_t i=0;i<nhistos; i++){ | |
597 | par1[i] = new TVectorD(3); | |
598 | par2[i] = new TVectorD(3); | |
599 | h1f[i] = new TH1F(Form("h1f%d",i),Form("h1f%d",i),20,-10,10); | |
600 | xTrue[i]= gRandom->Rndm(); | |
601 | gSystem->Sleep(2); | |
602 | sTrue[i]= .75+gRandom->Rndm()*.5; | |
603 | myg->SetParameters(1,xTrue[i],sTrue[i]); | |
604 | h1f[i]->FillRandom("myg"); | |
605 | } | |
606 | ||
607 | TStopwatch s; | |
608 | s.Start(); | |
609 | //standard gaus fit | |
610 | for (Int_t i=0; i<nhistos; i++){ | |
611 | h1f[i]->Fit(fit,"0q"); | |
612 | (*par1[i])(0) = fit->GetParameter(0); | |
613 | (*par1[i])(1) = fit->GetParameter(1); | |
614 | (*par1[i])(2) = fit->GetParameter(2); | |
615 | } | |
616 | s.Stop(); | |
617 | printf("Gaussian fit\t"); | |
618 | s.Print(); | |
619 | ||
620 | s.Start(); | |
621 | //TStatToolkit gaus fit | |
622 | for (Int_t i=0; i<nhistos; i++){ | |
623 | TStatToolkit::FitGaus(h1f[i]->GetArray()+1,h1f[i]->GetNbinsX(),h1f[i]->GetXaxis()->GetXmin(),h1f[i]->GetXaxis()->GetXmax(),par2[i],&dummy); | |
624 | } | |
625 | ||
626 | s.Stop(); | |
627 | printf("Parabolic fit\t"); | |
628 | s.Print(); | |
629 | //write stream | |
630 | for (Int_t i=0;i<nhistos; i++){ | |
631 | Float_t xt = xTrue[i]; | |
632 | Float_t st = sTrue[i]; | |
633 | (*pcstream)<<"data" | |
634 | <<"xTrue="<<xt | |
635 | <<"sTrue="<<st | |
636 | <<"pg.="<<(par1[i]) | |
637 | <<"pa.="<<(par2[i]) | |
638 | <<"\n"; | |
639 | } | |
640 | //delete pointers | |
641 | for (Int_t i=0;i<nhistos; i++){ | |
642 | delete par1[i]; | |
643 | delete par2[i]; | |
644 | delete h1f[i]; | |
645 | } | |
646 | delete pcstream; | |
647 | delete []h1f; | |
648 | delete []xTrue; | |
649 | delete []sTrue; | |
650 | // | |
651 | delete []par1; | |
652 | delete []par2; | |
653 | ||
654 | } | |
655 | ||
656 | ||
657 | ||
658 | TGraph2D * TStatToolkit::MakeStat2D(TH3 * his, Int_t delta0, Int_t delta1, Int_t type){ | |
659 | // | |
660 | // | |
661 | // | |
662 | // delta - number of bins to integrate | |
663 | // type - 0 - mean value | |
664 | ||
665 | TAxis * xaxis = his->GetXaxis(); | |
666 | TAxis * yaxis = his->GetYaxis(); | |
667 | // TAxis * zaxis = his->GetZaxis(); | |
668 | Int_t nbinx = xaxis->GetNbins(); | |
669 | Int_t nbiny = yaxis->GetNbins(); | |
670 | char name[1000]; | |
671 | Int_t icount=0; | |
672 | TGraph2D *graph = new TGraph2D(nbinx*nbiny); | |
673 | TF1 f1("f1","gaus"); | |
674 | for (Int_t ix=0; ix<nbinx;ix++) | |
675 | for (Int_t iy=0; iy<nbiny;iy++){ | |
676 | Float_t xcenter = xaxis->GetBinCenter(ix); | |
677 | Float_t ycenter = yaxis->GetBinCenter(iy); | |
cb1d20de | 678 | snprintf(name,1000,"%s_%d_%d",his->GetName(), ix,iy); |
21f3a443 | 679 | TH1 *projection = his->ProjectionZ(name,ix-delta0,ix+delta0,iy-delta1,iy+delta1); |
680 | Float_t stat= 0; | |
681 | if (type==0) stat = projection->GetMean(); | |
682 | if (type==1) stat = projection->GetRMS(); | |
683 | if (type==2 || type==3){ | |
684 | TVectorD vec(3); | |
685 | TStatToolkit::LTM((TH1F*)projection,&vec,0.7); | |
686 | if (type==2) stat= vec[1]; | |
687 | if (type==3) stat= vec[0]; | |
688 | } | |
689 | if (type==4|| type==5){ | |
690 | projection->Fit(&f1); | |
691 | if (type==4) stat= f1.GetParameter(1); | |
692 | if (type==5) stat= f1.GetParameter(2); | |
693 | } | |
694 | //printf("%d\t%f\t%f\t%f\n", icount,xcenter, ycenter, stat); | |
695 | graph->SetPoint(icount,xcenter, ycenter, stat); | |
696 | icount++; | |
697 | } | |
698 | return graph; | |
699 | } | |
700 | ||
701 | TGraph * TStatToolkit::MakeStat1D(TH3 * his, Int_t delta1, Int_t type){ | |
702 | // | |
703 | // | |
704 | // | |
705 | // delta - number of bins to integrate | |
706 | // type - 0 - mean value | |
707 | ||
708 | TAxis * xaxis = his->GetXaxis(); | |
709 | TAxis * yaxis = his->GetYaxis(); | |
710 | // TAxis * zaxis = his->GetZaxis(); | |
711 | Int_t nbinx = xaxis->GetNbins(); | |
712 | Int_t nbiny = yaxis->GetNbins(); | |
713 | char name[1000]; | |
714 | Int_t icount=0; | |
715 | TGraph *graph = new TGraph(nbinx); | |
716 | TF1 f1("f1","gaus"); | |
717 | for (Int_t ix=0; ix<nbinx;ix++){ | |
718 | Float_t xcenter = xaxis->GetBinCenter(ix); | |
719 | // Float_t ycenter = yaxis->GetBinCenter(iy); | |
cb1d20de | 720 | snprintf(name,1000,"%s_%d",his->GetName(), ix); |
21f3a443 | 721 | TH1 *projection = his->ProjectionZ(name,ix-delta1,ix+delta1,0,nbiny); |
722 | Float_t stat= 0; | |
723 | if (type==0) stat = projection->GetMean(); | |
724 | if (type==1) stat = projection->GetRMS(); | |
725 | if (type==2 || type==3){ | |
726 | TVectorD vec(3); | |
727 | TStatToolkit::LTM((TH1F*)projection,&vec,0.7); | |
728 | if (type==2) stat= vec[1]; | |
729 | if (type==3) stat= vec[0]; | |
730 | } | |
731 | if (type==4|| type==5){ | |
732 | projection->Fit(&f1); | |
733 | if (type==4) stat= f1.GetParameter(1); | |
734 | if (type==5) stat= f1.GetParameter(2); | |
735 | } | |
736 | //printf("%d\t%f\t%f\t%f\n", icount,xcenter, ycenter, stat); | |
737 | graph->SetPoint(icount,xcenter, stat); | |
738 | icount++; | |
739 | } | |
740 | return graph; | |
741 | } | |
742 | ||
743 | ||
744 | ||
745 | ||
746 | ||
88b1c775 | 747 | 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 | 748 | // |
749 | // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts | |
750 | // returns chi2, fitParam and covMatrix | |
751 | // returns TString with fitted formula | |
752 | // | |
dd46129c | 753 | |
21f3a443 | 754 | TString formulaStr(formula); |
755 | TString drawStr(drawCommand); | |
756 | TString cutStr(cuts); | |
dd46129c | 757 | TString ferr("1"); |
758 | ||
759 | TString strVal(drawCommand); | |
760 | if (strVal.Contains(":")){ | |
761 | TObjArray* valTokens = strVal.Tokenize(":"); | |
762 | drawStr = valTokens->At(0)->GetName(); | |
763 | ferr = valTokens->At(1)->GetName(); | |
09d5920f | 764 | delete valTokens; |
dd46129c | 765 | } |
766 | ||
21f3a443 | 767 | |
768 | formulaStr.ReplaceAll("++", "~"); | |
769 | TObjArray* formulaTokens = formulaStr.Tokenize("~"); | |
770 | Int_t dim = formulaTokens->GetEntriesFast(); | |
771 | ||
772 | fitParam.ResizeTo(dim); | |
773 | covMatrix.ResizeTo(dim,dim); | |
774 | ||
775 | TLinearFitter* fitter = new TLinearFitter(dim+1, Form("hyp%d",dim)); | |
776 | fitter->StoreData(kTRUE); | |
777 | fitter->ClearPoints(); | |
778 | ||
779 | Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start, start); | |
09d5920f | 780 | if (entries == -1) { |
781 | delete formulaTokens; | |
782 | return new TString("An ERROR has occured during fitting!"); | |
783 | } | |
bd7b4d18 | 784 | Double_t **values = new Double_t*[dim+1] ; |
785 | for (Int_t i=0; i<dim+1; i++) values[i]=NULL; | |
dd46129c | 786 | // |
787 | entries = tree->Draw(ferr.Data(), cutStr.Data(), "goff", stop-start, start); | |
b8072cce | 788 | if (entries == -1) { |
09d5920f | 789 | delete formulaTokens; |
b8072cce | 790 | delete []values; |
791 | return new TString("An ERROR has occured during fitting!"); | |
792 | } | |
dd46129c | 793 | Double_t *errors = new Double_t[entries]; |
794 | memcpy(errors, tree->GetV1(), entries*sizeof(Double_t)); | |
21f3a443 | 795 | |
796 | for (Int_t i = 0; i < dim + 1; i++){ | |
797 | Int_t centries = 0; | |
798 | if (i < dim) centries = tree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start); | |
799 | else centries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start); | |
800 | ||
b8072cce | 801 | if (entries != centries) { |
802 | delete []errors; | |
803 | delete []values; | |
804 | return new TString("An ERROR has occured during fitting!"); | |
805 | } | |
21f3a443 | 806 | values[i] = new Double_t[entries]; |
807 | memcpy(values[i], tree->GetV1(), entries*sizeof(Double_t)); | |
808 | } | |
809 | ||
810 | // add points to the fitter | |
811 | for (Int_t i = 0; i < entries; i++){ | |
812 | Double_t x[1000]; | |
813 | for (Int_t j=0; j<dim;j++) x[j]=values[j][i]; | |
dd46129c | 814 | fitter->AddPoint(x, values[dim][i], errors[i]); |
21f3a443 | 815 | } |
816 | ||
817 | fitter->Eval(); | |
2c629c56 | 818 | if (frac>0.5 && frac<1){ |
819 | fitter->EvalRobust(frac); | |
88b1c775 | 820 | }else{ |
821 | if (fix0) { | |
822 | fitter->FixParameter(0,0); | |
823 | fitter->Eval(); | |
824 | } | |
2c629c56 | 825 | } |
21f3a443 | 826 | fitter->GetParameters(fitParam); |
827 | fitter->GetCovarianceMatrix(covMatrix); | |
828 | chi2 = fitter->GetChisquare(); | |
b8072cce | 829 | npoints = entries; |
21f3a443 | 830 | TString *preturnFormula = new TString(Form("( %f+",fitParam[0])), &returnFormula = *preturnFormula; |
831 | ||
832 | for (Int_t iparam = 0; iparam < dim; iparam++) { | |
833 | returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam+1])); | |
834 | if (iparam < dim-1) returnFormula.Append("+"); | |
835 | } | |
836 | returnFormula.Append(" )"); | |
4d61c301 | 837 | |
838 | ||
b8072cce | 839 | for (Int_t j=0; j<dim+1;j++) delete [] values[j]; |
4d61c301 | 840 | |
841 | ||
cb1d20de | 842 | delete formulaTokens; |
843 | delete fitter; | |
844 | delete[] values; | |
b8072cce | 845 | delete[] errors; |
cb1d20de | 846 | return preturnFormula; |
847 | } | |
848 | ||
849 | 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){ | |
850 | // | |
851 | // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts | |
852 | // returns chi2, fitParam and covMatrix | |
853 | // returns TString with fitted formula | |
854 | // | |
855 | ||
856 | TString formulaStr(formula); | |
857 | TString drawStr(drawCommand); | |
858 | TString cutStr(cuts); | |
859 | TString ferr("1"); | |
860 | ||
861 | TString strVal(drawCommand); | |
862 | if (strVal.Contains(":")){ | |
863 | TObjArray* valTokens = strVal.Tokenize(":"); | |
864 | drawStr = valTokens->At(0)->GetName(); | |
865 | ferr = valTokens->At(1)->GetName(); | |
09d5920f | 866 | delete valTokens; |
cb1d20de | 867 | } |
868 | ||
869 | ||
870 | formulaStr.ReplaceAll("++", "~"); | |
871 | TObjArray* formulaTokens = formulaStr.Tokenize("~"); | |
872 | Int_t dim = formulaTokens->GetEntriesFast(); | |
873 | ||
874 | fitParam.ResizeTo(dim); | |
875 | covMatrix.ResizeTo(dim,dim); | |
876 | ||
877 | TLinearFitter* fitter = new TLinearFitter(dim+1, Form("hyp%d",dim)); | |
878 | fitter->StoreData(kTRUE); | |
879 | fitter->ClearPoints(); | |
880 | ||
881 | Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start, start); | |
09d5920f | 882 | if (entries == -1) { |
883 | delete formulaTokens; | |
884 | return new TString("An ERROR has occured during fitting!"); | |
885 | } | |
cb1d20de | 886 | Double_t **values = new Double_t*[dim+1] ; |
bd7b4d18 | 887 | for (Int_t i=0; i<dim+1; i++) values[i]=NULL; |
cb1d20de | 888 | // |
889 | entries = tree->Draw(ferr.Data(), cutStr.Data(), "goff", stop-start, start); | |
b8072cce | 890 | if (entries == -1) { |
09d5920f | 891 | delete formulaTokens; |
b8072cce | 892 | delete [] values; |
893 | return new TString("An ERROR has occured during fitting!"); | |
894 | } | |
cb1d20de | 895 | Double_t *errors = new Double_t[entries]; |
896 | memcpy(errors, tree->GetV1(), entries*sizeof(Double_t)); | |
897 | ||
898 | for (Int_t i = 0; i < dim + 1; i++){ | |
899 | Int_t centries = 0; | |
900 | if (i < dim) centries = tree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start); | |
901 | else centries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start); | |
902 | ||
b8072cce | 903 | if (entries != centries) { |
904 | delete []errors; | |
905 | delete []values; | |
09d5920f | 906 | delete formulaTokens; |
b8072cce | 907 | return new TString("An ERROR has occured during fitting!"); |
908 | } | |
cb1d20de | 909 | values[i] = new Double_t[entries]; |
910 | memcpy(values[i], tree->GetV1(), entries*sizeof(Double_t)); | |
911 | } | |
912 | ||
913 | // add points to the fitter | |
914 | for (Int_t i = 0; i < entries; i++){ | |
915 | Double_t x[1000]; | |
916 | for (Int_t j=0; j<dim;j++) x[j]=values[j][i]; | |
917 | fitter->AddPoint(x, values[dim][i], errors[i]); | |
918 | } | |
919 | if (constrain>0){ | |
920 | for (Int_t i = 0; i < dim; i++){ | |
921 | Double_t x[1000]; | |
922 | for (Int_t j=0; j<dim;j++) if (i!=j) x[j]=0; | |
923 | x[i]=1.; | |
924 | fitter->AddPoint(x, 0, constrain); | |
925 | } | |
926 | } | |
927 | ||
928 | ||
929 | fitter->Eval(); | |
930 | if (frac>0.5 && frac<1){ | |
931 | fitter->EvalRobust(frac); | |
932 | } | |
933 | fitter->GetParameters(fitParam); | |
934 | fitter->GetCovarianceMatrix(covMatrix); | |
935 | chi2 = fitter->GetChisquare(); | |
936 | npoints = entries; | |
cb1d20de | 937 | |
938 | TString *preturnFormula = new TString(Form("( %f+",fitParam[0])), &returnFormula = *preturnFormula; | |
939 | ||
940 | for (Int_t iparam = 0; iparam < dim; iparam++) { | |
941 | returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam+1])); | |
942 | if (iparam < dim-1) returnFormula.Append("+"); | |
943 | } | |
944 | returnFormula.Append(" )"); | |
945 | ||
b8072cce | 946 | for (Int_t j=0; j<dim+1;j++) delete [] values[j]; |
cb1d20de | 947 | |
948 | ||
949 | ||
950 | delete formulaTokens; | |
951 | delete fitter; | |
952 | delete[] values; | |
b8072cce | 953 | delete[] errors; |
cb1d20de | 954 | return preturnFormula; |
955 | } | |
956 | ||
957 | ||
958 | ||
959 | 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){ | |
960 | // | |
961 | // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts | |
962 | // returns chi2, fitParam and covMatrix | |
963 | // returns TString with fitted formula | |
964 | // | |
965 | ||
966 | TString formulaStr(formula); | |
967 | TString drawStr(drawCommand); | |
968 | TString cutStr(cuts); | |
969 | TString ferr("1"); | |
970 | ||
971 | TString strVal(drawCommand); | |
972 | if (strVal.Contains(":")){ | |
973 | TObjArray* valTokens = strVal.Tokenize(":"); | |
974 | drawStr = valTokens->At(0)->GetName(); | |
09d5920f | 975 | ferr = valTokens->At(1)->GetName(); |
976 | delete valTokens; | |
cb1d20de | 977 | } |
978 | ||
979 | ||
980 | formulaStr.ReplaceAll("++", "~"); | |
981 | TObjArray* formulaTokens = formulaStr.Tokenize("~"); | |
982 | Int_t dim = formulaTokens->GetEntriesFast(); | |
983 | ||
984 | fitParam.ResizeTo(dim); | |
985 | covMatrix.ResizeTo(dim,dim); | |
986 | TString fitString="x0"; | |
987 | for (Int_t i=1; i<dim; i++) fitString+=Form("++x%d",i); | |
988 | TLinearFitter* fitter = new TLinearFitter(dim, fitString.Data()); | |
989 | fitter->StoreData(kTRUE); | |
990 | fitter->ClearPoints(); | |
991 | ||
992 | Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start, start); | |
09d5920f | 993 | if (entries == -1) { |
994 | delete formulaTokens; | |
995 | return new TString("An ERROR has occured during fitting!"); | |
996 | } | |
cb1d20de | 997 | Double_t **values = new Double_t*[dim+1] ; |
bd7b4d18 | 998 | for (Int_t i=0; i<dim+1; i++) values[i]=NULL; |
cb1d20de | 999 | // |
1000 | entries = tree->Draw(ferr.Data(), cutStr.Data(), "goff", stop-start, start); | |
b8072cce | 1001 | if (entries == -1) { |
1002 | delete []values; | |
09d5920f | 1003 | delete formulaTokens; |
b8072cce | 1004 | return new TString("An ERROR has occured during fitting!"); |
1005 | } | |
cb1d20de | 1006 | Double_t *errors = new Double_t[entries]; |
1007 | memcpy(errors, tree->GetV1(), entries*sizeof(Double_t)); | |
1008 | ||
1009 | for (Int_t i = 0; i < dim + 1; i++){ | |
1010 | Int_t centries = 0; | |
1011 | if (i < dim) centries = tree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start); | |
1012 | else centries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start); | |
1013 | ||
b8072cce | 1014 | if (entries != centries) { |
1015 | delete []errors; | |
1016 | delete []values; | |
09d5920f | 1017 | delete formulaTokens; |
b8072cce | 1018 | return new TString("An ERROR has occured during fitting!"); |
1019 | } | |
cb1d20de | 1020 | values[i] = new Double_t[entries]; |
1021 | memcpy(values[i], tree->GetV1(), entries*sizeof(Double_t)); | |
1022 | } | |
1023 | ||
1024 | // add points to the fitter | |
1025 | for (Int_t i = 0; i < entries; i++){ | |
1026 | Double_t x[1000]; | |
1027 | for (Int_t j=0; j<dim;j++) x[j]=values[j][i]; | |
1028 | fitter->AddPoint(x, values[dim][i], errors[i]); | |
1029 | } | |
1030 | ||
1031 | fitter->Eval(); | |
1032 | if (frac>0.5 && frac<1){ | |
1033 | fitter->EvalRobust(frac); | |
1034 | } | |
1035 | fitter->GetParameters(fitParam); | |
1036 | fitter->GetCovarianceMatrix(covMatrix); | |
1037 | chi2 = fitter->GetChisquare(); | |
1038 | npoints = entries; | |
cb1d20de | 1039 | |
1040 | TString *preturnFormula = new TString("("), &returnFormula = *preturnFormula; | |
1041 | ||
1042 | for (Int_t iparam = 0; iparam < dim; iparam++) { | |
1043 | returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam])); | |
1044 | if (iparam < dim-1) returnFormula.Append("+"); | |
1045 | } | |
1046 | returnFormula.Append(" )"); | |
1047 | ||
1048 | ||
b8072cce | 1049 | for (Int_t j=0; j<dim+1;j++) delete [] values[j]; |
cb1d20de | 1050 | |
21f3a443 | 1051 | delete formulaTokens; |
1052 | delete fitter; | |
1053 | delete[] values; | |
b8072cce | 1054 | delete[] errors; |
21f3a443 | 1055 | return preturnFormula; |
1056 | } | |
7c9cf6e4 | 1057 | |
1058 | ||
1059 | ||
1060 | ||
1061 | ||
3d7cc0b4 | 1062 | Int_t TStatToolkit::GetFitIndex(const TString fString, const TString subString){ |
7c9cf6e4 | 1063 | // |
1064 | // fitString - ++ separated list of fits | |
1065 | // substring - ++ separated list of the requiered substrings | |
1066 | // | |
1067 | // return the last occurance of substring in fit string | |
1068 | // | |
1069 | TObjArray *arrFit = fString.Tokenize("++"); | |
1070 | TObjArray *arrSub = subString.Tokenize("++"); | |
1071 | Int_t index=-1; | |
1072 | for (Int_t i=0; i<arrFit->GetEntries(); i++){ | |
1073 | Bool_t isOK=kTRUE; | |
1074 | TString str =arrFit->At(i)->GetName(); | |
1075 | for (Int_t isub=0; isub<arrSub->GetEntries(); isub++){ | |
1076 | if (str.Contains(arrSub->At(isub)->GetName())==0) isOK=kFALSE; | |
1077 | } | |
1078 | if (isOK) index=i; | |
1079 | } | |
09d5920f | 1080 | delete arrFit; |
1081 | delete arrSub; | |
7c9cf6e4 | 1082 | return index; |
1083 | } | |
1084 | ||
1085 | ||
3d7cc0b4 | 1086 | TString TStatToolkit::FilterFit(const TString &input, const TString filter, TVectorD ¶m, TMatrixD & covar){ |
7c9cf6e4 | 1087 | // |
1088 | // Filter fit expression make sub-fit | |
1089 | // | |
1090 | TObjArray *array0= input.Tokenize("++"); | |
1091 | TObjArray *array1= filter.Tokenize("++"); | |
1092 | //TString *presult=new TString("(0"); | |
1093 | TString result="(0.0"; | |
1094 | for (Int_t i=0; i<array0->GetEntries(); i++){ | |
1095 | Bool_t isOK=kTRUE; | |
1096 | TString str(array0->At(i)->GetName()); | |
1097 | for (Int_t j=0; j<array1->GetEntries(); j++){ | |
1098 | if (str.Contains(array1->At(j)->GetName())==0) isOK=kFALSE; | |
1099 | } | |
1100 | if (isOK) { | |
1101 | result+="+"+str; | |
1102 | result+=Form("*(%f)",param[i+1]); | |
1103 | printf("%f\t%f\t%s\n",param[i+1], TMath::Sqrt(covar(i+1,i+1)),str.Data()); | |
1104 | } | |
1105 | } | |
1106 | result+="-0.)"; | |
09d5920f | 1107 | delete array0; |
1108 | delete array1; | |
7c9cf6e4 | 1109 | return result; |
1110 | } | |
1111 | ||
1112 | void TStatToolkit::Update1D(Double_t delta, Double_t sigma, Int_t s1, TMatrixD &vecXk, TMatrixD &covXk){ | |
1113 | // | |
1114 | // Update parameters and covariance - with one measurement | |
1115 | // Input: | |
1116 | // vecXk - input vector - Updated in function | |
1117 | // covXk - covariance matrix - Updated in function | |
1118 | // delta, sigma, s1 - new measurement, rms of new measurement and the index of measurement | |
1119 | const Int_t knMeas=1; | |
1120 | Int_t knElem=vecXk.GetNrows(); | |
1121 | ||
1122 | TMatrixD mat1(knElem,knElem); // update covariance matrix | |
1123 | TMatrixD matHk(1,knElem); // vector to mesurement | |
1124 | TMatrixD vecYk(knMeas,1); // Innovation or measurement residual | |
1125 | TMatrixD matHkT(knElem,knMeas); // helper matrix Hk transpose | |
1126 | TMatrixD matSk(knMeas,knMeas); // Innovation (or residual) covariance | |
1127 | TMatrixD matKk(knElem,knMeas); // Optimal Kalman gain | |
1128 | TMatrixD covXk2(knElem,knElem); // helper matrix | |
1129 | TMatrixD covXk3(knElem,knElem); // helper matrix | |
1130 | TMatrixD vecZk(1,1); | |
1131 | TMatrixD measR(1,1); | |
1132 | vecZk(0,0)=delta; | |
1133 | measR(0,0)=sigma*sigma; | |
1134 | // | |
1135 | // reset matHk | |
1136 | for (Int_t iel=0;iel<knElem;iel++) | |
1137 | for (Int_t ip=0;ip<knMeas;ip++) matHk(ip,iel)=0; | |
1138 | //mat1 | |
1139 | for (Int_t iel=0;iel<knElem;iel++) { | |
1140 | for (Int_t jel=0;jel<knElem;jel++) mat1(iel,jel)=0; | |
1141 | mat1(iel,iel)=1; | |
1142 | } | |
1143 | // | |
1144 | matHk(0, s1)=1; | |
1145 | vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual | |
1146 | matHkT=matHk.T(); matHk.T(); | |
1147 | matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance | |
1148 | matSk.Invert(); | |
1149 | matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain | |
1150 | vecXk += matKk*vecYk; // updated vector | |
1151 | covXk2= (mat1-(matKk*matHk)); | |
1152 | covXk3 = covXk2*covXk; | |
1153 | covXk = covXk3; | |
1154 | Int_t nrows=covXk3.GetNrows(); | |
1155 | ||
1156 | for (Int_t irow=0; irow<nrows; irow++) | |
1157 | for (Int_t icol=0; icol<nrows; icol++){ | |
1158 | // rounding problems - make matrix again symteric | |
1159 | covXk(irow,icol)=(covXk3(irow,icol)+covXk3(icol,irow))*0.5; | |
1160 | } | |
1161 | } | |
1162 | ||
1163 | ||
1164 | ||
3d7cc0b4 | 1165 | void TStatToolkit::Constrain1D(const TString &input, const TString filter, TVectorD ¶m, TMatrixD & covar, Double_t mean, Double_t sigma){ |
7c9cf6e4 | 1166 | // |
1167 | // constrain linear fit | |
1168 | // input - string description of fit function | |
1169 | // filter - string filter to select sub fits | |
1170 | // param,covar - parameters and covariance matrix of the fit | |
1171 | // mean,sigma - new measurement uning which the fit is updated | |
1172 | // | |
ae45c94d | 1173 | |
7c9cf6e4 | 1174 | TObjArray *array0= input.Tokenize("++"); |
1175 | TObjArray *array1= filter.Tokenize("++"); | |
1176 | TMatrixD paramM(param.GetNrows(),1); | |
1177 | for (Int_t i=0; i<=array0->GetEntries(); i++){paramM(i,0)=param(i);} | |
1178 | ||
ae45c94d | 1179 | if (filter.Length()==0){ |
1180 | TStatToolkit::Update1D(mean, sigma, 0, paramM, covar);// | |
1181 | }else{ | |
1182 | for (Int_t i=0; i<array0->GetEntries(); i++){ | |
1183 | Bool_t isOK=kTRUE; | |
1184 | TString str(array0->At(i)->GetName()); | |
1185 | for (Int_t j=0; j<array1->GetEntries(); j++){ | |
1186 | if (str.Contains(array1->At(j)->GetName())==0) isOK=kFALSE; | |
1187 | } | |
1188 | if (isOK) { | |
1189 | TStatToolkit::Update1D(mean, sigma, i+1, paramM, covar);// | |
1190 | } | |
7c9cf6e4 | 1191 | } |
1192 | } | |
1193 | for (Int_t i=0; i<=array0->GetEntries(); i++){ | |
1194 | param(i)=paramM(i,0); | |
1195 | } | |
09d5920f | 1196 | delete array0; |
1197 | delete array1; | |
7c9cf6e4 | 1198 | } |
1199 | ||
ae45c94d | 1200 | TString TStatToolkit::MakeFitString(const TString &input, const TVectorD ¶m, const TMatrixD & covar, Bool_t verbose){ |
7c9cf6e4 | 1201 | // |
1202 | // | |
1203 | // | |
1204 | TObjArray *array0= input.Tokenize("++"); | |
ae45c94d | 1205 | TString result=Form("(%f",param[0]); |
1206 | printf("%f\t%f\t\n", param[0], TMath::Sqrt(covar(0,0))); | |
7c9cf6e4 | 1207 | for (Int_t i=0; i<array0->GetEntries(); i++){ |
1208 | TString str(array0->At(i)->GetName()); | |
1209 | result+="+"+str; | |
1210 | result+=Form("*(%f)",param[i+1]); | |
ae45c94d | 1211 | if (verbose) printf("%f\t%f\t%s\n", param[i+1], TMath::Sqrt(covar(i+1,i+1)),str.Data()); |
7c9cf6e4 | 1212 | } |
1213 | result+="-0.)"; | |
09d5920f | 1214 | delete array0; |
7c9cf6e4 | 1215 | return result; |
1216 | } | |
df0a2a0a | 1217 | |
1a8f4649 | 1218 | TGraphErrors * TStatToolkit::MakeGraphErrors(TTree * tree, const char * expr, const char * cut, Int_t mstyle, Int_t mcolor, Float_t msize, Float_t offset){ |
1219 | // | |
1220 | // Query a graph errors | |
1221 | // return TGraphErrors specified by expr and cut | |
1222 | // Example usage TStatToolkit::MakeGraphError(tree,"Y:X:ErrY","X>0", 25,2,0.4) | |
1223 | // tree - tree with variable | |
1224 | // expr - examp | |
1225 | const Int_t entries = tree->Draw(expr,cut,"goff"); | |
1226 | if (entries<=0) { | |
1227 | TStatToolkit t; | |
1228 | t.Error("TStatToolkit::MakeGraphError",Form("Empty or Not valid expression (%s) or cut *%s)", expr,cut)); | |
1229 | return 0; | |
1230 | } | |
1231 | if ( tree->GetV2()==0){ | |
1232 | TStatToolkit t; | |
1233 | t.Error("TStatToolkit::MakeGraphError",Form("Not valid expression (%s) ", expr)); | |
1234 | return 0; | |
1235 | } | |
1236 | TGraphErrors * graph=0; | |
1237 | if ( tree->GetV3()!=0){ | |
1238 | graph = new TGraphErrors (entries, tree->GetV2(),tree->GetV1(),0,tree->GetV3()); | |
1239 | }else{ | |
1240 | graph = new TGraphErrors (entries, tree->GetV2(),tree->GetV1(),0,0); | |
1241 | } | |
1242 | graph->SetMarkerStyle(mstyle); | |
1243 | graph->SetMarkerColor(mcolor); | |
1244 | graph->SetLineColor(mcolor); | |
1245 | if (msize>0) graph->SetMarkerSize(msize); | |
1246 | for(Int_t i=0;i<graph->GetN();i++) graph->GetX()[i]+=offset; | |
1247 | return graph; | |
1248 | ||
1249 | } | |
1250 | ||
df0a2a0a | 1251 | |
377a7d60 | 1252 | TGraph * TStatToolkit::MakeGraphSparse(TTree * tree, const char * expr, const char * cut, Int_t mstyle, Int_t mcolor, Float_t msize, Float_t offset){ |
df0a2a0a | 1253 | // |
1254 | // Make a sparse draw of the variables | |
d02b8756 | 1255 | // Format of expr : Var:Run or Var:Run:ErrorY or Var:Run:ErrorY:ErrorX |
1256 | // offset : points can slightly be shifted in x for better visibility with more graphs | |
1257 | // | |
1258 | // Written by Weilin.Yu | |
1259 | // updated & merged with QA-code by Patrick Reichelt | |
1260 | // | |
1261 | const Int_t entries = tree->Draw(expr,cut,"goff"); | |
8fb3bea1 | 1262 | if (entries<=0) { |
1263 | TStatToolkit t; | |
d02b8756 | 1264 | t.Error("TStatToolkit::MakeGraphSparse",Form("Empty or Not valid expression (%s) or cut (%s)", expr, cut)); |
8fb3bea1 | 1265 | return 0; |
1266 | } | |
df0a2a0a | 1267 | // TGraph * graph = (TGraph*)gPad->GetPrimitive("Graph"); // 2D |
d02b8756 | 1268 | |
1269 | Double_t *graphY, *graphX; | |
1270 | graphY = tree->GetV1(); | |
1271 | graphX = tree->GetV2(); | |
1272 | ||
1273 | // sort according to run number | |
eda18694 | 1274 | Int_t *index = new Int_t[entries*4]; |
d02b8756 | 1275 | TMath::Sort(entries,graphX,index,kFALSE); |
df0a2a0a | 1276 | |
d02b8756 | 1277 | // define arrays for the new graph |
1278 | Double_t *unsortedX = new Double_t[entries]; | |
1279 | Int_t *runNumber = new Int_t[entries]; | |
df0a2a0a | 1280 | Double_t count = 0.5; |
d02b8756 | 1281 | |
1282 | // evaluate arrays for the new graph according to the run-number | |
baa0041d | 1283 | Int_t icount=0; |
d02b8756 | 1284 | //first entry |
1285 | unsortedX[index[0]] = count; | |
1286 | runNumber[0] = graphX[index[0]]; | |
1287 | // loop the rest of entries | |
1288 | for(Int_t i=1;i<entries;i++) | |
1289 | { | |
1290 | if(graphX[index[i]]==graphX[index[i-1]]) | |
1291 | unsortedX[index[i]] = count; | |
1292 | else if(graphX[index[i]]!=graphX[index[i-1]]){ | |
df0a2a0a | 1293 | count++; |
baa0041d | 1294 | icount++; |
d02b8756 | 1295 | unsortedX[index[i]] = count; |
1296 | runNumber[icount]=graphX[index[i]]; | |
df0a2a0a | 1297 | } |
1298 | } | |
d02b8756 | 1299 | |
1300 | // count the number of xbins (run-wise) for the new graph | |
df0a2a0a | 1301 | const Int_t newNbins = int(count+0.5); |
1302 | Double_t *newBins = new Double_t[newNbins+1]; | |
1303 | for(Int_t i=0; i<=count+1;i++){ | |
1304 | newBins[i] = i; | |
1305 | } | |
d02b8756 | 1306 | |
1307 | // define and fill the new graph | |
b3453fe7 | 1308 | TGraph *graphNew = 0; |
d02b8756 | 1309 | if (tree->GetV3()) { |
1310 | if (tree->GetV4()) { | |
1311 | graphNew = new TGraphErrors(entries,unsortedX,graphY,tree->GetV4(),tree->GetV3()); | |
1312 | } | |
1313 | else { graphNew = new TGraphErrors(entries,unsortedX,graphY,0,tree->GetV3()); } | |
1314 | } | |
1315 | else { graphNew = new TGraphErrors(entries,unsortedX,graphY,0,0); } | |
1316 | // with "Set(...)", the x-axis is being sorted | |
df0a2a0a | 1317 | graphNew->GetXaxis()->Set(newNbins,newBins); |
d02b8756 | 1318 | |
1319 | // set the bins for the x-axis, apply shifting of points | |
df0a2a0a | 1320 | Char_t xName[50]; |
df0a2a0a | 1321 | for(Int_t i=0;i<count;i++){ |
d02b8756 | 1322 | snprintf(xName,50,"%d",runNumber[i]); |
df0a2a0a | 1323 | graphNew->GetXaxis()->SetBinLabel(i+1,xName); |
d02b8756 | 1324 | graphNew->GetX()[i]+=offset; |
df0a2a0a | 1325 | } |
d02b8756 | 1326 | |
df0a2a0a | 1327 | graphNew->GetHistogram()->SetTitle(""); |
d02b8756 | 1328 | graphNew->SetMarkerStyle(mstyle); |
b3453fe7 | 1329 | graphNew->SetMarkerColor(mcolor); |
70989f8d | 1330 | if (msize>0) graphNew->SetMarkerSize(msize); |
d02b8756 | 1331 | delete [] unsortedX; |
1332 | delete [] runNumber; | |
df0a2a0a | 1333 | delete [] index; |
1334 | delete [] newBins; | |
d02b8756 | 1335 | // |
df0a2a0a | 1336 | return graphNew; |
1337 | } | |
1338 | ||
377a7d60 | 1339 | |
1340 | ||
1341 | // | |
d02b8756 | 1342 | // functions used for the trending |
377a7d60 | 1343 | // |
1344 | ||
1345 | Int_t TStatToolkit::MakeStatAlias(TTree * tree, const char * expr, const char * cut, const char * alias) | |
1346 | { | |
1347 | // | |
1348 | // Add alias using statistical values of a given variable. | |
1349 | // (by MI, Patrick Reichelt) | |
1350 | // | |
1351 | // tree - input tree | |
1352 | // expr - variable expression | |
1353 | // cut - selection criteria | |
1354 | // Output - return number of entries used to define variable | |
1355 | // In addition mean, rms, median, and robust mean and rms (choosing fraction of data with smallest RMS) | |
1356 | // | |
d02b8756 | 1357 | /* Example usage: |
1358 | 1.) create the robust estimators for variable expr="QA.TPC.CPass1.meanTPCncl" and create a corresponding | |
1359 | aliases with the prefix alias[0]="ncl", calculated using fraction alias[1]="0.90" | |
1360 | ||
1361 | TStatToolkit::MakeStatAlias(tree,"QA.TPC.CPass1.meanTPCncl","QA.TPC.CPass1.status>0","ncl:0.9"); | |
1362 | root [4] tree->GetListOfAliases().Print() | |
1363 | OBJ: TNamed ncl_Median (130.964333+0) | |
1364 | OBJ: TNamed ncl_Mean (122.120387+0) | |
1365 | OBJ: TNamed ncl_RMS (33.509623+0) | |
1366 | OBJ: TNamed ncl_Mean90 (131.503862+0) | |
1367 | OBJ: TNamed ncl_RMS90 (3.738260+0) | |
1368 | */ | |
377a7d60 | 1369 | // |
d02b8756 | 1370 | Int_t entries = tree->Draw(expr,cut,"goff"); |
377a7d60 | 1371 | if (entries<=1){ |
1372 | printf("Expression or cut not valid:\t%s\t%s\n", expr, cut); | |
1373 | return 0; | |
1374 | } | |
1375 | // | |
1376 | TObjArray* oaAlias = TString(alias).Tokenize(":"); | |
1377 | if (oaAlias->GetEntries()<2) return 0; | |
1378 | Float_t entryFraction = atof( oaAlias->At(1)->GetName() ); | |
1379 | // | |
1380 | Double_t median = TMath::Median(entries,tree->GetV1()); | |
d02b8756 | 1381 | Double_t mean = TMath::Mean(entries,tree->GetV1()); |
377a7d60 | 1382 | Double_t rms = TMath::RMS(entries,tree->GetV1()); |
1383 | Double_t meanEF=0, rmsEF=0; | |
1384 | TStatToolkit::EvaluateUni(entries, tree->GetV1(), meanEF, rmsEF, entries*entryFraction); | |
1385 | // | |
d02b8756 | 1386 | tree->SetAlias(Form("%s_Median",oaAlias->At(0)->GetName()), Form("(%f+0)",median)); |
377a7d60 | 1387 | tree->SetAlias(Form("%s_Mean",oaAlias->At(0)->GetName()), Form("(%f+0)",mean)); |
1388 | tree->SetAlias(Form("%s_RMS",oaAlias->At(0)->GetName()), Form("(%f+0)",rms)); | |
377a7d60 | 1389 | tree->SetAlias(Form("%s_Mean%d",oaAlias->At(0)->GetName(),Int_t(entryFraction*100)), Form("(%f+0)",meanEF)); |
1390 | tree->SetAlias(Form("%s_RMS%d",oaAlias->At(0)->GetName(),Int_t(entryFraction*100)), Form("(%f+0)",rmsEF)); | |
1391 | delete oaAlias; | |
1392 | return entries; | |
1393 | } | |
1394 | ||
1395 | Int_t TStatToolkit::SetStatusAlias(TTree * tree, const char * expr, const char * cut, const char * alias) | |
1396 | { | |
1397 | // | |
1398 | // Add alias to trending tree using statistical values of a given variable. | |
1399 | // (by MI, Patrick Reichelt) | |
1400 | // | |
1401 | // format of expr : varname (e.g. meanTPCncl) | |
1402 | // format of cut : char like in TCut | |
1403 | // format of alias: alias:query:entryFraction(EF) (fraction of entries used for uniformity evaluation) | |
1404 | // e.g.: varname_Out:(abs(varname-meanEF)>6.*rmsEF):0.8 | |
d02b8756 | 1405 | // available internal variables are: 'varname, Median, Mean, MeanEF, RMS, RMSEF' |
1406 | // in the alias, 'varname' will be replaced by its content, and 'EF' by the percentage (e.g. MeanEF -> Mean80) | |
377a7d60 | 1407 | // |
1408 | /* Example usage: | |
d02b8756 | 1409 | 1.) Define robust mean (possible, but easier done with TStatToolkit::MakeStatAlias(...)) |
1410 | TStatToolkit::SetStatusAlias(tree, "meanTPCnclF", "meanTPCnclF>0", "meanTPCnclF_MeanEF:MeanEF:0.80") ; | |
377a7d60 | 1411 | root [10] tree->GetListOfAliases()->Print() |
1412 | Collection name='TList', class='TList', size=1 | |
d02b8756 | 1413 | OBJ: TNamed meanTPCnclF_Mean80 0.899308 |
377a7d60 | 1414 | 2.) create alias outlyers - 6 sigma cut |
d02b8756 | 1415 | TStatToolkit::SetStatusAlias(tree, "meanTPCnclF", "meanTPCnclF>0", "meanTPCnclF_Out:(abs(meanTPCnclF-MeanEF)>6.*RMSEF):0.8") |
1416 | meanTPCnclF_Out ==> (abs(meanTPCnclF-0.899308)>6.*0.016590) | |
377a7d60 | 1417 | 3.) the same functionality as in 2.) |
d02b8756 | 1418 | TStatToolkit::SetStatusAlias(tree, "meanTPCnclF", "meanTPCnclF>0", "varname_Out2:(abs(varname-MeanEF)>6.*RMSEF):0.8") |
1419 | meanTPCnclF_Out2 ==> (abs(meanTPCnclF-0.899308)>6.*0.016590) | |
377a7d60 | 1420 | */ |
1421 | // | |
d02b8756 | 1422 | Int_t entries = tree->Draw(expr,cut,"goff"); |
1423 | if (entries<=1){ | |
1424 | printf("Expression or cut not valid:\t%s\t%s\n", expr, cut); | |
1425 | return 0; | |
1426 | } | |
1427 | // | |
377a7d60 | 1428 | TObjArray* oaVar = TString(expr).Tokenize(":"); |
1429 | char varname[50]; | |
377a7d60 | 1430 | snprintf(varname,50,"%s", oaVar->At(0)->GetName()); |
d02b8756 | 1431 | // |
377a7d60 | 1432 | TObjArray* oaAlias = TString(alias).Tokenize(":"); |
d02b8756 | 1433 | if (oaAlias->GetEntries()<3) return 0; |
377a7d60 | 1434 | Float_t entryFraction = atof( oaAlias->At(2)->GetName() ); |
1435 | // | |
377a7d60 | 1436 | Double_t median = TMath::Median(entries,tree->GetV1()); |
d02b8756 | 1437 | Double_t mean = TMath::Mean(entries,tree->GetV1()); |
377a7d60 | 1438 | Double_t rms = TMath::RMS(entries,tree->GetV1()); |
1439 | Double_t meanEF=0, rmsEF=0; | |
1440 | TStatToolkit::EvaluateUni(entries, tree->GetV1(), meanEF, rmsEF, entries*entryFraction); | |
377a7d60 | 1441 | // |
1442 | TString sAlias( oaAlias->At(0)->GetName() ); | |
1443 | sAlias.ReplaceAll("varname",varname); | |
d02b8756 | 1444 | sAlias.ReplaceAll("MeanEF", Form("Mean%1.0f",entryFraction*100) ); |
1445 | sAlias.ReplaceAll("RMSEF", Form("RMS%1.0f",entryFraction*100) ); | |
377a7d60 | 1446 | TString sQuery( oaAlias->At(1)->GetName() ); |
1447 | sQuery.ReplaceAll("varname",varname); | |
1448 | sQuery.ReplaceAll("MeanEF", Form("%f",meanEF) ); | |
d02b8756 | 1449 | sQuery.ReplaceAll("RMSEF", Form("%f",rmsEF) ); //make sure to replace 'RMSEF' before 'RMS'... |
377a7d60 | 1450 | sQuery.ReplaceAll("Median", Form("%f",median) ); |
d02b8756 | 1451 | sQuery.ReplaceAll("Mean", Form("%f",mean) ); |
377a7d60 | 1452 | sQuery.ReplaceAll("RMS", Form("%f",rms) ); |
d02b8756 | 1453 | printf("define alias:\t%s = %s\n", sAlias.Data(), sQuery.Data()); |
377a7d60 | 1454 | // |
1455 | char query[200]; | |
1456 | char aname[200]; | |
1457 | snprintf(query,200,"%s", sQuery.Data()); | |
1458 | snprintf(aname,200,"%s", sAlias.Data()); | |
1459 | tree->SetAlias(aname, query); | |
d02b8756 | 1460 | delete oaVar; |
1461 | delete oaAlias; | |
377a7d60 | 1462 | return entries; |
1463 | } | |
1464 | ||
1465 | TMultiGraph* TStatToolkit::MakeStatusMultGr(TTree * tree, const char * expr, const char * cut, const char * alias, Int_t igr) | |
1466 | { | |
1467 | // | |
1468 | // Compute a trending multigraph that shows for which runs a variable has outliers. | |
1469 | // (by MI, Patrick Reichelt) | |
1470 | // | |
1471 | // format of expr : varname:xaxis (e.g. meanTPCncl:run) | |
1472 | // format of cut : char like in TCut | |
1473 | // format of alias: (1):(varname_Out==0):(varname_Out)[:(varname_Warning):...] | |
1474 | // in the alias, 'varname' will be replaced by its content (e.g. varname_Out -> meanTPCncl_Out) | |
d02b8756 | 1475 | // note: the aliases 'varname_Out' etc have to be defined by function TStatToolkit::SetStatusAlias(...) |
377a7d60 | 1476 | // counter igr is used to shift the multigraph in y when filling a TObjArray. |
1477 | // | |
1478 | TObjArray* oaVar = TString(expr).Tokenize(":"); | |
d02b8756 | 1479 | if (oaVar->GetEntries()<2) return 0; |
377a7d60 | 1480 | char varname[50]; |
1481 | char var_x[50]; | |
1482 | snprintf(varname,50,"%s", oaVar->At(0)->GetName()); | |
1483 | snprintf(var_x ,50,"%s", oaVar->At(1)->GetName()); | |
d02b8756 | 1484 | // |
377a7d60 | 1485 | TString sAlias(alias); |
1486 | sAlias.ReplaceAll("varname",varname); | |
1487 | TObjArray* oaAlias = TString(sAlias.Data()).Tokenize(":"); | |
d02b8756 | 1488 | if (oaAlias->GetEntries()<3) return 0; |
377a7d60 | 1489 | // |
1490 | char query[200]; | |
1491 | TMultiGraph* multGr = new TMultiGraph(); | |
1492 | Int_t marArr[6] = {24+igr%2, 20+igr%2, 20+igr%2, 20+igr%2, 22, 23}; | |
1493 | Int_t colArr[6] = {kBlack, kBlack, kRed, kOrange, kMagenta, kViolet}; | |
1494 | Double_t sizArr[6] = {1.2, 1.1, 1.0, 1.0, 1, 1}; | |
1495 | const Int_t ngr = oaAlias->GetEntriesFast(); | |
1496 | for (Int_t i=0; i<ngr; i++){ | |
1497 | if (i==2) continue; // the Fatal(Out) graph will be added in the end to be plotted on top! | |
1498 | snprintf(query,200, "%f*(%s-0.5):%s", 1.+igr, oaAlias->At(i)->GetName(), var_x); | |
d02b8756 | 1499 | multGr->Add( (TGraphErrors*) TStatToolkit::MakeGraphSparse(tree,query,cut,marArr[i],colArr[i],sizArr[i]) ); |
377a7d60 | 1500 | } |
1501 | snprintf(query,200, "%f*(%s-0.5):%s", 1.+igr, oaAlias->At(2)->GetName(), var_x); | |
d02b8756 | 1502 | multGr->Add( (TGraphErrors*) TStatToolkit::MakeGraphSparse(tree,query,cut,marArr[2],colArr[2],sizArr[2]) ); |
377a7d60 | 1503 | // |
1504 | multGr->SetName(varname); | |
1505 | multGr->SetTitle(varname); // used for y-axis labels. // details to be included! | |
d02b8756 | 1506 | delete oaVar; |
1507 | delete oaAlias; | |
377a7d60 | 1508 | return multGr; |
1509 | } | |
1510 | ||
1511 | ||
1512 | void TStatToolkit::AddStatusPad(TCanvas* c1, Float_t padratio, Float_t bottommargin) | |
1513 | { | |
1514 | // | |
1515 | // add pad to bottom of canvas for Status graphs (by Patrick Reichelt) | |
1516 | // call function "DrawStatusGraphs(...)" afterwards | |
1517 | // | |
1518 | TCanvas* c1_clone = (TCanvas*) c1->Clone("c1_clone"); | |
1519 | c1->Clear(); | |
1520 | // produce new pads | |
1521 | c1->cd(); | |
1522 | TPad* pad1 = new TPad("pad1", "pad1", 0., padratio, 1., 1.); | |
1523 | pad1->Draw(); | |
1524 | pad1->SetNumber(1); // so it can be called via "c1->cd(1);" | |
1525 | c1->cd(); | |
1526 | TPad* pad2 = new TPad("pad2", "pad2", 0., 0., 1., padratio); | |
1527 | pad2->Draw(); | |
1528 | pad2->SetNumber(2); | |
1529 | // draw original canvas into first pad | |
1530 | c1->cd(1); | |
1531 | c1_clone->DrawClonePad(); | |
1532 | pad1->SetBottomMargin(0.001); | |
1533 | pad1->SetRightMargin(0.01); | |
1534 | // set up second pad | |
1535 | c1->cd(2); | |
1536 | pad2->SetGrid(3); | |
1537 | pad2->SetTopMargin(0); | |
1538 | pad2->SetBottomMargin(bottommargin); // for the long x-axis labels (runnumbers) | |
1539 | pad2->SetRightMargin(0.01); | |
1540 | } | |
1541 | ||
1542 | ||
1543 | void TStatToolkit::DrawStatusGraphs(TObjArray* oaMultGr) | |
1544 | { | |
1545 | // | |
1546 | // draw Status graphs into active pad of canvas (by MI, Patrick Reichelt) | |
1547 | // ...into bottom pad, if called after "AddStatusPad(...)" | |
1548 | // | |
1549 | const Int_t nvars = oaMultGr->GetEntriesFast(); | |
1550 | TGraph* grAxis = (TGraph*) ((TMultiGraph*) oaMultGr->At(0))->GetListOfGraphs()->At(0); | |
1551 | grAxis->SetMaximum(0.5*nvars+0.5); | |
1552 | grAxis->SetMinimum(0); | |
1553 | grAxis->GetYaxis()->SetLabelSize(0); | |
1554 | Int_t entries = grAxis->GetN(); | |
1555 | printf("entries (via GetN()) = %d\n",entries); | |
1556 | grAxis->GetXaxis()->SetLabelSize(5.7*TMath::Min(TMath::Max(5./entries,0.01),0.03)); | |
1557 | grAxis->GetXaxis()->LabelsOption("v"); | |
1558 | grAxis->Draw("ap"); | |
1559 | // | |
1560 | // draw multigraphs & names of status variables on the y axis | |
1561 | for (Int_t i=0; i<nvars; i++){ | |
1562 | ((TMultiGraph*) oaMultGr->At(i))->Draw("p"); | |
1563 | TLatex* ylabel = new TLatex(-0.1, 0.5*i+0.5, ((TMultiGraph*) oaMultGr->At(i))->GetTitle()); | |
1564 | ylabel->SetTextAlign(32); //hor:right & vert:centered | |
1565 | ylabel->SetTextSize(0.025/gPad->GetHNDC()); | |
1566 | ylabel->Draw(); | |
1567 | } | |
1568 | } | |
376089b6 | 1569 | |
1570 | ||
1571 | void TStatToolkit::DrawHistogram(TTree * tree, const char* drawCommand, const char* cuts, const char* histoname, const char* histotitle, Int_t nsigma, Float_t fraction ) | |
1572 | { | |
1573 | // | |
1574 | // Draw histogram from TTree with robust range | |
1575 | // Only for 1D so far! | |
1576 | // | |
1577 | // Parameters: | |
1578 | // - histoname: name of histogram | |
1579 | // - histotitle: title of histgram | |
1580 | // - fraction: fraction of data to define the robust mean | |
1581 | // - nsigma: nsigma value for range | |
1582 | // | |
1583 | ||
1584 | TString drawStr(drawCommand); | |
1585 | TString cutStr(cuts); | |
1586 | Int_t dim = 1; | |
1587 | ||
376089b6 | 1588 | if(!tree) { |
1589 | cerr<<" Tree pointer is NULL!"<<endl; | |
1590 | return; | |
1591 | } | |
1592 | ||
3240a856 | 1593 | // get entries |
376089b6 | 1594 | Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff"); |
1595 | if (entries == -1) { | |
1596 | cerr<<"TTree draw returns -1"<<endl; | |
1597 | return; | |
1598 | } | |
1599 | ||
3240a856 | 1600 | // get dimension |
1601 | if(tree->GetV1()) dim = 1; | |
1602 | if(tree->GetV2()) dim = 2; | |
1603 | if(tree->GetV3()) dim = 3; | |
1604 | if(dim > 2){ | |
1605 | cerr<<"TTree has more than 2 dimensions (not yet supported)"<<endl; | |
1606 | return; | |
1607 | } | |
1608 | ||
1609 | // draw robust | |
1610 | Double_t meanX, rmsX=0; | |
1611 | Double_t meanY, rmsY=0; | |
1612 | TStatToolkit::EvaluateUni(entries, tree->GetV1(),meanX,rmsX, fraction*entries); | |
1613 | if(dim==2){ | |
1614 | TStatToolkit::EvaluateUni(entries, tree->GetV1(),meanY,rmsY, fraction*entries); | |
1615 | TStatToolkit::EvaluateUni(entries, tree->GetV2(),meanX,rmsX, fraction*entries); | |
1616 | } | |
1617 | TH1* hOut; | |
1618 | if(dim==1){ | |
1619 | hOut = new TH1F(histoname, histotitle, 200, meanX-nsigma*rmsX, meanX+nsigma*rmsX); | |
1620 | for (Int_t i=0; i<entries; i++) hOut->Fill(tree->GetV1()[i]); | |
1621 | hOut->GetXaxis()->SetTitle(tree->GetHistogram()->GetXaxis()->GetTitle()); | |
1622 | hOut->Draw(); | |
1623 | } | |
1624 | else if(dim==2){ | |
1625 | hOut = new TH2F(histoname, histotitle, 200, meanX-nsigma*rmsX, meanX+nsigma*rmsX,200, meanY-nsigma*rmsY, meanY+nsigma*rmsY); | |
1626 | for (Int_t i=0; i<entries; i++) hOut->Fill(tree->GetV2()[i],tree->GetV1()[i]); | |
1627 | hOut->GetXaxis()->SetTitle(tree->GetHistogram()->GetXaxis()->GetTitle()); | |
1628 | hOut->GetYaxis()->SetTitle(tree->GetHistogram()->GetYaxis()->GetTitle()); | |
1629 | hOut->Draw("colz"); | |
1630 | } | |
376089b6 | 1631 | |
1632 | } |