]>
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 | // | |
8fae6121 MK |
22 | // |
23 | ///////////////////////////////////////////////////////////////////////// | |
21f3a443 | 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 | // includes neccessary for test functions | |
43 | // | |
44 | #include "TSystem.h" | |
45 | #include "TRandom.h" | |
46 | #include "TStopwatch.h" | |
47 | #include "TTreeStream.h" | |
48 | ||
49 | #include "TStatToolkit.h" | |
50 | ||
e3ad0e02 | 51 | |
52 | using std::cout; | |
53 | using std::cerr; | |
54 | using std::endl; | |
55 | ||
21f3a443 | 56 | |
57 | ClassImp(TStatToolkit) // Class implementation to enable ROOT I/O | |
58 | ||
59 | TStatToolkit::TStatToolkit() : TObject() | |
60 | { | |
61 | // | |
62 | // Default constructor | |
63 | // | |
64 | } | |
65 | /////////////////////////////////////////////////////////////////////////// | |
66 | TStatToolkit::~TStatToolkit() | |
67 | { | |
68 | // | |
69 | // Destructor | |
70 | // | |
71 | } | |
72 | ||
73 | ||
74 | //_____________________________________________________________________________ | |
75 | void TStatToolkit::EvaluateUni(Int_t nvectors, Double_t *data, Double_t &mean | |
76 | , Double_t &sigma, Int_t hh) | |
77 | { | |
78 | // | |
79 | // Robust estimator in 1D case MI version - (faster than ROOT version) | |
80 | // | |
81 | // For the univariate case | |
82 | // estimates of location and scatter are returned in mean and sigma parameters | |
83 | // the algorithm works on the same principle as in multivariate case - | |
84 | // it finds a subset of size hh with smallest sigma, and then returns mean and | |
85 | // sigma of this subset | |
86 | // | |
87 | ||
88 | if (hh==0) | |
89 | hh=(nvectors+2)/2; | |
90 | 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}; | |
91 | Int_t *index=new Int_t[nvectors]; | |
92 | TMath::Sort(nvectors, data, index, kFALSE); | |
93 | ||
94 | Int_t nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11); | |
95 | Double_t factor = faclts[TMath::Max(0,nquant-1)]; | |
96 | ||
97 | Double_t sumx =0; | |
98 | Double_t sumx2 =0; | |
99 | Int_t bestindex = -1; | |
100 | Double_t bestmean = 0; | |
101 | Double_t bestsigma = (data[index[nvectors-1]]-data[index[0]]+1.); // maximal possible sigma | |
102 | bestsigma *=bestsigma; | |
103 | ||
104 | for (Int_t i=0; i<hh; i++){ | |
105 | sumx += data[index[i]]; | |
106 | sumx2 += data[index[i]]*data[index[i]]; | |
107 | } | |
108 | ||
109 | Double_t norm = 1./Double_t(hh); | |
bd7b4d18 | 110 | Double_t norm2 = (hh-1)>0 ? 1./Double_t(hh-1):1; |
21f3a443 | 111 | for (Int_t i=hh; i<nvectors; i++){ |
112 | Double_t cmean = sumx*norm; | |
113 | Double_t csigma = (sumx2 - hh*cmean*cmean)*norm2; | |
114 | if (csigma<bestsigma){ | |
115 | bestmean = cmean; | |
116 | bestsigma = csigma; | |
117 | bestindex = i-hh; | |
118 | } | |
119 | ||
120 | sumx += data[index[i]]-data[index[i-hh]]; | |
121 | sumx2 += data[index[i]]*data[index[i]]-data[index[i-hh]]*data[index[i-hh]]; | |
122 | } | |
123 | ||
124 | Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma)); | |
125 | mean = bestmean; | |
126 | sigma = bstd; | |
127 | delete [] index; | |
128 | ||
129 | } | |
130 | ||
131 | ||
132 | ||
133 | void TStatToolkit::EvaluateUniExternal(Int_t nvectors, Double_t *data, Double_t &mean, Double_t &sigma, Int_t hh, Float_t externalfactor) | |
134 | { | |
135 | // Modified version of ROOT robust EvaluateUni | |
136 | // robust estimator in 1D case MI version | |
137 | // added external factor to include precision of external measurement | |
138 | // | |
139 | ||
140 | if (hh==0) | |
141 | hh=(nvectors+2)/2; | |
142 | 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}; | |
143 | Int_t *index=new Int_t[nvectors]; | |
144 | TMath::Sort(nvectors, data, index, kFALSE); | |
145 | // | |
146 | Int_t nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11); | |
147 | Double_t factor = faclts[0]; | |
148 | if (nquant>0){ | |
149 | // fix proper normalization - Anja | |
150 | factor = faclts[nquant-1]; | |
151 | } | |
152 | ||
153 | // | |
154 | // | |
155 | Double_t sumx =0; | |
156 | Double_t sumx2 =0; | |
157 | Int_t bestindex = -1; | |
158 | Double_t bestmean = 0; | |
159 | Double_t bestsigma = -1; | |
160 | for (Int_t i=0; i<hh; i++){ | |
161 | sumx += data[index[i]]; | |
162 | sumx2 += data[index[i]]*data[index[i]]; | |
163 | } | |
164 | // | |
165 | Double_t kfactor = 2.*externalfactor - externalfactor*externalfactor; | |
166 | Double_t norm = 1./Double_t(hh); | |
167 | for (Int_t i=hh; i<nvectors; i++){ | |
168 | Double_t cmean = sumx*norm; | |
169 | Double_t csigma = (sumx2*norm - cmean*cmean*kfactor); | |
170 | if (csigma<bestsigma || bestsigma<0){ | |
171 | bestmean = cmean; | |
172 | bestsigma = csigma; | |
173 | bestindex = i-hh; | |
174 | } | |
175 | // | |
176 | // | |
177 | sumx += data[index[i]]-data[index[i-hh]]; | |
178 | sumx2 += data[index[i]]*data[index[i]]-data[index[i-hh]]*data[index[i-hh]]; | |
179 | } | |
180 | ||
181 | Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma)); | |
182 | mean = bestmean; | |
183 | sigma = bstd; | |
184 | delete [] index; | |
185 | } | |
186 | ||
187 | ||
188 | //_____________________________________________________________________________ | |
189 | Int_t TStatToolkit::Freq(Int_t n, const Int_t *inlist | |
190 | , Int_t *outlist, Bool_t down) | |
191 | { | |
192 | // | |
193 | // Sort eleements according occurancy | |
194 | // The size of output array has is 2*n | |
195 | // | |
196 | ||
197 | Int_t * sindexS = new Int_t[n]; // temp array for sorting | |
198 | Int_t * sindexF = new Int_t[2*n]; | |
b8072cce | 199 | for (Int_t i=0;i<n;i++) sindexS[i]=0; |
200 | for (Int_t i=0;i<2*n;i++) sindexF[i]=0; | |
21f3a443 | 201 | // |
202 | TMath::Sort(n,inlist, sindexS, down); | |
203 | Int_t last = inlist[sindexS[0]]; | |
204 | Int_t val = last; | |
205 | sindexF[0] = 1; | |
206 | sindexF[0+n] = last; | |
207 | Int_t countPos = 0; | |
208 | // | |
209 | // find frequency | |
210 | for(Int_t i=1;i<n; i++){ | |
211 | val = inlist[sindexS[i]]; | |
212 | if (last == val) sindexF[countPos]++; | |
213 | else{ | |
214 | countPos++; | |
215 | sindexF[countPos+n] = val; | |
216 | sindexF[countPos]++; | |
217 | last =val; | |
218 | } | |
219 | } | |
220 | if (last==val) countPos++; | |
221 | // sort according frequency | |
222 | TMath::Sort(countPos, sindexF, sindexS, kTRUE); | |
223 | for (Int_t i=0;i<countPos;i++){ | |
224 | outlist[2*i ] = sindexF[sindexS[i]+n]; | |
225 | outlist[2*i+1] = sindexF[sindexS[i]]; | |
226 | } | |
227 | delete [] sindexS; | |
228 | delete [] sindexF; | |
229 | ||
230 | return countPos; | |
231 | ||
232 | } | |
233 | ||
234 | //___TStatToolkit__________________________________________________________________________ | |
3d7cc0b4 | 235 | void TStatToolkit::TruncatedMean(const TH1 * his, TVectorD *param, Float_t down, Float_t up, Bool_t verbose){ |
21f3a443 | 236 | // |
237 | // | |
238 | // | |
239 | Int_t nbins = his->GetNbinsX(); | |
240 | Float_t nentries = his->GetEntries(); | |
241 | Float_t sum =0; | |
242 | Float_t mean = 0; | |
243 | Float_t sigma2 = 0; | |
244 | Float_t ncumul=0; | |
245 | for (Int_t ibin=1;ibin<nbins; ibin++){ | |
246 | ncumul+= his->GetBinContent(ibin); | |
247 | Float_t fraction = Float_t(ncumul)/Float_t(nentries); | |
248 | if (fraction>down && fraction<up){ | |
249 | sum+=his->GetBinContent(ibin); | |
250 | mean+=his->GetBinCenter(ibin)*his->GetBinContent(ibin); | |
251 | sigma2+=his->GetBinCenter(ibin)*his->GetBinCenter(ibin)*his->GetBinContent(ibin); | |
252 | } | |
253 | } | |
254 | mean/=sum; | |
255 | sigma2= TMath::Sqrt(TMath::Abs(sigma2/sum-mean*mean)); | |
256 | if (param){ | |
257 | (*param)[0] = his->GetMaximum(); | |
258 | (*param)[1] = mean; | |
259 | (*param)[2] = sigma2; | |
260 | ||
261 | } | |
262 | if (verbose) printf("Mean\t%f\t Sigma2\t%f\n", mean,sigma2); | |
263 | } | |
264 | ||
32d8f622 | 265 | void TStatToolkit::LTM(TH1 * his, TVectorD *param , Float_t fraction, Bool_t verbose){ |
21f3a443 | 266 | // |
32d8f622 | 267 | // LTM : Trimmed mean on histogram - Modified version for binned data |
21f3a443 | 268 | // |
32d8f622 | 269 | // Robust statistic to estimate properties of the distribution |
270 | // See http://en.wikipedia.org/w/index.php?title=Trimmed_estimator&oldid=582847999 | |
271 | // | |
272 | // New faster version is under preparation | |
273 | // | |
274 | if (!param) return; | |
275 | (*param)[0]=0; | |
276 | (*param)[1]=0; | |
277 | (*param)[2]=0; | |
21f3a443 | 278 | Int_t nbins = his->GetNbinsX(); |
279 | Int_t nentries = (Int_t)his->GetEntries(); | |
32d8f622 | 280 | if (nentries<=0) return; |
21f3a443 | 281 | Double_t *data = new Double_t[nentries]; |
282 | Int_t npoints=0; | |
283 | for (Int_t ibin=1;ibin<nbins; ibin++){ | |
8ecd39c2 | 284 | Double_t entriesI = his->GetBinContent(ibin); |
285 | //Double_t xcenter= his->GetBinCenter(ibin); | |
286 | Double_t x0 = his->GetXaxis()->GetBinLowEdge(ibin); | |
287 | Double_t w = his->GetXaxis()->GetBinWidth(ibin); | |
21f3a443 | 288 | for (Int_t ic=0; ic<entriesI; ic++){ |
289 | if (npoints<nentries){ | |
8ecd39c2 | 290 | data[npoints]= x0+w*Double_t((ic+0.5)/entriesI); |
21f3a443 | 291 | npoints++; |
292 | } | |
293 | } | |
294 | } | |
32d8f622 | 295 | Double_t mean, sigma; |
21f3a443 | 296 | Int_t npoints2=TMath::Min(Int_t(fraction*Float_t(npoints)),npoints-1); |
297 | npoints2=TMath::Max(Int_t(0.5*Float_t(npoints)),npoints2); | |
298 | TStatToolkit::EvaluateUni(npoints, data, mean,sigma,npoints2); | |
299 | delete [] data; | |
300 | if (verbose) printf("Mean\t%f\t Sigma2\t%f\n", mean,sigma);if (param){ | |
301 | (*param)[0] = his->GetMaximum(); | |
302 | (*param)[1] = mean; | |
303 | (*param)[2] = sigma; | |
304 | } | |
305 | } | |
306 | ||
8fae6121 MK |
307 | |
308 | void TStatToolkit::MedianFilter(TH1 * his1D, Int_t nmedian){ | |
309 | // | |
310 | // Algorithm to filter histogram | |
311 | // author: marian.ivanov@cern.ch | |
312 | // Details of algorithm: | |
313 | // http://en.wikipedia.org/w/index.php?title=Median_filter&oldid=582191524 | |
314 | // Input parameters: | |
315 | // his1D - input histogam - to be modiefied by Medianfilter | |
316 | // nmendian - number of bins in median filter | |
317 | // | |
318 | Int_t nbins = his1D->GetNbinsX(); | |
319 | TVectorD vectorH(nbins); | |
320 | for (Int_t ibin=0; ibin<nbins; ibin++) vectorH[ibin]=his1D->GetBinContent(ibin+1); | |
321 | for (Int_t ibin=0; ibin<nbins; ibin++) { | |
322 | Int_t index0=ibin-nmedian; | |
323 | Int_t index1=ibin+nmedian; | |
324 | if (index0<0) {index1+=-index0; index0=0;} | |
325 | if (index1>=nbins) {index0-=index1-nbins+1; index1=nbins-1;} | |
326 | Double_t value= TMath::Median(index1-index0,&(vectorH.GetMatrixArray()[index0])); | |
327 | his1D->SetBinContent(ibin+1, value); | |
328 | } | |
329 | } | |
330 | ||
331 | Bool_t TStatToolkit::LTMHisto(TH1 *his1D, TVectorD ¶ms , Float_t fraction){ | |
332 | // | |
333 | // LTM : Trimmed mean on histogram - Modified version for binned data | |
8ecd39c2 | 334 | // |
8fae6121 | 335 | // Robust statistic to estimate properties of the distribution |
8ecd39c2 | 336 | // To handle binning error special treatment |
337 | // for definition of unbinned data see: | |
338 | // http://en.wikipedia.org/w/index.php?title=Trimmed_estimator&oldid=582847999 | |
8fae6121 | 339 | // |
8ecd39c2 | 340 | // Function parameters: |
8fae6121 MK |
341 | // his1D - input histogram |
342 | // params - vector with parameters | |
343 | // - 0 - area | |
344 | // - 1 - mean | |
8ecd39c2 | 345 | // - 2 - rms |
8fae6121 | 346 | // - 3 - error estimate of mean |
8ecd39c2 | 347 | // - 4 - error estimate of RMS |
8fae6121 MK |
348 | // - 5 - first accepted bin position |
349 | // - 6 - last accepted bin position | |
350 | // | |
351 | Int_t nbins = his1D->GetNbinsX(); | |
352 | Int_t nentries = (Int_t)his1D->GetEntries(); | |
8ecd39c2 | 353 | const Double_t kEpsilon=0.0000000001; |
354 | ||
8fae6121 | 355 | if (nentries<=0) return 0; |
8ecd39c2 | 356 | if (fraction>1) fraction=0; |
357 | if (fraction<0) return 0; | |
8fae6121 MK |
358 | TVectorD vectorX(nbins); |
359 | TVectorD vectorMean(nbins); | |
360 | TVectorD vectorRMS(nbins); | |
361 | Double_t sumCont=0; | |
8ecd39c2 | 362 | for (Int_t ibin0=1; ibin0<=nbins; ibin0++) sumCont+=his1D->GetBinContent(ibin0); |
8fae6121 MK |
363 | // |
364 | Double_t minRMS=his1D->GetRMS()*10000; | |
365 | Int_t maxBin=0; | |
366 | // | |
367 | for (Int_t ibin0=1; ibin0<nbins; ibin0++){ | |
368 | Double_t sum0=0, sum1=0, sum2=0; | |
369 | Int_t ibin1=ibin0; | |
370 | for ( ibin1=ibin0; ibin1<nbins; ibin1++){ | |
371 | Double_t cont=his1D->GetBinContent(ibin1); | |
372 | Double_t x= his1D->GetBinCenter(ibin1); | |
373 | sum0+=cont; | |
374 | sum1+=cont*x; | |
375 | sum2+=cont*x*x; | |
8ecd39c2 | 376 | if ( (ibin0!=ibin1) && sum0>=fraction*sumCont) break; |
8fae6121 MK |
377 | } |
378 | vectorX[ibin0]=his1D->GetBinCenter(ibin0); | |
8ecd39c2 | 379 | if (sum0<fraction*sumCont) continue; |
380 | // | |
381 | // substract fractions of bin0 and bin1 to keep sum0=fration*sumCont | |
382 | // | |
383 | Double_t diff = sum0-fraction*sumCont; | |
5d1391ec | 384 | Double_t mean = (sum0>0) ? sum1/sum0:0; |
8ecd39c2 | 385 | // |
386 | Double_t x0=his1D->GetBinCenter(ibin0); | |
387 | Double_t x1=his1D->GetBinCenter(ibin1); | |
388 | Double_t y0=his1D->GetBinContent(ibin0); | |
389 | Double_t y1=his1D->GetBinContent(ibin1); | |
390 | // | |
391 | Double_t d = y0+y1-diff; //enties to keep | |
392 | Double_t w0=0,w1=0; | |
393 | if (y0<=kEpsilon&&y1>kEpsilon){ | |
394 | w1=d/y1; | |
395 | } | |
396 | if (y1<=kEpsilon&&y0>kEpsilon){ | |
397 | w0=d/y0; | |
398 | } | |
399 | if (y0>kEpsilon && y1>kEpsilon && x1>x0 ){ | |
400 | w0 = (d*(x1-mean))/((x1-x0)*y0); | |
401 | w1 = (d-y0*w0)/y1; | |
402 | // | |
403 | if (w0>1) {w1+=(w0-1)*y0/y1; w0=1;} | |
404 | if (w1>1) {w0+=(w1-1)*y1/y0; w1=1;} | |
405 | } | |
406 | if ( (x1>x0) &&TMath::Abs(y0*w0+y1*w1-d)>kEpsilon*sum0){ | |
407 | printf(" TStatToolkit::LTMHisto error\n"); | |
408 | } | |
409 | sum0-=y0+y1; | |
410 | sum1-=y0*x0; | |
411 | sum1-=y1*x1; | |
412 | sum2-=y0*x0*x0; | |
413 | sum2-=y1*x1*x1; | |
414 | // | |
415 | Double_t xx0=his1D->GetXaxis()->GetBinUpEdge(ibin0)-0.5*w0*his1D->GetBinWidth(ibin0); | |
416 | Double_t xx1=his1D->GetXaxis()->GetBinLowEdge(ibin1)+0.5*w1*his1D->GetBinWidth(ibin1); | |
417 | sum0+=y0*w0+y1*w1; | |
418 | sum1+=y0*w0*xx0; | |
419 | sum1+=y1*w1*xx1; | |
420 | sum2+=y0*w0*xx0*xx0; | |
421 | sum2+=y1*w1*xx1*xx1; | |
422 | ||
423 | // | |
424 | // choose the bin with smallest rms | |
425 | // | |
426 | if (sum0>0){ | |
8fae6121 | 427 | vectorMean[ibin0]=sum1/sum0; |
8ecd39c2 | 428 | vectorRMS[ibin0]=TMath::Sqrt(TMath::Abs(sum2/sum0-vectorMean[ibin0]*vectorMean[ibin0])); |
8fae6121 MK |
429 | if (vectorRMS[ibin0]<minRMS){ |
430 | minRMS=vectorRMS[ibin0]; | |
431 | params[0]=sum0; | |
432 | params[1]=vectorMean[ibin0]; | |
433 | params[2]=vectorRMS[ibin0]; | |
434 | params[3]=vectorRMS[ibin0]/TMath::Sqrt(sumCont*fraction); | |
8ecd39c2 | 435 | params[4]=0; // what is the formula for error of RMS??? |
8fae6121 MK |
436 | params[5]=ibin0; |
437 | params[6]=ibin1; | |
438 | params[7]=his1D->GetBinCenter(ibin0); | |
439 | params[8]=his1D->GetBinCenter(ibin1); | |
440 | maxBin=ibin0; | |
441 | } | |
442 | }else{ | |
443 | break; | |
444 | } | |
445 | } | |
446 | return kTRUE; | |
8fae6121 MK |
447 | } |
448 | ||
449 | ||
94a43b22 | 450 | Double_t TStatToolkit::FitGaus(TH1* his, TVectorD *param, TMatrixD */*matrix*/, Float_t xmin, Float_t xmax, Bool_t verbose){ |
21f3a443 | 451 | // |
452 | // Fit histogram with gaussian function | |
453 | // | |
454 | // Prameters: | |
455 | // return value- chi2 - if negative ( not enough points) | |
456 | // his - input histogram | |
457 | // param - vector with parameters | |
458 | // xmin, xmax - range to fit - if xmin=xmax=0 - the full histogram range used | |
459 | // Fitting: | |
460 | // 1. Step - make logarithm | |
461 | // 2. Linear fit (parabola) - more robust - always converge | |
462 | // 3. In case of small statistic bins are averaged | |
463 | // | |
464 | static TLinearFitter fitter(3,"pol2"); | |
465 | TVectorD par(3); | |
466 | TVectorD sigma(3); | |
467 | TMatrixD mat(3,3); | |
468 | if (his->GetMaximum()<4) return -1; | |
469 | if (his->GetEntries()<12) return -1; | |
470 | if (his->GetRMS()<mat.GetTol()) return -1; | |
471 | Float_t maxEstimate = his->GetEntries()*his->GetBinWidth(1)/TMath::Sqrt((TMath::TwoPi()*his->GetRMS())); | |
472 | Int_t dsmooth = TMath::Nint(6./TMath::Sqrt(maxEstimate)); | |
473 | ||
474 | if (maxEstimate<1) return -1; | |
475 | Int_t nbins = his->GetNbinsX(); | |
476 | Int_t npoints=0; | |
477 | // | |
478 | ||
479 | ||
480 | if (xmin>=xmax){ | |
481 | xmin = his->GetXaxis()->GetXmin(); | |
482 | xmax = his->GetXaxis()->GetXmax(); | |
483 | } | |
484 | for (Int_t iter=0; iter<2; iter++){ | |
485 | fitter.ClearPoints(); | |
486 | npoints=0; | |
487 | for (Int_t ibin=1;ibin<nbins+1; ibin++){ | |
488 | Int_t countB=1; | |
489 | Float_t entriesI = his->GetBinContent(ibin); | |
490 | for (Int_t delta = -dsmooth; delta<=dsmooth; delta++){ | |
491 | if (ibin+delta>1 &&ibin+delta<nbins-1){ | |
492 | entriesI += his->GetBinContent(ibin+delta); | |
493 | countB++; | |
494 | } | |
495 | } | |
496 | entriesI/=countB; | |
497 | Double_t xcenter= his->GetBinCenter(ibin); | |
498 | if (xcenter<xmin || xcenter>xmax) continue; | |
499 | Double_t error=1./TMath::Sqrt(countB); | |
500 | Float_t cont=2; | |
501 | if (iter>0){ | |
502 | if (par[0]+par[1]*xcenter+par[2]*xcenter*xcenter>20) return 0; | |
503 | cont = TMath::Exp(par[0]+par[1]*xcenter+par[2]*xcenter*xcenter); | |
504 | if (cont>1.) error = 1./TMath::Sqrt(cont*Float_t(countB)); | |
505 | } | |
506 | if (entriesI>1&&cont>1){ | |
507 | fitter.AddPoint(&xcenter,TMath::Log(Float_t(entriesI)),error); | |
508 | npoints++; | |
509 | } | |
510 | } | |
511 | if (npoints>3){ | |
512 | fitter.Eval(); | |
513 | fitter.GetParameters(par); | |
514 | }else{ | |
515 | break; | |
516 | } | |
517 | } | |
518 | if (npoints<=3){ | |
519 | return -1; | |
520 | } | |
521 | fitter.GetParameters(par); | |
522 | fitter.GetCovarianceMatrix(mat); | |
523 | if (TMath::Abs(par[1])<mat.GetTol()) return -1; | |
524 | if (TMath::Abs(par[2])<mat.GetTol()) return -1; | |
525 | Double_t chi2 = fitter.GetChisquare()/Float_t(npoints); | |
526 | //fitter.GetParameters(); | |
527 | if (!param) param = new TVectorD(3); | |
cb1d20de | 528 | // if (!matrix) matrix = new TMatrixD(3,3); // Covariance matrix to be implemented |
21f3a443 | 529 | (*param)[1] = par[1]/(-2.*par[2]); |
530 | (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2])); | |
531 | (*param)[0] = TMath::Exp(par[0]+ par[1]* (*param)[1] + par[2]*(*param)[1]*(*param)[1]); | |
532 | if (verbose){ | |
533 | par.Print(); | |
534 | mat.Print(); | |
535 | param->Print(); | |
536 | printf("Chi2=%f\n",chi2); | |
537 | TF1 * f1= new TF1("f1","[0]*exp(-(x-[1])^2/(2*[2]*[2]))",his->GetXaxis()->GetXmin(),his->GetXaxis()->GetXmax()); | |
538 | f1->SetParameter(0, (*param)[0]); | |
539 | f1->SetParameter(1, (*param)[1]); | |
540 | f1->SetParameter(2, (*param)[2]); | |
541 | f1->Draw("same"); | |
542 | } | |
543 | return chi2; | |
544 | } | |
545 | ||
cb1d20de | 546 | Double_t TStatToolkit::FitGaus(Float_t *arr, Int_t nBins, Float_t xMin, Float_t xMax, TVectorD *param, TMatrixD */*matrix*/, Bool_t verbose){ |
21f3a443 | 547 | // |
548 | // Fit histogram with gaussian function | |
549 | // | |
550 | // Prameters: | |
551 | // nbins: size of the array and number of histogram bins | |
552 | // xMin, xMax: histogram range | |
553 | // param: paramters of the fit (0-Constant, 1-Mean, 2-Sigma) | |
554 | // matrix: covariance matrix -- not implemented yet, pass dummy matrix!!! | |
555 | // | |
556 | // Return values: | |
557 | // >0: the chi2 returned by TLinearFitter | |
558 | // -3: only three points have been used for the calculation - no fitter was used | |
559 | // -2: only two points have been used for the calculation - center of gravity was uesed for calculation | |
560 | // -1: only one point has been used for the calculation - center of gravity was uesed for calculation | |
561 | // -4: invalid result!! | |
562 | // | |
563 | // Fitting: | |
564 | // 1. Step - make logarithm | |
565 | // 2. Linear fit (parabola) - more robust - always converge | |
566 | // | |
567 | static TLinearFitter fitter(3,"pol2"); | |
568 | static TMatrixD mat(3,3); | |
569 | static Double_t kTol = mat.GetTol(); | |
570 | fitter.StoreData(kFALSE); | |
571 | fitter.ClearPoints(); | |
572 | TVectorD par(3); | |
573 | TVectorD sigma(3); | |
3d7cc0b4 | 574 | TMatrixD matA(3,3); |
21f3a443 | 575 | TMatrixD b(3,1); |
576 | Float_t rms = TMath::RMS(nBins,arr); | |
577 | Float_t max = TMath::MaxElement(nBins,arr); | |
578 | Float_t binWidth = (xMax-xMin)/(Float_t)nBins; | |
579 | ||
580 | Float_t meanCOG = 0; | |
581 | Float_t rms2COG = 0; | |
582 | Float_t sumCOG = 0; | |
583 | ||
584 | Float_t entries = 0; | |
585 | Int_t nfilled=0; | |
586 | ||
587 | for (Int_t i=0; i<nBins; i++){ | |
588 | entries+=arr[i]; | |
589 | if (arr[i]>0) nfilled++; | |
590 | } | |
591 | ||
592 | if (max<4) return -4; | |
593 | if (entries<12) return -4; | |
594 | if (rms<kTol) return -4; | |
595 | ||
596 | Int_t npoints=0; | |
597 | // | |
598 | ||
599 | // | |
600 | for (Int_t ibin=0;ibin<nBins; ibin++){ | |
601 | Float_t entriesI = arr[ibin]; | |
602 | if (entriesI>1){ | |
603 | Double_t xcenter = xMin+(ibin+0.5)*binWidth; | |
604 | ||
605 | Float_t error = 1./TMath::Sqrt(entriesI); | |
606 | Float_t val = TMath::Log(Float_t(entriesI)); | |
607 | fitter.AddPoint(&xcenter,val,error); | |
608 | if (npoints<3){ | |
3d7cc0b4 | 609 | matA(npoints,0)=1; |
610 | matA(npoints,1)=xcenter; | |
611 | matA(npoints,2)=xcenter*xcenter; | |
21f3a443 | 612 | b(npoints,0)=val; |
613 | meanCOG+=xcenter*entriesI; | |
614 | rms2COG +=xcenter*entriesI*xcenter; | |
615 | sumCOG +=entriesI; | |
616 | } | |
617 | npoints++; | |
618 | } | |
619 | } | |
620 | ||
621 | ||
622 | Double_t chi2 = 0; | |
623 | if (npoints>=3){ | |
624 | if ( npoints == 3 ){ | |
625 | //analytic calculation of the parameters for three points | |
3d7cc0b4 | 626 | matA.Invert(); |
21f3a443 | 627 | TMatrixD res(1,3); |
3d7cc0b4 | 628 | res.Mult(matA,b); |
21f3a443 | 629 | par[0]=res(0,0); |
630 | par[1]=res(0,1); | |
631 | par[2]=res(0,2); | |
632 | chi2 = -3.; | |
633 | } else { | |
634 | // use fitter for more than three points | |
635 | fitter.Eval(); | |
636 | fitter.GetParameters(par); | |
637 | fitter.GetCovarianceMatrix(mat); | |
638 | chi2 = fitter.GetChisquare()/Float_t(npoints); | |
639 | } | |
640 | if (TMath::Abs(par[1])<kTol) return -4; | |
641 | if (TMath::Abs(par[2])<kTol) return -4; | |
642 | ||
643 | if (!param) param = new TVectorD(3); | |
cb1d20de | 644 | //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 | 645 | |
646 | (*param)[1] = par[1]/(-2.*par[2]); | |
647 | (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2])); | |
648 | Double_t lnparam0 = par[0]+ par[1]* (*param)[1] + par[2]*(*param)[1]*(*param)[1]; | |
649 | if ( lnparam0>307 ) return -4; | |
650 | (*param)[0] = TMath::Exp(lnparam0); | |
651 | if (verbose){ | |
652 | par.Print(); | |
653 | mat.Print(); | |
654 | param->Print(); | |
655 | printf("Chi2=%f\n",chi2); | |
656 | TF1 * f1= new TF1("f1","[0]*exp(-(x-[1])^2/(2*[2]*[2]))",xMin,xMax); | |
657 | f1->SetParameter(0, (*param)[0]); | |
658 | f1->SetParameter(1, (*param)[1]); | |
659 | f1->SetParameter(2, (*param)[2]); | |
660 | f1->Draw("same"); | |
661 | } | |
662 | return chi2; | |
663 | } | |
664 | ||
665 | if (npoints == 2){ | |
666 | //use center of gravity for 2 points | |
667 | meanCOG/=sumCOG; | |
668 | rms2COG /=sumCOG; | |
669 | (*param)[0] = max; | |
670 | (*param)[1] = meanCOG; | |
671 | (*param)[2] = TMath::Sqrt(TMath::Abs(meanCOG*meanCOG-rms2COG)); | |
672 | chi2=-2.; | |
673 | } | |
674 | if ( npoints == 1 ){ | |
675 | meanCOG/=sumCOG; | |
676 | (*param)[0] = max; | |
677 | (*param)[1] = meanCOG; | |
678 | (*param)[2] = binWidth/TMath::Sqrt(12); | |
679 | chi2=-1.; | |
680 | } | |
681 | return chi2; | |
682 | ||
683 | } | |
684 | ||
685 | ||
3d7cc0b4 | 686 | Float_t TStatToolkit::GetCOG(const Short_t *arr, Int_t nBins, Float_t xMin, Float_t xMax, Float_t *rms, Float_t *sum) |
21f3a443 | 687 | { |
688 | // | |
689 | // calculate center of gravity rms and sum for array 'arr' with nBins an a x range xMin to xMax | |
690 | // return COG; in case of failure return xMin | |
691 | // | |
692 | Float_t meanCOG = 0; | |
693 | Float_t rms2COG = 0; | |
694 | Float_t sumCOG = 0; | |
695 | Int_t npoints = 0; | |
696 | ||
697 | Float_t binWidth = (xMax-xMin)/(Float_t)nBins; | |
698 | ||
699 | for (Int_t ibin=0; ibin<nBins; ibin++){ | |
700 | Float_t entriesI = (Float_t)arr[ibin]; | |
701 | Double_t xcenter = xMin+(ibin+0.5)*binWidth; | |
702 | if ( entriesI>0 ){ | |
703 | meanCOG += xcenter*entriesI; | |
704 | rms2COG += xcenter*entriesI*xcenter; | |
705 | sumCOG += entriesI; | |
706 | npoints++; | |
707 | } | |
708 | } | |
709 | if ( sumCOG == 0 ) return xMin; | |
710 | meanCOG/=sumCOG; | |
711 | ||
712 | if ( rms ){ | |
713 | rms2COG /=sumCOG; | |
714 | (*rms) = TMath::Sqrt(TMath::Abs(meanCOG*meanCOG-rms2COG)); | |
715 | if ( npoints == 1 ) (*rms) = binWidth/TMath::Sqrt(12); | |
716 | } | |
717 | ||
718 | if ( sum ) | |
719 | (*sum) = sumCOG; | |
720 | ||
721 | return meanCOG; | |
722 | } | |
723 | ||
724 | ||
725 | ||
726 | /////////////////////////////////////////////////////////////// | |
727 | ////////////// TEST functions ///////////////////////// | |
728 | /////////////////////////////////////////////////////////////// | |
729 | ||
730 | ||
731 | ||
732 | ||
733 | ||
734 | void TStatToolkit::TestGausFit(Int_t nhistos){ | |
735 | // | |
736 | // Test performance of the parabolic - gaussian fit - compare it with | |
737 | // ROOT gauss fit | |
738 | // nhistos - number of histograms to be used for test | |
739 | // | |
32d8f622 | 740 | TTreeSRedirector *pcstream = new TTreeSRedirector("fitdebug.root","recreate"); |
21f3a443 | 741 | |
742 | Float_t *xTrue = new Float_t[nhistos]; | |
743 | Float_t *sTrue = new Float_t[nhistos]; | |
744 | TVectorD **par1 = new TVectorD*[nhistos]; | |
745 | TVectorD **par2 = new TVectorD*[nhistos]; | |
746 | TMatrixD dummy(3,3); | |
747 | ||
748 | ||
749 | TH1F **h1f = new TH1F*[nhistos]; | |
750 | TF1 *myg = new TF1("myg","gaus"); | |
751 | TF1 *fit = new TF1("fit","gaus"); | |
752 | gRandom->SetSeed(0); | |
753 | ||
754 | //init | |
755 | for (Int_t i=0;i<nhistos; i++){ | |
756 | par1[i] = new TVectorD(3); | |
757 | par2[i] = new TVectorD(3); | |
758 | h1f[i] = new TH1F(Form("h1f%d",i),Form("h1f%d",i),20,-10,10); | |
759 | xTrue[i]= gRandom->Rndm(); | |
760 | gSystem->Sleep(2); | |
761 | sTrue[i]= .75+gRandom->Rndm()*.5; | |
762 | myg->SetParameters(1,xTrue[i],sTrue[i]); | |
763 | h1f[i]->FillRandom("myg"); | |
764 | } | |
765 | ||
32d8f622 | 766 | TStopwatch s; |
21f3a443 | 767 | s.Start(); |
768 | //standard gaus fit | |
769 | for (Int_t i=0; i<nhistos; i++){ | |
770 | h1f[i]->Fit(fit,"0q"); | |
771 | (*par1[i])(0) = fit->GetParameter(0); | |
772 | (*par1[i])(1) = fit->GetParameter(1); | |
773 | (*par1[i])(2) = fit->GetParameter(2); | |
774 | } | |
775 | s.Stop(); | |
776 | printf("Gaussian fit\t"); | |
777 | s.Print(); | |
778 | ||
779 | s.Start(); | |
780 | //TStatToolkit gaus fit | |
781 | for (Int_t i=0; i<nhistos; i++){ | |
782 | TStatToolkit::FitGaus(h1f[i]->GetArray()+1,h1f[i]->GetNbinsX(),h1f[i]->GetXaxis()->GetXmin(),h1f[i]->GetXaxis()->GetXmax(),par2[i],&dummy); | |
783 | } | |
784 | ||
785 | s.Stop(); | |
786 | printf("Parabolic fit\t"); | |
787 | s.Print(); | |
788 | //write stream | |
789 | for (Int_t i=0;i<nhistos; i++){ | |
790 | Float_t xt = xTrue[i]; | |
791 | Float_t st = sTrue[i]; | |
792 | (*pcstream)<<"data" | |
793 | <<"xTrue="<<xt | |
794 | <<"sTrue="<<st | |
795 | <<"pg.="<<(par1[i]) | |
796 | <<"pa.="<<(par2[i]) | |
797 | <<"\n"; | |
798 | } | |
799 | //delete pointers | |
800 | for (Int_t i=0;i<nhistos; i++){ | |
801 | delete par1[i]; | |
802 | delete par2[i]; | |
803 | delete h1f[i]; | |
804 | } | |
805 | delete pcstream; | |
806 | delete []h1f; | |
807 | delete []xTrue; | |
808 | delete []sTrue; | |
809 | // | |
810 | delete []par1; | |
811 | delete []par2; | |
812 | ||
813 | } | |
814 | ||
815 | ||
816 | ||
817 | TGraph2D * TStatToolkit::MakeStat2D(TH3 * his, Int_t delta0, Int_t delta1, Int_t type){ | |
818 | // | |
819 | // | |
820 | // | |
821 | // delta - number of bins to integrate | |
822 | // type - 0 - mean value | |
823 | ||
824 | TAxis * xaxis = his->GetXaxis(); | |
825 | TAxis * yaxis = his->GetYaxis(); | |
826 | // TAxis * zaxis = his->GetZaxis(); | |
827 | Int_t nbinx = xaxis->GetNbins(); | |
828 | Int_t nbiny = yaxis->GetNbins(); | |
829 | char name[1000]; | |
830 | Int_t icount=0; | |
831 | TGraph2D *graph = new TGraph2D(nbinx*nbiny); | |
832 | TF1 f1("f1","gaus"); | |
833 | for (Int_t ix=0; ix<nbinx;ix++) | |
834 | for (Int_t iy=0; iy<nbiny;iy++){ | |
835 | Float_t xcenter = xaxis->GetBinCenter(ix); | |
836 | Float_t ycenter = yaxis->GetBinCenter(iy); | |
cb1d20de | 837 | snprintf(name,1000,"%s_%d_%d",his->GetName(), ix,iy); |
21f3a443 | 838 | TH1 *projection = his->ProjectionZ(name,ix-delta0,ix+delta0,iy-delta1,iy+delta1); |
839 | Float_t stat= 0; | |
840 | if (type==0) stat = projection->GetMean(); | |
841 | if (type==1) stat = projection->GetRMS(); | |
842 | if (type==2 || type==3){ | |
8fae6121 | 843 | TVectorD vec(10); |
21f3a443 | 844 | TStatToolkit::LTM((TH1F*)projection,&vec,0.7); |
845 | if (type==2) stat= vec[1]; | |
846 | if (type==3) stat= vec[0]; | |
847 | } | |
848 | if (type==4|| type==5){ | |
849 | projection->Fit(&f1); | |
850 | if (type==4) stat= f1.GetParameter(1); | |
851 | if (type==5) stat= f1.GetParameter(2); | |
852 | } | |
853 | //printf("%d\t%f\t%f\t%f\n", icount,xcenter, ycenter, stat); | |
854 | graph->SetPoint(icount,xcenter, ycenter, stat); | |
855 | icount++; | |
856 | } | |
857 | return graph; | |
858 | } | |
859 | ||
32d8f622 | 860 | TGraphErrors * TStatToolkit::MakeStat1D(TH2 * his, Int_t deltaBin, Double_t fraction, Int_t returnType, Int_t markerStyle, Int_t markerColor){ |
861 | // | |
862 | // function to retrieve the "mean and RMS estimate" of 2D histograms | |
863 | // | |
864 | // Robust statistic to estimate properties of the distribution | |
865 | // See http://en.wikipedia.org/wiki/Trimmed_estimator | |
866 | // | |
867 | // deltaBin - number of bins to integrate (bin+-deltaBin) | |
868 | // fraction - fraction of values for the LTM and for the gauss fit | |
869 | // returnType - | |
870 | // 0 - mean value | |
871 | // 1 - RMS | |
872 | // 2 - LTM mean | |
873 | // 3 - LTM sigma | |
874 | // 4 - Gaus fit mean - on LTM range | |
875 | // 5 - Gaus fit sigma - on LTM range | |
876 | // | |
21f3a443 | 877 | TAxis * xaxis = his->GetXaxis(); |
21f3a443 | 878 | Int_t nbinx = xaxis->GetNbins(); |
21f3a443 | 879 | char name[1000]; |
880 | Int_t icount=0; | |
32d8f622 | 881 | // |
882 | TVectorD vecX(nbinx); | |
883 | TVectorD vecXErr(nbinx); | |
884 | TVectorD vecY(nbinx); | |
885 | TVectorD vecYErr(nbinx); | |
886 | // | |
21f3a443 | 887 | TF1 f1("f1","gaus"); |
8fae6121 | 888 | TVectorD vecLTM(10); |
32d8f622 | 889 | |
8fae6121 MK |
890 | for (Int_t jx=1; jx<=nbinx;jx++){ |
891 | Int_t ix=jx-1; | |
892 | Float_t xcenter = xaxis->GetBinCenter(jx); | |
cb1d20de | 893 | snprintf(name,1000,"%s_%d",his->GetName(), ix); |
8fae6121 | 894 | TH1 *projection = his->ProjectionY(name,TMath::Max(jx-deltaBin,1),TMath::Min(jx+deltaBin,nbinx)); |
32d8f622 | 895 | Double_t stat= 0; |
896 | Double_t err =0; | |
8fae6121 | 897 | TStatToolkit::LTMHisto((TH1F*)projection,vecLTM,fraction); |
32d8f622 | 898 | // |
899 | if (returnType==0) { | |
900 | stat = projection->GetMean(); | |
901 | err = projection->GetMeanError(); | |
21f3a443 | 902 | } |
32d8f622 | 903 | if (returnType==1) { |
904 | stat = projection->GetRMS(); | |
905 | err = projection->GetRMSError(); | |
21f3a443 | 906 | } |
32d8f622 | 907 | if (returnType==2 || returnType==3){ |
8ecd39c2 | 908 | if (returnType==2) {stat= vecLTM[1]; err =projection->GetRMSError();} |
909 | if (returnType==3) {stat= vecLTM[2]; err =projection->GetRMSError();} | |
32d8f622 | 910 | } |
911 | if (returnType==4|| returnType==5){ | |
8fae6121 | 912 | projection->Fit(&f1,"QN","QN", vecLTM[7], vecLTM[8]); |
32d8f622 | 913 | if (returnType==4) { |
914 | stat= f1.GetParameter(1); | |
915 | err=f1.GetParError(1); | |
916 | } | |
917 | if (returnType==5) { | |
918 | stat= f1.GetParameter(2); | |
919 | err=f1.GetParError(2); | |
920 | } | |
921 | } | |
922 | vecX[icount]=xcenter; | |
923 | vecY[icount]=stat; | |
924 | vecYErr[icount]=err; | |
21f3a443 | 925 | icount++; |
32d8f622 | 926 | delete projection; |
21f3a443 | 927 | } |
32d8f622 | 928 | TGraphErrors *graph = new TGraphErrors(icount,vecX.GetMatrixArray(), vecY.GetMatrixArray(),0, vecYErr.GetMatrixArray()); |
929 | graph->SetMarkerStyle(markerStyle); | |
930 | graph->SetMarkerColor(markerColor); | |
21f3a443 | 931 | return graph; |
932 | } | |
933 | ||
934 | ||
935 | ||
936 | ||
937 | ||
88b1c775 | 938 | 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 | 939 | // |
940 | // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts | |
941 | // returns chi2, fitParam and covMatrix | |
942 | // returns TString with fitted formula | |
943 | // | |
dd46129c | 944 | |
21f3a443 | 945 | TString formulaStr(formula); |
946 | TString drawStr(drawCommand); | |
947 | TString cutStr(cuts); | |
dd46129c | 948 | TString ferr("1"); |
949 | ||
950 | TString strVal(drawCommand); | |
951 | if (strVal.Contains(":")){ | |
952 | TObjArray* valTokens = strVal.Tokenize(":"); | |
953 | drawStr = valTokens->At(0)->GetName(); | |
954 | ferr = valTokens->At(1)->GetName(); | |
09d5920f | 955 | delete valTokens; |
dd46129c | 956 | } |
957 | ||
21f3a443 | 958 | |
959 | formulaStr.ReplaceAll("++", "~"); | |
960 | TObjArray* formulaTokens = formulaStr.Tokenize("~"); | |
961 | Int_t dim = formulaTokens->GetEntriesFast(); | |
962 | ||
963 | fitParam.ResizeTo(dim); | |
964 | covMatrix.ResizeTo(dim,dim); | |
965 | ||
966 | TLinearFitter* fitter = new TLinearFitter(dim+1, Form("hyp%d",dim)); | |
967 | fitter->StoreData(kTRUE); | |
968 | fitter->ClearPoints(); | |
969 | ||
970 | Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start, start); | |
09d5920f | 971 | if (entries == -1) { |
972 | delete formulaTokens; | |
5d1391ec | 973 | return new TString(TString::Format("ERROR expr: %s\t%s\tEntries==0",drawStr.Data(),cutStr.Data())); |
09d5920f | 974 | } |
bd7b4d18 | 975 | Double_t **values = new Double_t*[dim+1] ; |
976 | for (Int_t i=0; i<dim+1; i++) values[i]=NULL; | |
dd46129c | 977 | // |
978 | entries = tree->Draw(ferr.Data(), cutStr.Data(), "goff", stop-start, start); | |
b8072cce | 979 | if (entries == -1) { |
09d5920f | 980 | delete formulaTokens; |
b8072cce | 981 | delete []values; |
5d1391ec | 982 | return new TString(TString::Format("ERROR error part: %s\t%s\tEntries==0",ferr.Data(),cutStr.Data())); |
b8072cce | 983 | } |
dd46129c | 984 | Double_t *errors = new Double_t[entries]; |
985 | memcpy(errors, tree->GetV1(), entries*sizeof(Double_t)); | |
21f3a443 | 986 | |
987 | for (Int_t i = 0; i < dim + 1; i++){ | |
988 | Int_t centries = 0; | |
989 | if (i < dim) centries = tree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start); | |
990 | else centries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start); | |
991 | ||
b8072cce | 992 | if (entries != centries) { |
993 | delete []errors; | |
994 | delete []values; | |
5d1391ec | 995 | return new TString(TString::Format("ERROR: %s\t%s\tEntries==%d\tEntries2=%d\n",drawStr.Data(),cutStr.Data(),entries,centries)); |
b8072cce | 996 | } |
21f3a443 | 997 | values[i] = new Double_t[entries]; |
998 | memcpy(values[i], tree->GetV1(), entries*sizeof(Double_t)); | |
999 | } | |
1000 | ||
1001 | // add points to the fitter | |
1002 | for (Int_t i = 0; i < entries; i++){ | |
1003 | Double_t x[1000]; | |
1004 | for (Int_t j=0; j<dim;j++) x[j]=values[j][i]; | |
dd46129c | 1005 | fitter->AddPoint(x, values[dim][i], errors[i]); |
21f3a443 | 1006 | } |
1007 | ||
1008 | fitter->Eval(); | |
2c629c56 | 1009 | if (frac>0.5 && frac<1){ |
1010 | fitter->EvalRobust(frac); | |
88b1c775 | 1011 | }else{ |
1012 | if (fix0) { | |
1013 | fitter->FixParameter(0,0); | |
1014 | fitter->Eval(); | |
1015 | } | |
2c629c56 | 1016 | } |
21f3a443 | 1017 | fitter->GetParameters(fitParam); |
1018 | fitter->GetCovarianceMatrix(covMatrix); | |
1019 | chi2 = fitter->GetChisquare(); | |
b8072cce | 1020 | npoints = entries; |
21f3a443 | 1021 | TString *preturnFormula = new TString(Form("( %f+",fitParam[0])), &returnFormula = *preturnFormula; |
1022 | ||
1023 | for (Int_t iparam = 0; iparam < dim; iparam++) { | |
1024 | returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam+1])); | |
1025 | if (iparam < dim-1) returnFormula.Append("+"); | |
1026 | } | |
1027 | returnFormula.Append(" )"); | |
4d61c301 | 1028 | |
1029 | ||
b8072cce | 1030 | for (Int_t j=0; j<dim+1;j++) delete [] values[j]; |
4d61c301 | 1031 | |
1032 | ||
cb1d20de | 1033 | delete formulaTokens; |
1034 | delete fitter; | |
1035 | delete[] values; | |
b8072cce | 1036 | delete[] errors; |
cb1d20de | 1037 | return preturnFormula; |
1038 | } | |
1039 | ||
1040 | 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){ | |
1041 | // | |
1042 | // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts | |
1043 | // returns chi2, fitParam and covMatrix | |
1044 | // returns TString with fitted formula | |
1045 | // | |
1046 | ||
1047 | TString formulaStr(formula); | |
1048 | TString drawStr(drawCommand); | |
1049 | TString cutStr(cuts); | |
1050 | TString ferr("1"); | |
1051 | ||
1052 | TString strVal(drawCommand); | |
1053 | if (strVal.Contains(":")){ | |
1054 | TObjArray* valTokens = strVal.Tokenize(":"); | |
1055 | drawStr = valTokens->At(0)->GetName(); | |
1056 | ferr = valTokens->At(1)->GetName(); | |
09d5920f | 1057 | delete valTokens; |
cb1d20de | 1058 | } |
1059 | ||
1060 | ||
1061 | formulaStr.ReplaceAll("++", "~"); | |
1062 | TObjArray* formulaTokens = formulaStr.Tokenize("~"); | |
1063 | Int_t dim = formulaTokens->GetEntriesFast(); | |
1064 | ||
1065 | fitParam.ResizeTo(dim); | |
1066 | covMatrix.ResizeTo(dim,dim); | |
1067 | ||
1068 | TLinearFitter* fitter = new TLinearFitter(dim+1, Form("hyp%d",dim)); | |
1069 | fitter->StoreData(kTRUE); | |
1070 | fitter->ClearPoints(); | |
1071 | ||
1072 | Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start, start); | |
09d5920f | 1073 | if (entries == -1) { |
1074 | delete formulaTokens; | |
1075 | return new TString("An ERROR has occured during fitting!"); | |
1076 | } | |
cb1d20de | 1077 | Double_t **values = new Double_t*[dim+1] ; |
bd7b4d18 | 1078 | for (Int_t i=0; i<dim+1; i++) values[i]=NULL; |
cb1d20de | 1079 | // |
1080 | entries = tree->Draw(ferr.Data(), cutStr.Data(), "goff", stop-start, start); | |
b8072cce | 1081 | if (entries == -1) { |
09d5920f | 1082 | delete formulaTokens; |
b8072cce | 1083 | delete [] values; |
1084 | return new TString("An ERROR has occured during fitting!"); | |
1085 | } | |
cb1d20de | 1086 | Double_t *errors = new Double_t[entries]; |
1087 | memcpy(errors, tree->GetV1(), entries*sizeof(Double_t)); | |
1088 | ||
1089 | for (Int_t i = 0; i < dim + 1; i++){ | |
1090 | Int_t centries = 0; | |
1091 | if (i < dim) centries = tree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start); | |
1092 | else centries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start); | |
1093 | ||
b8072cce | 1094 | if (entries != centries) { |
1095 | delete []errors; | |
1096 | delete []values; | |
09d5920f | 1097 | delete formulaTokens; |
b8072cce | 1098 | return new TString("An ERROR has occured during fitting!"); |
1099 | } | |
cb1d20de | 1100 | values[i] = new Double_t[entries]; |
1101 | memcpy(values[i], tree->GetV1(), entries*sizeof(Double_t)); | |
1102 | } | |
1103 | ||
1104 | // add points to the fitter | |
1105 | for (Int_t i = 0; i < entries; i++){ | |
1106 | Double_t x[1000]; | |
1107 | for (Int_t j=0; j<dim;j++) x[j]=values[j][i]; | |
1108 | fitter->AddPoint(x, values[dim][i], errors[i]); | |
1109 | } | |
1110 | if (constrain>0){ | |
1111 | for (Int_t i = 0; i < dim; i++){ | |
1112 | Double_t x[1000]; | |
1113 | for (Int_t j=0; j<dim;j++) if (i!=j) x[j]=0; | |
1114 | x[i]=1.; | |
1115 | fitter->AddPoint(x, 0, constrain); | |
1116 | } | |
1117 | } | |
1118 | ||
1119 | ||
1120 | fitter->Eval(); | |
1121 | if (frac>0.5 && frac<1){ | |
1122 | fitter->EvalRobust(frac); | |
1123 | } | |
1124 | fitter->GetParameters(fitParam); | |
1125 | fitter->GetCovarianceMatrix(covMatrix); | |
1126 | chi2 = fitter->GetChisquare(); | |
1127 | npoints = entries; | |
cb1d20de | 1128 | |
1129 | TString *preturnFormula = new TString(Form("( %f+",fitParam[0])), &returnFormula = *preturnFormula; | |
1130 | ||
1131 | for (Int_t iparam = 0; iparam < dim; iparam++) { | |
1132 | returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam+1])); | |
1133 | if (iparam < dim-1) returnFormula.Append("+"); | |
1134 | } | |
1135 | returnFormula.Append(" )"); | |
1136 | ||
b8072cce | 1137 | for (Int_t j=0; j<dim+1;j++) delete [] values[j]; |
cb1d20de | 1138 | |
1139 | ||
1140 | ||
1141 | delete formulaTokens; | |
1142 | delete fitter; | |
1143 | delete[] values; | |
b8072cce | 1144 | delete[] errors; |
cb1d20de | 1145 | return preturnFormula; |
1146 | } | |
1147 | ||
1148 | ||
1149 | ||
1150 | 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){ | |
1151 | // | |
1152 | // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts | |
1153 | // returns chi2, fitParam and covMatrix | |
1154 | // returns TString with fitted formula | |
1155 | // | |
1156 | ||
1157 | TString formulaStr(formula); | |
1158 | TString drawStr(drawCommand); | |
1159 | TString cutStr(cuts); | |
1160 | TString ferr("1"); | |
1161 | ||
1162 | TString strVal(drawCommand); | |
1163 | if (strVal.Contains(":")){ | |
1164 | TObjArray* valTokens = strVal.Tokenize(":"); | |
1165 | drawStr = valTokens->At(0)->GetName(); | |
09d5920f | 1166 | ferr = valTokens->At(1)->GetName(); |
1167 | delete valTokens; | |
cb1d20de | 1168 | } |
1169 | ||
1170 | ||
1171 | formulaStr.ReplaceAll("++", "~"); | |
1172 | TObjArray* formulaTokens = formulaStr.Tokenize("~"); | |
1173 | Int_t dim = formulaTokens->GetEntriesFast(); | |
1174 | ||
1175 | fitParam.ResizeTo(dim); | |
1176 | covMatrix.ResizeTo(dim,dim); | |
1177 | TString fitString="x0"; | |
1178 | for (Int_t i=1; i<dim; i++) fitString+=Form("++x%d",i); | |
1179 | TLinearFitter* fitter = new TLinearFitter(dim, fitString.Data()); | |
1180 | fitter->StoreData(kTRUE); | |
1181 | fitter->ClearPoints(); | |
1182 | ||
1183 | Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start, start); | |
09d5920f | 1184 | if (entries == -1) { |
1185 | delete formulaTokens; | |
1186 | return new TString("An ERROR has occured during fitting!"); | |
1187 | } | |
cb1d20de | 1188 | Double_t **values = new Double_t*[dim+1] ; |
bd7b4d18 | 1189 | for (Int_t i=0; i<dim+1; i++) values[i]=NULL; |
cb1d20de | 1190 | // |
1191 | entries = tree->Draw(ferr.Data(), cutStr.Data(), "goff", stop-start, start); | |
b8072cce | 1192 | if (entries == -1) { |
1193 | delete []values; | |
09d5920f | 1194 | delete formulaTokens; |
b8072cce | 1195 | return new TString("An ERROR has occured during fitting!"); |
1196 | } | |
cb1d20de | 1197 | Double_t *errors = new Double_t[entries]; |
1198 | memcpy(errors, tree->GetV1(), entries*sizeof(Double_t)); | |
1199 | ||
1200 | for (Int_t i = 0; i < dim + 1; i++){ | |
1201 | Int_t centries = 0; | |
1202 | if (i < dim) centries = tree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start); | |
1203 | else centries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start); | |
1204 | ||
b8072cce | 1205 | if (entries != centries) { |
1206 | delete []errors; | |
1207 | delete []values; | |
09d5920f | 1208 | delete formulaTokens; |
b8072cce | 1209 | return new TString("An ERROR has occured during fitting!"); |
1210 | } | |
cb1d20de | 1211 | values[i] = new Double_t[entries]; |
1212 | memcpy(values[i], tree->GetV1(), entries*sizeof(Double_t)); | |
1213 | } | |
1214 | ||
1215 | // add points to the fitter | |
1216 | for (Int_t i = 0; i < entries; i++){ | |
1217 | Double_t x[1000]; | |
1218 | for (Int_t j=0; j<dim;j++) x[j]=values[j][i]; | |
1219 | fitter->AddPoint(x, values[dim][i], errors[i]); | |
1220 | } | |
1221 | ||
1222 | fitter->Eval(); | |
1223 | if (frac>0.5 && frac<1){ | |
1224 | fitter->EvalRobust(frac); | |
1225 | } | |
1226 | fitter->GetParameters(fitParam); | |
1227 | fitter->GetCovarianceMatrix(covMatrix); | |
1228 | chi2 = fitter->GetChisquare(); | |
1229 | npoints = entries; | |
cb1d20de | 1230 | |
1231 | TString *preturnFormula = new TString("("), &returnFormula = *preturnFormula; | |
1232 | ||
1233 | for (Int_t iparam = 0; iparam < dim; iparam++) { | |
1234 | returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam])); | |
1235 | if (iparam < dim-1) returnFormula.Append("+"); | |
1236 | } | |
1237 | returnFormula.Append(" )"); | |
1238 | ||
1239 | ||
b8072cce | 1240 | for (Int_t j=0; j<dim+1;j++) delete [] values[j]; |
cb1d20de | 1241 | |
21f3a443 | 1242 | delete formulaTokens; |
1243 | delete fitter; | |
1244 | delete[] values; | |
b8072cce | 1245 | delete[] errors; |
21f3a443 | 1246 | return preturnFormula; |
1247 | } | |
7c9cf6e4 | 1248 | |
1249 | ||
1250 | ||
1251 | ||
1252 | ||
3d7cc0b4 | 1253 | Int_t TStatToolkit::GetFitIndex(const TString fString, const TString subString){ |
7c9cf6e4 | 1254 | // |
1255 | // fitString - ++ separated list of fits | |
1256 | // substring - ++ separated list of the requiered substrings | |
1257 | // | |
1258 | // return the last occurance of substring in fit string | |
1259 | // | |
1260 | TObjArray *arrFit = fString.Tokenize("++"); | |
1261 | TObjArray *arrSub = subString.Tokenize("++"); | |
1262 | Int_t index=-1; | |
1263 | for (Int_t i=0; i<arrFit->GetEntries(); i++){ | |
1264 | Bool_t isOK=kTRUE; | |
1265 | TString str =arrFit->At(i)->GetName(); | |
1266 | for (Int_t isub=0; isub<arrSub->GetEntries(); isub++){ | |
1267 | if (str.Contains(arrSub->At(isub)->GetName())==0) isOK=kFALSE; | |
1268 | } | |
1269 | if (isOK) index=i; | |
1270 | } | |
09d5920f | 1271 | delete arrFit; |
1272 | delete arrSub; | |
7c9cf6e4 | 1273 | return index; |
1274 | } | |
1275 | ||
1276 | ||
3d7cc0b4 | 1277 | TString TStatToolkit::FilterFit(const TString &input, const TString filter, TVectorD ¶m, TMatrixD & covar){ |
7c9cf6e4 | 1278 | // |
1279 | // Filter fit expression make sub-fit | |
1280 | // | |
1281 | TObjArray *array0= input.Tokenize("++"); | |
1282 | TObjArray *array1= filter.Tokenize("++"); | |
1283 | //TString *presult=new TString("(0"); | |
1284 | TString result="(0.0"; | |
1285 | for (Int_t i=0; i<array0->GetEntries(); i++){ | |
1286 | Bool_t isOK=kTRUE; | |
1287 | TString str(array0->At(i)->GetName()); | |
1288 | for (Int_t j=0; j<array1->GetEntries(); j++){ | |
1289 | if (str.Contains(array1->At(j)->GetName())==0) isOK=kFALSE; | |
1290 | } | |
1291 | if (isOK) { | |
1292 | result+="+"+str; | |
1293 | result+=Form("*(%f)",param[i+1]); | |
1294 | printf("%f\t%f\t%s\n",param[i+1], TMath::Sqrt(covar(i+1,i+1)),str.Data()); | |
1295 | } | |
1296 | } | |
1297 | result+="-0.)"; | |
09d5920f | 1298 | delete array0; |
1299 | delete array1; | |
7c9cf6e4 | 1300 | return result; |
1301 | } | |
1302 | ||
1303 | void TStatToolkit::Update1D(Double_t delta, Double_t sigma, Int_t s1, TMatrixD &vecXk, TMatrixD &covXk){ | |
1304 | // | |
1305 | // Update parameters and covariance - with one measurement | |
1306 | // Input: | |
1307 | // vecXk - input vector - Updated in function | |
1308 | // covXk - covariance matrix - Updated in function | |
1309 | // delta, sigma, s1 - new measurement, rms of new measurement and the index of measurement | |
1310 | const Int_t knMeas=1; | |
1311 | Int_t knElem=vecXk.GetNrows(); | |
1312 | ||
1313 | TMatrixD mat1(knElem,knElem); // update covariance matrix | |
1314 | TMatrixD matHk(1,knElem); // vector to mesurement | |
1315 | TMatrixD vecYk(knMeas,1); // Innovation or measurement residual | |
1316 | TMatrixD matHkT(knElem,knMeas); // helper matrix Hk transpose | |
1317 | TMatrixD matSk(knMeas,knMeas); // Innovation (or residual) covariance | |
1318 | TMatrixD matKk(knElem,knMeas); // Optimal Kalman gain | |
1319 | TMatrixD covXk2(knElem,knElem); // helper matrix | |
1320 | TMatrixD covXk3(knElem,knElem); // helper matrix | |
1321 | TMatrixD vecZk(1,1); | |
1322 | TMatrixD measR(1,1); | |
1323 | vecZk(0,0)=delta; | |
1324 | measR(0,0)=sigma*sigma; | |
1325 | // | |
1326 | // reset matHk | |
1327 | for (Int_t iel=0;iel<knElem;iel++) | |
1328 | for (Int_t ip=0;ip<knMeas;ip++) matHk(ip,iel)=0; | |
1329 | //mat1 | |
1330 | for (Int_t iel=0;iel<knElem;iel++) { | |
1331 | for (Int_t jel=0;jel<knElem;jel++) mat1(iel,jel)=0; | |
1332 | mat1(iel,iel)=1; | |
1333 | } | |
1334 | // | |
1335 | matHk(0, s1)=1; | |
1336 | vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual | |
1337 | matHkT=matHk.T(); matHk.T(); | |
1338 | matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance | |
1339 | matSk.Invert(); | |
1340 | matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain | |
1341 | vecXk += matKk*vecYk; // updated vector | |
1342 | covXk2= (mat1-(matKk*matHk)); | |
1343 | covXk3 = covXk2*covXk; | |
1344 | covXk = covXk3; | |
1345 | Int_t nrows=covXk3.GetNrows(); | |
1346 | ||
1347 | for (Int_t irow=0; irow<nrows; irow++) | |
1348 | for (Int_t icol=0; icol<nrows; icol++){ | |
1349 | // rounding problems - make matrix again symteric | |
1350 | covXk(irow,icol)=(covXk3(irow,icol)+covXk3(icol,irow))*0.5; | |
1351 | } | |
1352 | } | |
1353 | ||
1354 | ||
1355 | ||
3d7cc0b4 | 1356 | void TStatToolkit::Constrain1D(const TString &input, const TString filter, TVectorD ¶m, TMatrixD & covar, Double_t mean, Double_t sigma){ |
7c9cf6e4 | 1357 | // |
1358 | // constrain linear fit | |
1359 | // input - string description of fit function | |
1360 | // filter - string filter to select sub fits | |
1361 | // param,covar - parameters and covariance matrix of the fit | |
1362 | // mean,sigma - new measurement uning which the fit is updated | |
1363 | // | |
ae45c94d | 1364 | |
7c9cf6e4 | 1365 | TObjArray *array0= input.Tokenize("++"); |
1366 | TObjArray *array1= filter.Tokenize("++"); | |
1367 | TMatrixD paramM(param.GetNrows(),1); | |
1368 | for (Int_t i=0; i<=array0->GetEntries(); i++){paramM(i,0)=param(i);} | |
1369 | ||
ae45c94d | 1370 | if (filter.Length()==0){ |
1371 | TStatToolkit::Update1D(mean, sigma, 0, paramM, covar);// | |
1372 | }else{ | |
1373 | for (Int_t i=0; i<array0->GetEntries(); i++){ | |
1374 | Bool_t isOK=kTRUE; | |
1375 | TString str(array0->At(i)->GetName()); | |
1376 | for (Int_t j=0; j<array1->GetEntries(); j++){ | |
1377 | if (str.Contains(array1->At(j)->GetName())==0) isOK=kFALSE; | |
1378 | } | |
1379 | if (isOK) { | |
1380 | TStatToolkit::Update1D(mean, sigma, i+1, paramM, covar);// | |
1381 | } | |
7c9cf6e4 | 1382 | } |
1383 | } | |
1384 | for (Int_t i=0; i<=array0->GetEntries(); i++){ | |
1385 | param(i)=paramM(i,0); | |
1386 | } | |
09d5920f | 1387 | delete array0; |
1388 | delete array1; | |
7c9cf6e4 | 1389 | } |
1390 | ||
ae45c94d | 1391 | TString TStatToolkit::MakeFitString(const TString &input, const TVectorD ¶m, const TMatrixD & covar, Bool_t verbose){ |
7c9cf6e4 | 1392 | // |
1393 | // | |
1394 | // | |
1395 | TObjArray *array0= input.Tokenize("++"); | |
ae45c94d | 1396 | TString result=Form("(%f",param[0]); |
1397 | printf("%f\t%f\t\n", param[0], TMath::Sqrt(covar(0,0))); | |
7c9cf6e4 | 1398 | for (Int_t i=0; i<array0->GetEntries(); i++){ |
1399 | TString str(array0->At(i)->GetName()); | |
1400 | result+="+"+str; | |
1401 | result+=Form("*(%f)",param[i+1]); | |
ae45c94d | 1402 | if (verbose) printf("%f\t%f\t%s\n", param[i+1], TMath::Sqrt(covar(i+1,i+1)),str.Data()); |
7c9cf6e4 | 1403 | } |
1404 | result+="-0.)"; | |
09d5920f | 1405 | delete array0; |
7c9cf6e4 | 1406 | return result; |
1407 | } | |
df0a2a0a | 1408 | |
1a8f4649 | 1409 | TGraphErrors * TStatToolkit::MakeGraphErrors(TTree * tree, const char * expr, const char * cut, Int_t mstyle, Int_t mcolor, Float_t msize, Float_t offset){ |
1410 | // | |
1411 | // Query a graph errors | |
1412 | // return TGraphErrors specified by expr and cut | |
1413 | // Example usage TStatToolkit::MakeGraphError(tree,"Y:X:ErrY","X>0", 25,2,0.4) | |
1414 | // tree - tree with variable | |
1415 | // expr - examp | |
1416 | const Int_t entries = tree->Draw(expr,cut,"goff"); | |
1417 | if (entries<=0) { | |
1418 | TStatToolkit t; | |
1419 | t.Error("TStatToolkit::MakeGraphError",Form("Empty or Not valid expression (%s) or cut *%s)", expr,cut)); | |
1420 | return 0; | |
1421 | } | |
1422 | if ( tree->GetV2()==0){ | |
1423 | TStatToolkit t; | |
1424 | t.Error("TStatToolkit::MakeGraphError",Form("Not valid expression (%s) ", expr)); | |
1425 | return 0; | |
1426 | } | |
1427 | TGraphErrors * graph=0; | |
1428 | if ( tree->GetV3()!=0){ | |
1429 | graph = new TGraphErrors (entries, tree->GetV2(),tree->GetV1(),0,tree->GetV3()); | |
1430 | }else{ | |
1431 | graph = new TGraphErrors (entries, tree->GetV2(),tree->GetV1(),0,0); | |
1432 | } | |
1433 | graph->SetMarkerStyle(mstyle); | |
1434 | graph->SetMarkerColor(mcolor); | |
1435 | graph->SetLineColor(mcolor); | |
5d6816d1 | 1436 | graph->SetTitle(expr); |
1437 | TString chstring(expr); | |
1438 | TObjArray *charray = chstring.Tokenize(":"); | |
1439 | graph->GetXaxis()->SetTitle(charray->At(1)->GetName()); | |
1440 | graph->GetYaxis()->SetTitle(charray->At(0)->GetName()); | |
1441 | delete charray; | |
1a8f4649 | 1442 | if (msize>0) graph->SetMarkerSize(msize); |
1443 | for(Int_t i=0;i<graph->GetN();i++) graph->GetX()[i]+=offset; | |
1444 | return graph; | |
1445 | ||
1446 | } | |
1447 | ||
df0a2a0a | 1448 | |
377a7d60 | 1449 | TGraph * TStatToolkit::MakeGraphSparse(TTree * tree, const char * expr, const char * cut, Int_t mstyle, Int_t mcolor, Float_t msize, Float_t offset){ |
df0a2a0a | 1450 | // |
1451 | // Make a sparse draw of the variables | |
d02b8756 | 1452 | // Format of expr : Var:Run or Var:Run:ErrorY or Var:Run:ErrorY:ErrorX |
1453 | // offset : points can slightly be shifted in x for better visibility with more graphs | |
1454 | // | |
1455 | // Written by Weilin.Yu | |
1456 | // updated & merged with QA-code by Patrick Reichelt | |
1457 | // | |
1458 | const Int_t entries = tree->Draw(expr,cut,"goff"); | |
8fb3bea1 | 1459 | if (entries<=0) { |
1460 | TStatToolkit t; | |
d02b8756 | 1461 | t.Error("TStatToolkit::MakeGraphSparse",Form("Empty or Not valid expression (%s) or cut (%s)", expr, cut)); |
8fb3bea1 | 1462 | return 0; |
1463 | } | |
df0a2a0a | 1464 | // TGraph * graph = (TGraph*)gPad->GetPrimitive("Graph"); // 2D |
d02b8756 | 1465 | |
1466 | Double_t *graphY, *graphX; | |
1467 | graphY = tree->GetV1(); | |
1468 | graphX = tree->GetV2(); | |
1469 | ||
1470 | // sort according to run number | |
eda18694 | 1471 | Int_t *index = new Int_t[entries*4]; |
d02b8756 | 1472 | TMath::Sort(entries,graphX,index,kFALSE); |
df0a2a0a | 1473 | |
d02b8756 | 1474 | // define arrays for the new graph |
1475 | Double_t *unsortedX = new Double_t[entries]; | |
1476 | Int_t *runNumber = new Int_t[entries]; | |
df0a2a0a | 1477 | Double_t count = 0.5; |
d02b8756 | 1478 | |
1479 | // evaluate arrays for the new graph according to the run-number | |
baa0041d | 1480 | Int_t icount=0; |
d02b8756 | 1481 | //first entry |
1482 | unsortedX[index[0]] = count; | |
1483 | runNumber[0] = graphX[index[0]]; | |
1484 | // loop the rest of entries | |
1485 | for(Int_t i=1;i<entries;i++) | |
1486 | { | |
1487 | if(graphX[index[i]]==graphX[index[i-1]]) | |
1488 | unsortedX[index[i]] = count; | |
1489 | else if(graphX[index[i]]!=graphX[index[i-1]]){ | |
df0a2a0a | 1490 | count++; |
baa0041d | 1491 | icount++; |
d02b8756 | 1492 | unsortedX[index[i]] = count; |
1493 | runNumber[icount]=graphX[index[i]]; | |
df0a2a0a | 1494 | } |
1495 | } | |
d02b8756 | 1496 | |
1497 | // count the number of xbins (run-wise) for the new graph | |
df0a2a0a | 1498 | const Int_t newNbins = int(count+0.5); |
1499 | Double_t *newBins = new Double_t[newNbins+1]; | |
1500 | for(Int_t i=0; i<=count+1;i++){ | |
1501 | newBins[i] = i; | |
1502 | } | |
d02b8756 | 1503 | |
1504 | // define and fill the new graph | |
b3453fe7 | 1505 | TGraph *graphNew = 0; |
d02b8756 | 1506 | if (tree->GetV3()) { |
1507 | if (tree->GetV4()) { | |
1508 | graphNew = new TGraphErrors(entries,unsortedX,graphY,tree->GetV4(),tree->GetV3()); | |
1509 | } | |
1510 | else { graphNew = new TGraphErrors(entries,unsortedX,graphY,0,tree->GetV3()); } | |
1511 | } | |
1512 | else { graphNew = new TGraphErrors(entries,unsortedX,graphY,0,0); } | |
1513 | // with "Set(...)", the x-axis is being sorted | |
df0a2a0a | 1514 | graphNew->GetXaxis()->Set(newNbins,newBins); |
d02b8756 | 1515 | |
1516 | // set the bins for the x-axis, apply shifting of points | |
df0a2a0a | 1517 | Char_t xName[50]; |
df0a2a0a | 1518 | for(Int_t i=0;i<count;i++){ |
d02b8756 | 1519 | snprintf(xName,50,"%d",runNumber[i]); |
df0a2a0a | 1520 | graphNew->GetXaxis()->SetBinLabel(i+1,xName); |
d02b8756 | 1521 | graphNew->GetX()[i]+=offset; |
df0a2a0a | 1522 | } |
d02b8756 | 1523 | |
df0a2a0a | 1524 | graphNew->GetHistogram()->SetTitle(""); |
d02b8756 | 1525 | graphNew->SetMarkerStyle(mstyle); |
fc03fa55 | 1526 | graphNew->SetMarkerColor(mcolor); graphNew->SetLineColor(mcolor); |
1527 | if (msize>0) { graphNew->SetMarkerSize(msize); graphNew->SetLineWidth(msize); } | |
d02b8756 | 1528 | delete [] unsortedX; |
1529 | delete [] runNumber; | |
df0a2a0a | 1530 | delete [] index; |
1531 | delete [] newBins; | |
5d6816d1 | 1532 | // |
1533 | graphNew->SetTitle(expr); | |
1534 | TString chstring(expr); | |
1535 | TObjArray *charray = chstring.Tokenize(":"); | |
1536 | graphNew->GetXaxis()->SetTitle(charray->At(1)->GetName()); | |
1537 | graphNew->GetYaxis()->SetTitle(charray->At(0)->GetName()); | |
1538 | delete charray; | |
df0a2a0a | 1539 | return graphNew; |
1540 | } | |
1541 | ||
377a7d60 | 1542 | |
1543 | ||
1544 | // | |
d02b8756 | 1545 | // functions used for the trending |
377a7d60 | 1546 | // |
1547 | ||
1548 | Int_t TStatToolkit::MakeStatAlias(TTree * tree, const char * expr, const char * cut, const char * alias) | |
1549 | { | |
1550 | // | |
1551 | // Add alias using statistical values of a given variable. | |
1552 | // (by MI, Patrick Reichelt) | |
1553 | // | |
1554 | // tree - input tree | |
1555 | // expr - variable expression | |
1556 | // cut - selection criteria | |
1557 | // Output - return number of entries used to define variable | |
1558 | // In addition mean, rms, median, and robust mean and rms (choosing fraction of data with smallest RMS) | |
1559 | // | |
d02b8756 | 1560 | /* Example usage: |
1561 | 1.) create the robust estimators for variable expr="QA.TPC.CPass1.meanTPCncl" and create a corresponding | |
1562 | aliases with the prefix alias[0]="ncl", calculated using fraction alias[1]="0.90" | |
1563 | ||
1564 | TStatToolkit::MakeStatAlias(tree,"QA.TPC.CPass1.meanTPCncl","QA.TPC.CPass1.status>0","ncl:0.9"); | |
1565 | root [4] tree->GetListOfAliases().Print() | |
1566 | OBJ: TNamed ncl_Median (130.964333+0) | |
1567 | OBJ: TNamed ncl_Mean (122.120387+0) | |
1568 | OBJ: TNamed ncl_RMS (33.509623+0) | |
1569 | OBJ: TNamed ncl_Mean90 (131.503862+0) | |
1570 | OBJ: TNamed ncl_RMS90 (3.738260+0) | |
1571 | */ | |
377a7d60 | 1572 | // |
d02b8756 | 1573 | Int_t entries = tree->Draw(expr,cut,"goff"); |
377a7d60 | 1574 | if (entries<=1){ |
1575 | printf("Expression or cut not valid:\t%s\t%s\n", expr, cut); | |
1576 | return 0; | |
1577 | } | |
1578 | // | |
1579 | TObjArray* oaAlias = TString(alias).Tokenize(":"); | |
fc03fa55 | 1580 | if (oaAlias->GetEntries()<2) { |
1581 | printf("Alias must have 2 arguments:\t%s\n", alias); | |
1582 | return 0; | |
1583 | } | |
377a7d60 | 1584 | Float_t entryFraction = atof( oaAlias->At(1)->GetName() ); |
1585 | // | |
1586 | Double_t median = TMath::Median(entries,tree->GetV1()); | |
d02b8756 | 1587 | Double_t mean = TMath::Mean(entries,tree->GetV1()); |
377a7d60 | 1588 | Double_t rms = TMath::RMS(entries,tree->GetV1()); |
1589 | Double_t meanEF=0, rmsEF=0; | |
1590 | TStatToolkit::EvaluateUni(entries, tree->GetV1(), meanEF, rmsEF, entries*entryFraction); | |
1591 | // | |
d02b8756 | 1592 | tree->SetAlias(Form("%s_Median",oaAlias->At(0)->GetName()), Form("(%f+0)",median)); |
377a7d60 | 1593 | tree->SetAlias(Form("%s_Mean",oaAlias->At(0)->GetName()), Form("(%f+0)",mean)); |
1594 | tree->SetAlias(Form("%s_RMS",oaAlias->At(0)->GetName()), Form("(%f+0)",rms)); | |
377a7d60 | 1595 | tree->SetAlias(Form("%s_Mean%d",oaAlias->At(0)->GetName(),Int_t(entryFraction*100)), Form("(%f+0)",meanEF)); |
1596 | tree->SetAlias(Form("%s_RMS%d",oaAlias->At(0)->GetName(),Int_t(entryFraction*100)), Form("(%f+0)",rmsEF)); | |
1597 | delete oaAlias; | |
1598 | return entries; | |
1599 | } | |
1600 | ||
1601 | Int_t TStatToolkit::SetStatusAlias(TTree * tree, const char * expr, const char * cut, const char * alias) | |
1602 | { | |
1603 | // | |
1604 | // Add alias to trending tree using statistical values of a given variable. | |
1605 | // (by MI, Patrick Reichelt) | |
1606 | // | |
1607 | // format of expr : varname (e.g. meanTPCncl) | |
1608 | // format of cut : char like in TCut | |
1609 | // format of alias: alias:query:entryFraction(EF) (fraction of entries used for uniformity evaluation) | |
1610 | // e.g.: varname_Out:(abs(varname-meanEF)>6.*rmsEF):0.8 | |
d02b8756 | 1611 | // available internal variables are: 'varname, Median, Mean, MeanEF, RMS, RMSEF' |
1612 | // in the alias, 'varname' will be replaced by its content, and 'EF' by the percentage (e.g. MeanEF -> Mean80) | |
377a7d60 | 1613 | // |
1614 | /* Example usage: | |
d02b8756 | 1615 | 1.) Define robust mean (possible, but easier done with TStatToolkit::MakeStatAlias(...)) |
1616 | TStatToolkit::SetStatusAlias(tree, "meanTPCnclF", "meanTPCnclF>0", "meanTPCnclF_MeanEF:MeanEF:0.80") ; | |
377a7d60 | 1617 | root [10] tree->GetListOfAliases()->Print() |
1618 | Collection name='TList', class='TList', size=1 | |
d02b8756 | 1619 | OBJ: TNamed meanTPCnclF_Mean80 0.899308 |
377a7d60 | 1620 | 2.) create alias outlyers - 6 sigma cut |
d02b8756 | 1621 | TStatToolkit::SetStatusAlias(tree, "meanTPCnclF", "meanTPCnclF>0", "meanTPCnclF_Out:(abs(meanTPCnclF-MeanEF)>6.*RMSEF):0.8") |
1622 | meanTPCnclF_Out ==> (abs(meanTPCnclF-0.899308)>6.*0.016590) | |
377a7d60 | 1623 | 3.) the same functionality as in 2.) |
d02b8756 | 1624 | TStatToolkit::SetStatusAlias(tree, "meanTPCnclF", "meanTPCnclF>0", "varname_Out2:(abs(varname-MeanEF)>6.*RMSEF):0.8") |
1625 | meanTPCnclF_Out2 ==> (abs(meanTPCnclF-0.899308)>6.*0.016590) | |
377a7d60 | 1626 | */ |
1627 | // | |
d02b8756 | 1628 | Int_t entries = tree->Draw(expr,cut,"goff"); |
e0bb28a2 | 1629 | if (entries<1){ |
d02b8756 | 1630 | printf("Expression or cut not valid:\t%s\t%s\n", expr, cut); |
1631 | return 0; | |
1632 | } | |
1633 | // | |
377a7d60 | 1634 | TObjArray* oaVar = TString(expr).Tokenize(":"); |
1635 | char varname[50]; | |
377a7d60 | 1636 | snprintf(varname,50,"%s", oaVar->At(0)->GetName()); |
fc03fa55 | 1637 | Float_t entryFraction = 0.8; |
d02b8756 | 1638 | // |
377a7d60 | 1639 | TObjArray* oaAlias = TString(alias).Tokenize(":"); |
fc03fa55 | 1640 | if (oaAlias->GetEntries()<2) { |
1641 | printf("Alias must have at least 2 arguments:\t%s\n", alias); | |
1642 | return 0; | |
1643 | } | |
1644 | else if (oaAlias->GetEntries()<3) { | |
1645 | //printf("Using default entryFraction if needed:\t%f\n", entryFraction); | |
1646 | } | |
1647 | else entryFraction = atof( oaAlias->At(2)->GetName() ); | |
377a7d60 | 1648 | // |
377a7d60 | 1649 | Double_t median = TMath::Median(entries,tree->GetV1()); |
d02b8756 | 1650 | Double_t mean = TMath::Mean(entries,tree->GetV1()); |
377a7d60 | 1651 | Double_t rms = TMath::RMS(entries,tree->GetV1()); |
1652 | Double_t meanEF=0, rmsEF=0; | |
1653 | TStatToolkit::EvaluateUni(entries, tree->GetV1(), meanEF, rmsEF, entries*entryFraction); | |
377a7d60 | 1654 | // |
1655 | TString sAlias( oaAlias->At(0)->GetName() ); | |
1656 | sAlias.ReplaceAll("varname",varname); | |
d02b8756 | 1657 | sAlias.ReplaceAll("MeanEF", Form("Mean%1.0f",entryFraction*100) ); |
1658 | sAlias.ReplaceAll("RMSEF", Form("RMS%1.0f",entryFraction*100) ); | |
377a7d60 | 1659 | TString sQuery( oaAlias->At(1)->GetName() ); |
1660 | sQuery.ReplaceAll("varname",varname); | |
1661 | sQuery.ReplaceAll("MeanEF", Form("%f",meanEF) ); | |
d02b8756 | 1662 | sQuery.ReplaceAll("RMSEF", Form("%f",rmsEF) ); //make sure to replace 'RMSEF' before 'RMS'... |
377a7d60 | 1663 | sQuery.ReplaceAll("Median", Form("%f",median) ); |
d02b8756 | 1664 | sQuery.ReplaceAll("Mean", Form("%f",mean) ); |
377a7d60 | 1665 | sQuery.ReplaceAll("RMS", Form("%f",rms) ); |
d02b8756 | 1666 | printf("define alias:\t%s = %s\n", sAlias.Data(), sQuery.Data()); |
377a7d60 | 1667 | // |
1668 | char query[200]; | |
1669 | char aname[200]; | |
1670 | snprintf(query,200,"%s", sQuery.Data()); | |
1671 | snprintf(aname,200,"%s", sAlias.Data()); | |
1672 | tree->SetAlias(aname, query); | |
d02b8756 | 1673 | delete oaVar; |
1674 | delete oaAlias; | |
377a7d60 | 1675 | return entries; |
1676 | } | |
1677 | ||
1678 | TMultiGraph* TStatToolkit::MakeStatusMultGr(TTree * tree, const char * expr, const char * cut, const char * alias, Int_t igr) | |
1679 | { | |
1680 | // | |
1681 | // Compute a trending multigraph that shows for which runs a variable has outliers. | |
1682 | // (by MI, Patrick Reichelt) | |
1683 | // | |
fc03fa55 | 1684 | // format of expr : varname:xaxis (e.g. meanTPCncl:run, but 'varname' can be any string that you need for seach-and-replace) |
377a7d60 | 1685 | // format of cut : char like in TCut |
fc03fa55 | 1686 | // format of alias: (1):(statisticOK):(varname_Warning):(varname_Out)[:(varname_PhysAcc):(varname_Extra)] |
1687 | // | |
1688 | // function MakeGraphSparse() is called for each alias argument, which will be used as tree expression. | |
1689 | // each alias argument is supposed to be a Boolean statement which can be evaluated as tree expression. | |
1690 | // the order of these criteria should be kept, as the marker styles and colors are chosen to be meaningful this way! | |
1691 | // 'statisticOK' could e.g. be an alias for '(meanTPCncl>0)'. | |
1692 | // if you dont need e.g. a 'warning' condition, then just replace it by (0). | |
377a7d60 | 1693 | // in the alias, 'varname' will be replaced by its content (e.g. varname_Out -> meanTPCncl_Out) |
d02b8756 | 1694 | // note: the aliases 'varname_Out' etc have to be defined by function TStatToolkit::SetStatusAlias(...) |
377a7d60 | 1695 | // counter igr is used to shift the multigraph in y when filling a TObjArray. |
1696 | // | |
92f227d3 | 1697 | // |
1698 | // To create the Status Bar, the following is done in principle. | |
1699 | // ( example current usage in $ALICE_ROOT/PWGPP/TPC/macros/drawPerformanceTPCQAMatchTrends.C and ./qaConfig.C. ) | |
1700 | // | |
1701 | // TStatToolkit::SetStatusAlias(tree, "meanTPCncl", "", "varname_Out:(abs(varname-MeanEF)>6.*RMSEF):0.8"); | |
1702 | // TStatToolkit::SetStatusAlias(tree, "tpcItsMatchA", "", "varname_Out:(abs(varname-MeanEF)>6.*RMSEF):0.8"); | |
1703 | // TStatToolkit::SetStatusAlias(tree, "meanTPCncl", "", "varname_Warning:(abs(varname-MeanEF)>3.*RMSEF):0.8"); | |
1704 | // TStatToolkit::SetStatusAlias(tree, "tpcItsMatchA", "", "varname_Warning:(abs(varname-MeanEF)>3.*RMSEF):0.8"); | |
1705 | // TObjArray* oaMultGr = new TObjArray(); int igr=0; | |
1706 | // oaMultGr->Add( TStatToolkit::MakeStatusMultGr(tree, "tpcItsMatchA:run", "", "(1):(meanTPCncl>0):(varname_Warning):(varname_Outlier):", igr) ); igr++; | |
1707 | // oaMultGr->Add( TStatToolkit::MakeStatusMultGr(tree, "meanTPCncl:run", "", "(1):(meanTPCncl>0):(varname_Warning):(varname_Outlier):", igr) ); igr++; | |
1708 | // TCanvas *c1 = new TCanvas("c1","c1"); | |
1709 | // TStatToolkit::AddStatusPad(c1, 0.30, 0.40); | |
1710 | // TStatToolkit::DrawStatusGraphs(oaMultGr); | |
1711 | ||
1712 | ||
377a7d60 | 1713 | TObjArray* oaVar = TString(expr).Tokenize(":"); |
fc03fa55 | 1714 | if (oaVar->GetEntries()<2) { |
1715 | printf("Expression has to be of type 'varname:xaxis':\t%s\n", expr); | |
1716 | return 0; | |
1717 | } | |
377a7d60 | 1718 | char varname[50]; |
1719 | char var_x[50]; | |
1720 | snprintf(varname,50,"%s", oaVar->At(0)->GetName()); | |
1721 | snprintf(var_x ,50,"%s", oaVar->At(1)->GetName()); | |
d02b8756 | 1722 | // |
377a7d60 | 1723 | TString sAlias(alias); |
1724 | sAlias.ReplaceAll("varname",varname); | |
1725 | TObjArray* oaAlias = TString(sAlias.Data()).Tokenize(":"); | |
fc03fa55 | 1726 | if (oaAlias->GetEntries()<2) { |
1727 | printf("Alias must have 2-6 arguments:\t%s\n", alias); | |
1728 | return 0; | |
1729 | } | |
377a7d60 | 1730 | char query[200]; |
1731 | TMultiGraph* multGr = new TMultiGraph(); | |
fc03fa55 | 1732 | Int_t marArr[6] = {24+igr%2, 20+igr%2, 20+igr%2, 20+igr%2, 20+igr%2, 20+igr%2}; |
1733 | Int_t colArr[6] = {kBlack, kBlack, kOrange, kRed, kGreen+1, kBlue}; | |
1734 | Double_t sizeArr[6] = {1.4, 1.1, 1.5, 1.1, 1.4, 0.8}; | |
1735 | Double_t shiftArr[6] = {0., 0., 0.25, 0.25, -0.25, -0.25}; | |
377a7d60 | 1736 | const Int_t ngr = oaAlias->GetEntriesFast(); |
1737 | for (Int_t i=0; i<ngr; i++){ | |
377a7d60 | 1738 | snprintf(query,200, "%f*(%s-0.5):%s", 1.+igr, oaAlias->At(i)->GetName(), var_x); |
fc03fa55 | 1739 | multGr->Add( (TGraphErrors*) TStatToolkit::MakeGraphSparse(tree,query,cut,marArr[i],colArr[i],sizeArr[i],shiftArr[i]) ); |
377a7d60 | 1740 | } |
377a7d60 | 1741 | // |
1742 | multGr->SetName(varname); | |
fc03fa55 | 1743 | multGr->SetTitle(varname); // used for y-axis labels of status bar, can be modified by calling function. |
d02b8756 | 1744 | delete oaVar; |
1745 | delete oaAlias; | |
377a7d60 | 1746 | return multGr; |
1747 | } | |
1748 | ||
1749 | ||
1750 | void TStatToolkit::AddStatusPad(TCanvas* c1, Float_t padratio, Float_t bottommargin) | |
1751 | { | |
1752 | // | |
1753 | // add pad to bottom of canvas for Status graphs (by Patrick Reichelt) | |
1754 | // call function "DrawStatusGraphs(...)" afterwards | |
1755 | // | |
1756 | TCanvas* c1_clone = (TCanvas*) c1->Clone("c1_clone"); | |
1757 | c1->Clear(); | |
1758 | // produce new pads | |
1759 | c1->cd(); | |
1760 | TPad* pad1 = new TPad("pad1", "pad1", 0., padratio, 1., 1.); | |
1761 | pad1->Draw(); | |
1762 | pad1->SetNumber(1); // so it can be called via "c1->cd(1);" | |
1763 | c1->cd(); | |
1764 | TPad* pad2 = new TPad("pad2", "pad2", 0., 0., 1., padratio); | |
1765 | pad2->Draw(); | |
1766 | pad2->SetNumber(2); | |
1767 | // draw original canvas into first pad | |
1768 | c1->cd(1); | |
1769 | c1_clone->DrawClonePad(); | |
1770 | pad1->SetBottomMargin(0.001); | |
1771 | pad1->SetRightMargin(0.01); | |
1772 | // set up second pad | |
1773 | c1->cd(2); | |
1774 | pad2->SetGrid(3); | |
1775 | pad2->SetTopMargin(0); | |
1776 | pad2->SetBottomMargin(bottommargin); // for the long x-axis labels (runnumbers) | |
1777 | pad2->SetRightMargin(0.01); | |
1778 | } | |
1779 | ||
1780 | ||
1781 | void TStatToolkit::DrawStatusGraphs(TObjArray* oaMultGr) | |
1782 | { | |
1783 | // | |
1784 | // draw Status graphs into active pad of canvas (by MI, Patrick Reichelt) | |
1785 | // ...into bottom pad, if called after "AddStatusPad(...)" | |
1786 | // | |
1787 | const Int_t nvars = oaMultGr->GetEntriesFast(); | |
1788 | TGraph* grAxis = (TGraph*) ((TMultiGraph*) oaMultGr->At(0))->GetListOfGraphs()->At(0); | |
1789 | grAxis->SetMaximum(0.5*nvars+0.5); | |
fc03fa55 | 1790 | grAxis->SetMinimum(0); |
377a7d60 | 1791 | grAxis->GetYaxis()->SetLabelSize(0); |
fc03fa55 | 1792 | grAxis->GetYaxis()->SetTitle(""); |
1793 | grAxis->SetTitle(""); | |
377a7d60 | 1794 | Int_t entries = grAxis->GetN(); |
377a7d60 | 1795 | grAxis->GetXaxis()->SetLabelSize(5.7*TMath::Min(TMath::Max(5./entries,0.01),0.03)); |
1796 | grAxis->GetXaxis()->LabelsOption("v"); | |
1797 | grAxis->Draw("ap"); | |
1798 | // | |
1799 | // draw multigraphs & names of status variables on the y axis | |
1800 | for (Int_t i=0; i<nvars; i++){ | |
1801 | ((TMultiGraph*) oaMultGr->At(i))->Draw("p"); | |
1802 | TLatex* ylabel = new TLatex(-0.1, 0.5*i+0.5, ((TMultiGraph*) oaMultGr->At(i))->GetTitle()); | |
1803 | ylabel->SetTextAlign(32); //hor:right & vert:centered | |
1804 | ylabel->SetTextSize(0.025/gPad->GetHNDC()); | |
1805 | ylabel->Draw(); | |
1806 | } | |
1807 | } | |
376089b6 | 1808 | |
1809 | ||
fc03fa55 | 1810 | TTree* TStatToolkit::WriteStatusToTree(TObject* oStatusGr) |
1811 | { | |
1812 | // | |
1813 | // Create Tree with Integers for each status variable flag (warning, outlier, physacc). | |
1814 | // (by Patrick Reichelt) | |
1815 | // | |
1816 | // input: either a TMultiGraph with status of single variable, which | |
1817 | // was computed by TStatToolkit::MakeStatusMultGr(), | |
1818 | // or a TObjArray which contains up to 10 of such variables. | |
92f227d3 | 1819 | // example: TTree* statusTree = WriteStatusToTree( TStatToolkit::MakeStatusMultGr(tree, "tpcItsMatch:run", "", sCriteria.Data(), 0) ); |
1820 | // or : TTree* statusTree = TStatToolkit::WriteStatusToTree(oaMultGr); | |
fc03fa55 | 1821 | // |
1822 | // output tree: 1=flag is true, 0=flag is false, -1=flag was not computed. | |
8e63b06b | 1823 | // To be rewritten to the pcstream |
fc03fa55 | 1824 | |
1825 | TObjArray* oaMultGr = NULL; | |
1826 | Bool_t needDeletion=kFALSE; | |
1827 | if (oStatusGr->IsA() == TObjArray::Class()) { | |
1828 | oaMultGr = (TObjArray*) oStatusGr; | |
1829 | } | |
1830 | else if (oStatusGr->IsA() == TMultiGraph::Class()) { | |
1831 | oaMultGr = new TObjArray(); needDeletion=kTRUE; | |
1832 | oaMultGr->Add((TMultiGraph*) oStatusGr); | |
1833 | } | |
1834 | else { | |
1835 | Printf("WriteStatusToTree(): Error! 'oStatusGr' must be a TMultiGraph or a TObjArray of them!"); | |
1836 | return 0; | |
1837 | } | |
1838 | // variables for output tree | |
1839 | const int nvarsMax=10; | |
1840 | const int ncritMax=5; | |
1841 | Int_t currentRun; | |
1842 | Int_t treevars[nvarsMax*ncritMax]; | |
1843 | TString varnames[nvarsMax*ncritMax]; | |
1844 | for (int i=0; i<nvarsMax*ncritMax; i++) treevars[i]=-1; | |
1845 | ||
1846 | Printf("WriteStatusToTree(): writing following variables to TTree (maybe only subset of listed criteria filled)"); | |
1847 | for (Int_t vari=0; vari<nvarsMax; vari++) | |
1848 | { | |
1849 | if (vari < oaMultGr->GetEntriesFast()) { | |
1850 | varnames[vari*ncritMax+0] = Form("%s_statisticOK", ((TMultiGraph*) oaMultGr->At(vari))->GetName()); | |
1851 | varnames[vari*ncritMax+1] = Form("%s_Warning", ((TMultiGraph*) oaMultGr->At(vari))->GetName()); | |
1852 | varnames[vari*ncritMax+2] = Form("%s_Outlier", ((TMultiGraph*) oaMultGr->At(vari))->GetName()); | |
1853 | varnames[vari*ncritMax+3] = Form("%s_PhysAcc", ((TMultiGraph*) oaMultGr->At(vari))->GetName()); | |
1854 | varnames[vari*ncritMax+4] = Form("%s_Extra", ((TMultiGraph*) oaMultGr->At(vari))->GetName()); | |
1855 | } | |
1856 | else { | |
1857 | varnames[vari*ncritMax+0] = Form("dummy"); | |
1858 | varnames[vari*ncritMax+1] = Form("dummy"); | |
1859 | varnames[vari*ncritMax+2] = Form("dummy"); | |
1860 | varnames[vari*ncritMax+3] = Form("dummy"); | |
1861 | varnames[vari*ncritMax+4] = Form("dummy"); | |
1862 | } | |
1863 | cout << " " << varnames[vari*ncritMax+0].Data() << " " << varnames[vari*ncritMax+1].Data() << " " << varnames[vari*ncritMax+2].Data() << " " << varnames[vari*ncritMax+3].Data() << " " << varnames[vari*ncritMax+4].Data() << endl; | |
1864 | } | |
1865 | ||
1866 | TTree* statusTree = new TTree("statusTree","statusTree"); | |
1867 | statusTree->Branch("run", ¤tRun ); | |
1868 | statusTree->Branch(varnames[ 0].Data(), &treevars[ 0]); | |
1869 | statusTree->Branch(varnames[ 1].Data(), &treevars[ 1]); | |
1870 | statusTree->Branch(varnames[ 2].Data(), &treevars[ 2]); | |
1871 | statusTree->Branch(varnames[ 3].Data(), &treevars[ 3]); | |
1872 | statusTree->Branch(varnames[ 4].Data(), &treevars[ 4]); | |
1873 | statusTree->Branch(varnames[ 5].Data(), &treevars[ 5]); | |
1874 | statusTree->Branch(varnames[ 6].Data(), &treevars[ 6]); | |
1875 | statusTree->Branch(varnames[ 7].Data(), &treevars[ 7]); | |
1876 | statusTree->Branch(varnames[ 8].Data(), &treevars[ 8]); | |
1877 | statusTree->Branch(varnames[ 9].Data(), &treevars[ 9]); | |
1878 | statusTree->Branch(varnames[10].Data(), &treevars[10]); | |
1879 | statusTree->Branch(varnames[11].Data(), &treevars[11]); | |
1880 | statusTree->Branch(varnames[12].Data(), &treevars[12]); | |
1881 | statusTree->Branch(varnames[13].Data(), &treevars[13]); | |
1882 | statusTree->Branch(varnames[14].Data(), &treevars[14]); | |
1883 | statusTree->Branch(varnames[15].Data(), &treevars[15]); | |
1884 | statusTree->Branch(varnames[16].Data(), &treevars[16]); | |
1885 | statusTree->Branch(varnames[17].Data(), &treevars[17]); | |
1886 | statusTree->Branch(varnames[18].Data(), &treevars[18]); | |
1887 | statusTree->Branch(varnames[19].Data(), &treevars[19]); | |
1888 | statusTree->Branch(varnames[20].Data(), &treevars[20]); | |
1889 | statusTree->Branch(varnames[21].Data(), &treevars[21]); | |
1890 | statusTree->Branch(varnames[22].Data(), &treevars[22]); | |
1891 | statusTree->Branch(varnames[23].Data(), &treevars[23]); | |
1892 | statusTree->Branch(varnames[24].Data(), &treevars[24]); | |
1893 | statusTree->Branch(varnames[25].Data(), &treevars[25]); | |
1894 | statusTree->Branch(varnames[26].Data(), &treevars[26]); | |
1895 | statusTree->Branch(varnames[27].Data(), &treevars[27]); | |
1896 | statusTree->Branch(varnames[28].Data(), &treevars[28]); | |
1897 | statusTree->Branch(varnames[29].Data(), &treevars[29]); | |
1898 | statusTree->Branch(varnames[30].Data(), &treevars[30]); | |
1899 | statusTree->Branch(varnames[31].Data(), &treevars[31]); | |
1900 | statusTree->Branch(varnames[32].Data(), &treevars[32]); | |
1901 | statusTree->Branch(varnames[33].Data(), &treevars[33]); | |
1902 | statusTree->Branch(varnames[34].Data(), &treevars[34]); | |
1903 | statusTree->Branch(varnames[35].Data(), &treevars[35]); | |
1904 | statusTree->Branch(varnames[36].Data(), &treevars[36]); | |
1905 | statusTree->Branch(varnames[37].Data(), &treevars[37]); | |
1906 | statusTree->Branch(varnames[38].Data(), &treevars[38]); | |
1907 | statusTree->Branch(varnames[39].Data(), &treevars[39]); | |
1908 | statusTree->Branch(varnames[40].Data(), &treevars[40]); | |
1909 | statusTree->Branch(varnames[41].Data(), &treevars[41]); | |
1910 | statusTree->Branch(varnames[42].Data(), &treevars[42]); | |
1911 | statusTree->Branch(varnames[43].Data(), &treevars[43]); | |
1912 | statusTree->Branch(varnames[44].Data(), &treevars[44]); | |
1913 | statusTree->Branch(varnames[45].Data(), &treevars[45]); | |
1914 | statusTree->Branch(varnames[46].Data(), &treevars[46]); | |
1915 | statusTree->Branch(varnames[47].Data(), &treevars[47]); | |
1916 | statusTree->Branch(varnames[48].Data(), &treevars[48]); | |
1917 | statusTree->Branch(varnames[49].Data(), &treevars[49]); | |
1918 | ||
1919 | // run loop | |
1920 | Double_t graphX; // x-position of marker (0.5, 1.5, ...) | |
1921 | Double_t graphY; // if >0 -> warning/outlier/physacc! if =-0.5 -> no warning/outlier/physacc | |
1922 | TList* arrRuns = (TList*) ((TGraph*) ((TMultiGraph*) oaMultGr->At(0))->GetListOfGraphs()->At(0))->GetXaxis()->GetLabels(); | |
1923 | //'TAxis->GetLabels()' returns THashList of TObjString, but using THashList gives compilation error "... incomplete type 'struct THashList' " | |
1924 | for (Int_t runi=0; runi<arrRuns->GetSize(); runi++) | |
1925 | { | |
1926 | currentRun = atoi( arrRuns->At(runi)->GetName() ); | |
1927 | //Printf(" runi=%2i, name: %s \t run number: %i", runi, arrRuns->At(runi)->GetName(), currentRun); | |
1928 | ||
1929 | // status variable loop | |
1930 | for (Int_t vari=0; vari<oaMultGr->GetEntriesFast(); vari++) | |
1931 | { | |
1932 | TMultiGraph* multGr = (TMultiGraph*) oaMultGr->At(vari); | |
1933 | ||
1934 | // criteria loop | |
1935 | // the order is given by TStatToolkit::MakeStatusMultGr(). | |
1936 | // criterion #1 is 'statisticOK' and mandatory, the rest is optional. (#0 is always True, thus skipped) | |
1937 | for (Int_t criti=1; criti<multGr->GetListOfGraphs()->GetEntries(); criti++) | |
1938 | { | |
1939 | TGraph* grCriterion = (TGraph*) multGr->GetListOfGraphs()->At(criti); | |
1940 | graphX = -1, graphY = -1; | |
1941 | grCriterion->GetPoint(runi, graphX, graphY); | |
1942 | treevars[(vari)*ncritMax+(criti-1)] = (graphY>0)?1:0; | |
1943 | } | |
1944 | } | |
1945 | statusTree->Fill(); | |
1946 | } | |
1947 | ||
1948 | if (needDeletion) delete oaMultGr; | |
1949 | ||
1950 | return statusTree; | |
1951 | } | |
1952 | ||
1953 | ||
8e63b06b | 1954 | void TStatToolkit::MakeSummaryTree(TTree* treeIn, TTreeSRedirector *pcstream, TObjString & sumID, TCut &selection){ |
1955 | // | |
1956 | // Make a summary tree for the input tree | |
1957 | // For the moment statistic works only for the primitive branches (Float/Double/Int) | |
1958 | // Extension recursive version planned for graphs a and histograms | |
1959 | // | |
1960 | // Following statistics are exctracted: | |
1961 | // - Standard: mean, meadian, rms | |
1962 | // - LTM robust statistic: mean60, rms60, mean90, rms90 | |
1963 | // Parameters: | |
1964 | // treeIn - input tree | |
1965 | // pctream - Output redirector | |
1966 | // sumID - ID as will be used in output tree | |
1967 | // selection - selection criteria define the set of entries used to evaluat statistic | |
1968 | // | |
1969 | TObjArray * brArray = treeIn->GetListOfBranches(); | |
1970 | Int_t tEntries= treeIn->GetEntries(); | |
1971 | Int_t nBranches=brArray->GetEntries(); | |
1972 | TString treeName = treeIn->GetName(); | |
1973 | treeName+="Summary"; | |
1974 | ||
1975 | (*pcstream)<<treeName.Data()<<"entries="<<tEntries; | |
1976 | (*pcstream)<<treeName.Data()<<"ID.="<<&sumID; | |
1977 | ||
1978 | TMatrixD valBranch(nBranches,7); | |
1979 | for (Int_t iBr=0; iBr<nBranches; iBr++){ | |
1980 | TString brName= brArray->At(iBr)->GetName(); | |
1981 | Int_t entries=treeIn->Draw(brArray->At(iBr)->GetName(),selection); | |
1982 | if (entries==0) continue; | |
1983 | Double_t median, mean, rms, mean60,rms60, mean90, rms90; | |
1984 | mean = TMath::Mean(entries,treeIn->GetV1()); | |
1985 | median= TMath::Median(entries,treeIn->GetV1()); | |
1986 | rms = TMath::RMS(entries,treeIn->GetV1()); | |
1987 | TStatToolkit::EvaluateUni(entries, treeIn->GetV1(), mean60,rms60,TMath::Min(TMath::Max(2., 0.60*entries),Double_t(entries))); | |
1988 | TStatToolkit::EvaluateUni(entries, treeIn->GetV1(), mean90,rms90,TMath::Min(TMath::Max(2., 0.90*entries),Double_t(entries))); | |
1989 | valBranch(iBr,0)=mean; | |
1990 | valBranch(iBr,1)=median; | |
1991 | valBranch(iBr,2)=rms; | |
1992 | valBranch(iBr,3)=mean60; | |
1993 | valBranch(iBr,4)=rms60; | |
1994 | valBranch(iBr,5)=mean90; | |
1995 | valBranch(iBr,6)=rms90; | |
1996 | (*pcstream)<<treeName.Data()<< | |
1997 | brName+"_Mean="<<valBranch(iBr,0)<< | |
1998 | brName+"_Median="<<valBranch(iBr,1)<< | |
1999 | brName+"_RMS="<<valBranch(iBr,2)<< | |
2000 | brName+"_Mean60="<<valBranch(iBr,3)<< | |
2001 | brName+"_RMS60="<<valBranch(iBr,4)<< | |
2002 | brName+"_Mean90="<<valBranch(iBr,5)<< | |
2003 | brName+"_RMS90="<<valBranch(iBr,6); | |
2004 | } | |
2005 | (*pcstream)<<treeName.Data()<<"\n"; | |
2006 | } | |
2007 | ||
2008 | ||
2009 | ||
fc03fa55 | 2010 | TMultiGraph* TStatToolkit::MakeStatusLines(TTree * tree, const char * expr, const char * cut, const char * alias) |
2011 | { | |
2012 | // | |
2013 | // Create status lines for trending using MakeGraphSparse(), very similar to MakeStatusMultGr(). | |
2014 | // (by Patrick Reichelt) | |
2015 | // | |
2016 | // format of expr : varname:xaxis (e.g. meanTPCncl:run, but 'varname' can be any string that you need for seach-and-replace) | |
2017 | // format of cut : char like in TCut | |
2018 | // format of alias: varname_OutlierMin:varname_OutlierMax:varname_WarningMin:varname_WarningMax:varname_PhysAccMin:varname_PhysAccMax:varname_RobustMean | |
2019 | // | |
2020 | TObjArray* oaVar = TString(expr).Tokenize(":"); | |
2021 | if (oaVar->GetEntries()<2) { | |
2022 | printf("Expression has to be of type 'varname:xaxis':\t%s\n", expr); | |
2023 | return 0; | |
2024 | } | |
2025 | char varname[50]; | |
2026 | char var_x[50]; | |
2027 | snprintf(varname,50,"%s", oaVar->At(0)->GetName()); | |
2028 | snprintf(var_x ,50,"%s", oaVar->At(1)->GetName()); | |
2029 | // | |
2030 | TString sAlias(alias); | |
2031 | if (sAlias.IsNull()) { // alias for default usage set here: | |
2032 | sAlias = "varname_OutlierMin:varname_OutlierMax:varname_WarningMin:varname_WarningMax:varname_PhysAccMin:varname_PhysAccMax:varname_RobustMean"; | |
2033 | } | |
2034 | sAlias.ReplaceAll("varname",varname); | |
2035 | TObjArray* oaAlias = TString(sAlias.Data()).Tokenize(":"); | |
2036 | if (oaAlias->GetEntries()<2) { | |
2037 | printf("Alias must have 2-7 arguments:\t%s\n", alias); | |
2038 | return 0; | |
2039 | } | |
2040 | char query[200]; | |
2041 | TMultiGraph* multGr = new TMultiGraph(); | |
2042 | Int_t colArr[7] = {kRed, kRed, kOrange, kOrange, kGreen+1, kGreen+1, kGray+2}; | |
2043 | const Int_t ngr = oaAlias->GetEntriesFast(); | |
2044 | for (Int_t i=0; i<ngr; i++){ | |
2045 | snprintf(query,200, "%s:%s", oaAlias->At(i)->GetName(), var_x); | |
2046 | multGr->Add( (TGraphErrors*) TStatToolkit::MakeGraphSparse(tree,query,cut,29,colArr[i],1.5) ); | |
2047 | } | |
2048 | // | |
2049 | multGr->SetName(varname); | |
2050 | multGr->SetTitle(varname); | |
2051 | delete oaVar; | |
2052 | delete oaAlias; | |
2053 | return multGr; | |
2054 | } | |
2055 | ||
2056 | ||
5d6816d1 | 2057 | TH1* TStatToolkit::DrawHistogram(TTree * tree, const char* drawCommand, const char* cuts, const char* histoname, const char* histotitle, Int_t nsigma, Float_t fraction ) |
376089b6 | 2058 | { |
2059 | // | |
2060 | // Draw histogram from TTree with robust range | |
2061 | // Only for 1D so far! | |
2062 | // | |
2063 | // Parameters: | |
2064 | // - histoname: name of histogram | |
2065 | // - histotitle: title of histgram | |
2066 | // - fraction: fraction of data to define the robust mean | |
2067 | // - nsigma: nsigma value for range | |
2068 | // | |
2069 | ||
2070 | TString drawStr(drawCommand); | |
2071 | TString cutStr(cuts); | |
2072 | Int_t dim = 1; | |
2073 | ||
376089b6 | 2074 | if(!tree) { |
2075 | cerr<<" Tree pointer is NULL!"<<endl; | |
5d6816d1 | 2076 | return 0; |
376089b6 | 2077 | } |
2078 | ||
3240a856 | 2079 | // get entries |
376089b6 | 2080 | Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff"); |
2081 | if (entries == -1) { | |
2082 | cerr<<"TTree draw returns -1"<<endl; | |
5d6816d1 | 2083 | return 0; |
376089b6 | 2084 | } |
2085 | ||
3240a856 | 2086 | // get dimension |
2087 | if(tree->GetV1()) dim = 1; | |
2088 | if(tree->GetV2()) dim = 2; | |
2089 | if(tree->GetV3()) dim = 3; | |
2090 | if(dim > 2){ | |
2091 | cerr<<"TTree has more than 2 dimensions (not yet supported)"<<endl; | |
5d6816d1 | 2092 | return 0; |
3240a856 | 2093 | } |
2094 | ||
2095 | // draw robust | |
2096 | Double_t meanX, rmsX=0; | |
2097 | Double_t meanY, rmsY=0; | |
2098 | TStatToolkit::EvaluateUni(entries, tree->GetV1(),meanX,rmsX, fraction*entries); | |
2099 | if(dim==2){ | |
2100 | TStatToolkit::EvaluateUni(entries, tree->GetV1(),meanY,rmsY, fraction*entries); | |
2101 | TStatToolkit::EvaluateUni(entries, tree->GetV2(),meanX,rmsX, fraction*entries); | |
2102 | } | |
fc03fa55 | 2103 | TH1* hOut=NULL; |
3240a856 | 2104 | if(dim==1){ |
2105 | hOut = new TH1F(histoname, histotitle, 200, meanX-nsigma*rmsX, meanX+nsigma*rmsX); | |
2106 | for (Int_t i=0; i<entries; i++) hOut->Fill(tree->GetV1()[i]); | |
2107 | hOut->GetXaxis()->SetTitle(tree->GetHistogram()->GetXaxis()->GetTitle()); | |
2108 | hOut->Draw(); | |
2109 | } | |
2110 | else if(dim==2){ | |
2111 | hOut = new TH2F(histoname, histotitle, 200, meanX-nsigma*rmsX, meanX+nsigma*rmsX,200, meanY-nsigma*rmsY, meanY+nsigma*rmsY); | |
2112 | for (Int_t i=0; i<entries; i++) hOut->Fill(tree->GetV2()[i],tree->GetV1()[i]); | |
2113 | hOut->GetXaxis()->SetTitle(tree->GetHistogram()->GetXaxis()->GetTitle()); | |
2114 | hOut->GetYaxis()->SetTitle(tree->GetHistogram()->GetYaxis()->GetTitle()); | |
2115 | hOut->Draw("colz"); | |
2116 | } | |
5d6816d1 | 2117 | return hOut; |
376089b6 | 2118 | } |
5d1391ec | 2119 | |
2120 | void TStatToolkit::CheckTreeAliases(TTree * tree, Int_t ncheck){ | |
2121 | // | |
2122 | // Check consistency of tree aliases | |
2123 | // | |
2124 | Int_t nCheck=100; | |
2125 | TList * aliases = (TList*)tree->GetListOfAliases(); | |
2126 | Int_t entries = aliases->GetEntries(); | |
2127 | for (Int_t i=0; i<entries; i++){ | |
2128 | TObject * object= aliases->At(i); | |
2129 | if (!object) continue; | |
2130 | Int_t ndraw=tree->Draw(aliases->At(i)->GetName(),"1","goff",nCheck); | |
2131 | if (ndraw==0){ | |
2132 | ::Error("Alias:\tProblem",aliases->At(i)->GetName()); | |
2133 | }else{ | |
2134 | ::Info("Alias:\tOK",aliases->At(i)->GetName()); | |
2135 | } | |
2136 | } | |
2137 | } |