]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STAT/TStatToolkit.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
PAR: pass proper defines for certain ROOT features
[u/mrichter/AliRoot.git] / STAT / TStatToolkit.cxx
CommitLineData
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
52using std::cout;
53using std::cerr;
54using std::endl;
55
21f3a443 56
57ClassImp(TStatToolkit) // Class implementation to enable ROOT I/O
58
59TStatToolkit::TStatToolkit() : TObject()
60{
61 //
62 // Default constructor
63 //
64}
65///////////////////////////////////////////////////////////////////////////
66TStatToolkit::~TStatToolkit()
67{
68 //
69 // Destructor
70 //
71}
72
73
74//_____________________________________________________________________________
75void 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
133void 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//_____________________________________________________________________________
189Int_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 235void 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 265void 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
308void 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
331Bool_t TStatToolkit::LTMHisto(TH1 *his1D, TVectorD &params , 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 450Double_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 546Double_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 686Float_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
734void 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
817TGraph2D * 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 860TGraphErrors * 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 938TString* 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
1040TString* 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
1150TString* 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 1253Int_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 1277TString TStatToolkit::FilterFit(const TString &input, const TString filter, TVectorD &param, 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
1303void 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 1356void TStatToolkit::Constrain1D(const TString &input, const TString filter, TVectorD &param, 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 1391TString TStatToolkit::MakeFitString(const TString &input, const TVectorD &param, 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 1409TGraphErrors * 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 1449TGraph * 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
1548Int_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
1601Int_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
1678TMultiGraph* 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
1750void 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
1781void 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 1810TTree* 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", &currentRun );
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 1954void 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 2010TMultiGraph* 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 2057TH1* 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
2120void 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}