]> 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
fine tuning of TOF tail (developing task)
[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 //
1697 TObjArray* oaVar = TString(expr).Tokenize(":");
fc03fa55 1698 if (oaVar->GetEntries()<2) {
1699 printf("Expression has to be of type 'varname:xaxis':\t%s\n", expr);
1700 return 0;
1701 }
377a7d60 1702 char varname[50];
1703 char var_x[50];
1704 snprintf(varname,50,"%s", oaVar->At(0)->GetName());
1705 snprintf(var_x ,50,"%s", oaVar->At(1)->GetName());
d02b8756 1706 //
377a7d60 1707 TString sAlias(alias);
1708 sAlias.ReplaceAll("varname",varname);
1709 TObjArray* oaAlias = TString(sAlias.Data()).Tokenize(":");
fc03fa55 1710 if (oaAlias->GetEntries()<2) {
1711 printf("Alias must have 2-6 arguments:\t%s\n", alias);
1712 return 0;
1713 }
377a7d60 1714 char query[200];
1715 TMultiGraph* multGr = new TMultiGraph();
fc03fa55 1716 Int_t marArr[6] = {24+igr%2, 20+igr%2, 20+igr%2, 20+igr%2, 20+igr%2, 20+igr%2};
1717 Int_t colArr[6] = {kBlack, kBlack, kOrange, kRed, kGreen+1, kBlue};
1718 Double_t sizeArr[6] = {1.4, 1.1, 1.5, 1.1, 1.4, 0.8};
1719 Double_t shiftArr[6] = {0., 0., 0.25, 0.25, -0.25, -0.25};
377a7d60 1720 const Int_t ngr = oaAlias->GetEntriesFast();
1721 for (Int_t i=0; i<ngr; i++){
377a7d60 1722 snprintf(query,200, "%f*(%s-0.5):%s", 1.+igr, oaAlias->At(i)->GetName(), var_x);
fc03fa55 1723 multGr->Add( (TGraphErrors*) TStatToolkit::MakeGraphSparse(tree,query,cut,marArr[i],colArr[i],sizeArr[i],shiftArr[i]) );
377a7d60 1724 }
377a7d60 1725 //
1726 multGr->SetName(varname);
fc03fa55 1727 multGr->SetTitle(varname); // used for y-axis labels of status bar, can be modified by calling function.
d02b8756 1728 delete oaVar;
1729 delete oaAlias;
377a7d60 1730 return multGr;
1731}
1732
1733
1734void TStatToolkit::AddStatusPad(TCanvas* c1, Float_t padratio, Float_t bottommargin)
1735{
1736 //
1737 // add pad to bottom of canvas for Status graphs (by Patrick Reichelt)
1738 // call function "DrawStatusGraphs(...)" afterwards
1739 //
1740 TCanvas* c1_clone = (TCanvas*) c1->Clone("c1_clone");
1741 c1->Clear();
1742 // produce new pads
1743 c1->cd();
1744 TPad* pad1 = new TPad("pad1", "pad1", 0., padratio, 1., 1.);
1745 pad1->Draw();
1746 pad1->SetNumber(1); // so it can be called via "c1->cd(1);"
1747 c1->cd();
1748 TPad* pad2 = new TPad("pad2", "pad2", 0., 0., 1., padratio);
1749 pad2->Draw();
1750 pad2->SetNumber(2);
1751 // draw original canvas into first pad
1752 c1->cd(1);
1753 c1_clone->DrawClonePad();
1754 pad1->SetBottomMargin(0.001);
1755 pad1->SetRightMargin(0.01);
1756 // set up second pad
1757 c1->cd(2);
1758 pad2->SetGrid(3);
1759 pad2->SetTopMargin(0);
1760 pad2->SetBottomMargin(bottommargin); // for the long x-axis labels (runnumbers)
1761 pad2->SetRightMargin(0.01);
1762}
1763
1764
1765void TStatToolkit::DrawStatusGraphs(TObjArray* oaMultGr)
1766{
1767 //
1768 // draw Status graphs into active pad of canvas (by MI, Patrick Reichelt)
1769 // ...into bottom pad, if called after "AddStatusPad(...)"
1770 //
1771 const Int_t nvars = oaMultGr->GetEntriesFast();
1772 TGraph* grAxis = (TGraph*) ((TMultiGraph*) oaMultGr->At(0))->GetListOfGraphs()->At(0);
1773 grAxis->SetMaximum(0.5*nvars+0.5);
fc03fa55 1774 grAxis->SetMinimum(0);
377a7d60 1775 grAxis->GetYaxis()->SetLabelSize(0);
fc03fa55 1776 grAxis->GetYaxis()->SetTitle("");
1777 grAxis->SetTitle("");
377a7d60 1778 Int_t entries = grAxis->GetN();
377a7d60 1779 grAxis->GetXaxis()->SetLabelSize(5.7*TMath::Min(TMath::Max(5./entries,0.01),0.03));
1780 grAxis->GetXaxis()->LabelsOption("v");
1781 grAxis->Draw("ap");
1782 //
1783 // draw multigraphs & names of status variables on the y axis
1784 for (Int_t i=0; i<nvars; i++){
1785 ((TMultiGraph*) oaMultGr->At(i))->Draw("p");
1786 TLatex* ylabel = new TLatex(-0.1, 0.5*i+0.5, ((TMultiGraph*) oaMultGr->At(i))->GetTitle());
1787 ylabel->SetTextAlign(32); //hor:right & vert:centered
1788 ylabel->SetTextSize(0.025/gPad->GetHNDC());
1789 ylabel->Draw();
1790 }
1791}
376089b6 1792
1793
fc03fa55 1794TTree* TStatToolkit::WriteStatusToTree(TObject* oStatusGr)
1795{
1796 //
1797 // Create Tree with Integers for each status variable flag (warning, outlier, physacc).
1798 // (by Patrick Reichelt)
1799 //
1800 // input: either a TMultiGraph with status of single variable, which
1801 // was computed by TStatToolkit::MakeStatusMultGr(),
1802 // or a TObjArray which contains up to 10 of such variables.
1803 //
1804 // output tree: 1=flag is true, 0=flag is false, -1=flag was not computed.
1805 //
1806
1807 TObjArray* oaMultGr = NULL;
1808 Bool_t needDeletion=kFALSE;
1809 if (oStatusGr->IsA() == TObjArray::Class()) {
1810 oaMultGr = (TObjArray*) oStatusGr;
1811 }
1812 else if (oStatusGr->IsA() == TMultiGraph::Class()) {
1813 oaMultGr = new TObjArray(); needDeletion=kTRUE;
1814 oaMultGr->Add((TMultiGraph*) oStatusGr);
1815 }
1816 else {
1817 Printf("WriteStatusToTree(): Error! 'oStatusGr' must be a TMultiGraph or a TObjArray of them!");
1818 return 0;
1819 }
1820 // variables for output tree
1821 const int nvarsMax=10;
1822 const int ncritMax=5;
1823 Int_t currentRun;
1824 Int_t treevars[nvarsMax*ncritMax];
1825 TString varnames[nvarsMax*ncritMax];
1826 for (int i=0; i<nvarsMax*ncritMax; i++) treevars[i]=-1;
1827
1828 Printf("WriteStatusToTree(): writing following variables to TTree (maybe only subset of listed criteria filled)");
1829 for (Int_t vari=0; vari<nvarsMax; vari++)
1830 {
1831 if (vari < oaMultGr->GetEntriesFast()) {
1832 varnames[vari*ncritMax+0] = Form("%s_statisticOK", ((TMultiGraph*) oaMultGr->At(vari))->GetName());
1833 varnames[vari*ncritMax+1] = Form("%s_Warning", ((TMultiGraph*) oaMultGr->At(vari))->GetName());
1834 varnames[vari*ncritMax+2] = Form("%s_Outlier", ((TMultiGraph*) oaMultGr->At(vari))->GetName());
1835 varnames[vari*ncritMax+3] = Form("%s_PhysAcc", ((TMultiGraph*) oaMultGr->At(vari))->GetName());
1836 varnames[vari*ncritMax+4] = Form("%s_Extra", ((TMultiGraph*) oaMultGr->At(vari))->GetName());
1837 }
1838 else {
1839 varnames[vari*ncritMax+0] = Form("dummy");
1840 varnames[vari*ncritMax+1] = Form("dummy");
1841 varnames[vari*ncritMax+2] = Form("dummy");
1842 varnames[vari*ncritMax+3] = Form("dummy");
1843 varnames[vari*ncritMax+4] = Form("dummy");
1844 }
1845 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;
1846 }
1847
1848 TTree* statusTree = new TTree("statusTree","statusTree");
1849 statusTree->Branch("run", &currentRun );
1850 statusTree->Branch(varnames[ 0].Data(), &treevars[ 0]);
1851 statusTree->Branch(varnames[ 1].Data(), &treevars[ 1]);
1852 statusTree->Branch(varnames[ 2].Data(), &treevars[ 2]);
1853 statusTree->Branch(varnames[ 3].Data(), &treevars[ 3]);
1854 statusTree->Branch(varnames[ 4].Data(), &treevars[ 4]);
1855 statusTree->Branch(varnames[ 5].Data(), &treevars[ 5]);
1856 statusTree->Branch(varnames[ 6].Data(), &treevars[ 6]);
1857 statusTree->Branch(varnames[ 7].Data(), &treevars[ 7]);
1858 statusTree->Branch(varnames[ 8].Data(), &treevars[ 8]);
1859 statusTree->Branch(varnames[ 9].Data(), &treevars[ 9]);
1860 statusTree->Branch(varnames[10].Data(), &treevars[10]);
1861 statusTree->Branch(varnames[11].Data(), &treevars[11]);
1862 statusTree->Branch(varnames[12].Data(), &treevars[12]);
1863 statusTree->Branch(varnames[13].Data(), &treevars[13]);
1864 statusTree->Branch(varnames[14].Data(), &treevars[14]);
1865 statusTree->Branch(varnames[15].Data(), &treevars[15]);
1866 statusTree->Branch(varnames[16].Data(), &treevars[16]);
1867 statusTree->Branch(varnames[17].Data(), &treevars[17]);
1868 statusTree->Branch(varnames[18].Data(), &treevars[18]);
1869 statusTree->Branch(varnames[19].Data(), &treevars[19]);
1870 statusTree->Branch(varnames[20].Data(), &treevars[20]);
1871 statusTree->Branch(varnames[21].Data(), &treevars[21]);
1872 statusTree->Branch(varnames[22].Data(), &treevars[22]);
1873 statusTree->Branch(varnames[23].Data(), &treevars[23]);
1874 statusTree->Branch(varnames[24].Data(), &treevars[24]);
1875 statusTree->Branch(varnames[25].Data(), &treevars[25]);
1876 statusTree->Branch(varnames[26].Data(), &treevars[26]);
1877 statusTree->Branch(varnames[27].Data(), &treevars[27]);
1878 statusTree->Branch(varnames[28].Data(), &treevars[28]);
1879 statusTree->Branch(varnames[29].Data(), &treevars[29]);
1880 statusTree->Branch(varnames[30].Data(), &treevars[30]);
1881 statusTree->Branch(varnames[31].Data(), &treevars[31]);
1882 statusTree->Branch(varnames[32].Data(), &treevars[32]);
1883 statusTree->Branch(varnames[33].Data(), &treevars[33]);
1884 statusTree->Branch(varnames[34].Data(), &treevars[34]);
1885 statusTree->Branch(varnames[35].Data(), &treevars[35]);
1886 statusTree->Branch(varnames[36].Data(), &treevars[36]);
1887 statusTree->Branch(varnames[37].Data(), &treevars[37]);
1888 statusTree->Branch(varnames[38].Data(), &treevars[38]);
1889 statusTree->Branch(varnames[39].Data(), &treevars[39]);
1890 statusTree->Branch(varnames[40].Data(), &treevars[40]);
1891 statusTree->Branch(varnames[41].Data(), &treevars[41]);
1892 statusTree->Branch(varnames[42].Data(), &treevars[42]);
1893 statusTree->Branch(varnames[43].Data(), &treevars[43]);
1894 statusTree->Branch(varnames[44].Data(), &treevars[44]);
1895 statusTree->Branch(varnames[45].Data(), &treevars[45]);
1896 statusTree->Branch(varnames[46].Data(), &treevars[46]);
1897 statusTree->Branch(varnames[47].Data(), &treevars[47]);
1898 statusTree->Branch(varnames[48].Data(), &treevars[48]);
1899 statusTree->Branch(varnames[49].Data(), &treevars[49]);
1900
1901 // run loop
1902 Double_t graphX; // x-position of marker (0.5, 1.5, ...)
1903 Double_t graphY; // if >0 -> warning/outlier/physacc! if =-0.5 -> no warning/outlier/physacc
1904 TList* arrRuns = (TList*) ((TGraph*) ((TMultiGraph*) oaMultGr->At(0))->GetListOfGraphs()->At(0))->GetXaxis()->GetLabels();
1905 //'TAxis->GetLabels()' returns THashList of TObjString, but using THashList gives compilation error "... incomplete type 'struct THashList' "
1906 for (Int_t runi=0; runi<arrRuns->GetSize(); runi++)
1907 {
1908 currentRun = atoi( arrRuns->At(runi)->GetName() );
1909 //Printf(" runi=%2i, name: %s \t run number: %i", runi, arrRuns->At(runi)->GetName(), currentRun);
1910
1911 // status variable loop
1912 for (Int_t vari=0; vari<oaMultGr->GetEntriesFast(); vari++)
1913 {
1914 TMultiGraph* multGr = (TMultiGraph*) oaMultGr->At(vari);
1915
1916 // criteria loop
1917 // the order is given by TStatToolkit::MakeStatusMultGr().
1918 // criterion #1 is 'statisticOK' and mandatory, the rest is optional. (#0 is always True, thus skipped)
1919 for (Int_t criti=1; criti<multGr->GetListOfGraphs()->GetEntries(); criti++)
1920 {
1921 TGraph* grCriterion = (TGraph*) multGr->GetListOfGraphs()->At(criti);
1922 graphX = -1, graphY = -1;
1923 grCriterion->GetPoint(runi, graphX, graphY);
1924 treevars[(vari)*ncritMax+(criti-1)] = (graphY>0)?1:0;
1925 }
1926 }
1927 statusTree->Fill();
1928 }
1929
1930 if (needDeletion) delete oaMultGr;
1931
1932 return statusTree;
1933}
1934
1935
1936TMultiGraph* TStatToolkit::MakeStatusLines(TTree * tree, const char * expr, const char * cut, const char * alias)
1937{
1938 //
1939 // Create status lines for trending using MakeGraphSparse(), very similar to MakeStatusMultGr().
1940 // (by Patrick Reichelt)
1941 //
1942 // format of expr : varname:xaxis (e.g. meanTPCncl:run, but 'varname' can be any string that you need for seach-and-replace)
1943 // format of cut : char like in TCut
1944 // format of alias: varname_OutlierMin:varname_OutlierMax:varname_WarningMin:varname_WarningMax:varname_PhysAccMin:varname_PhysAccMax:varname_RobustMean
1945 //
1946 TObjArray* oaVar = TString(expr).Tokenize(":");
1947 if (oaVar->GetEntries()<2) {
1948 printf("Expression has to be of type 'varname:xaxis':\t%s\n", expr);
1949 return 0;
1950 }
1951 char varname[50];
1952 char var_x[50];
1953 snprintf(varname,50,"%s", oaVar->At(0)->GetName());
1954 snprintf(var_x ,50,"%s", oaVar->At(1)->GetName());
1955 //
1956 TString sAlias(alias);
1957 if (sAlias.IsNull()) { // alias for default usage set here:
1958 sAlias = "varname_OutlierMin:varname_OutlierMax:varname_WarningMin:varname_WarningMax:varname_PhysAccMin:varname_PhysAccMax:varname_RobustMean";
1959 }
1960 sAlias.ReplaceAll("varname",varname);
1961 TObjArray* oaAlias = TString(sAlias.Data()).Tokenize(":");
1962 if (oaAlias->GetEntries()<2) {
1963 printf("Alias must have 2-7 arguments:\t%s\n", alias);
1964 return 0;
1965 }
1966 char query[200];
1967 TMultiGraph* multGr = new TMultiGraph();
1968 Int_t colArr[7] = {kRed, kRed, kOrange, kOrange, kGreen+1, kGreen+1, kGray+2};
1969 const Int_t ngr = oaAlias->GetEntriesFast();
1970 for (Int_t i=0; i<ngr; i++){
1971 snprintf(query,200, "%s:%s", oaAlias->At(i)->GetName(), var_x);
1972 multGr->Add( (TGraphErrors*) TStatToolkit::MakeGraphSparse(tree,query,cut,29,colArr[i],1.5) );
1973 }
1974 //
1975 multGr->SetName(varname);
1976 multGr->SetTitle(varname);
1977 delete oaVar;
1978 delete oaAlias;
1979 return multGr;
1980}
1981
1982
5d6816d1 1983TH1* TStatToolkit::DrawHistogram(TTree * tree, const char* drawCommand, const char* cuts, const char* histoname, const char* histotitle, Int_t nsigma, Float_t fraction )
376089b6 1984{
1985 //
1986 // Draw histogram from TTree with robust range
1987 // Only for 1D so far!
1988 //
1989 // Parameters:
1990 // - histoname: name of histogram
1991 // - histotitle: title of histgram
1992 // - fraction: fraction of data to define the robust mean
1993 // - nsigma: nsigma value for range
1994 //
1995
1996 TString drawStr(drawCommand);
1997 TString cutStr(cuts);
1998 Int_t dim = 1;
1999
376089b6 2000 if(!tree) {
2001 cerr<<" Tree pointer is NULL!"<<endl;
5d6816d1 2002 return 0;
376089b6 2003 }
2004
3240a856 2005 // get entries
376089b6 2006 Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff");
2007 if (entries == -1) {
2008 cerr<<"TTree draw returns -1"<<endl;
5d6816d1 2009 return 0;
376089b6 2010 }
2011
3240a856 2012 // get dimension
2013 if(tree->GetV1()) dim = 1;
2014 if(tree->GetV2()) dim = 2;
2015 if(tree->GetV3()) dim = 3;
2016 if(dim > 2){
2017 cerr<<"TTree has more than 2 dimensions (not yet supported)"<<endl;
5d6816d1 2018 return 0;
3240a856 2019 }
2020
2021 // draw robust
2022 Double_t meanX, rmsX=0;
2023 Double_t meanY, rmsY=0;
2024 TStatToolkit::EvaluateUni(entries, tree->GetV1(),meanX,rmsX, fraction*entries);
2025 if(dim==2){
2026 TStatToolkit::EvaluateUni(entries, tree->GetV1(),meanY,rmsY, fraction*entries);
2027 TStatToolkit::EvaluateUni(entries, tree->GetV2(),meanX,rmsX, fraction*entries);
2028 }
fc03fa55 2029 TH1* hOut=NULL;
3240a856 2030 if(dim==1){
2031 hOut = new TH1F(histoname, histotitle, 200, meanX-nsigma*rmsX, meanX+nsigma*rmsX);
2032 for (Int_t i=0; i<entries; i++) hOut->Fill(tree->GetV1()[i]);
2033 hOut->GetXaxis()->SetTitle(tree->GetHistogram()->GetXaxis()->GetTitle());
2034 hOut->Draw();
2035 }
2036 else if(dim==2){
2037 hOut = new TH2F(histoname, histotitle, 200, meanX-nsigma*rmsX, meanX+nsigma*rmsX,200, meanY-nsigma*rmsY, meanY+nsigma*rmsY);
2038 for (Int_t i=0; i<entries; i++) hOut->Fill(tree->GetV2()[i],tree->GetV1()[i]);
2039 hOut->GetXaxis()->SetTitle(tree->GetHistogram()->GetXaxis()->GetTitle());
2040 hOut->GetYaxis()->SetTitle(tree->GetHistogram()->GetYaxis()->GetTitle());
2041 hOut->Draw("colz");
2042 }
5d6816d1 2043 return hOut;
376089b6 2044}
5d1391ec 2045
2046void TStatToolkit::CheckTreeAliases(TTree * tree, Int_t ncheck){
2047 //
2048 // Check consistency of tree aliases
2049 //
2050 Int_t nCheck=100;
2051 TList * aliases = (TList*)tree->GetListOfAliases();
2052 Int_t entries = aliases->GetEntries();
2053 for (Int_t i=0; i<entries; i++){
2054 TObject * object= aliases->At(i);
2055 if (!object) continue;
2056 Int_t ndraw=tree->Draw(aliases->At(i)->GetName(),"1","goff",nCheck);
2057 if (ndraw==0){
2058 ::Error("Alias:\tProblem",aliases->At(i)->GetName());
2059 }else{
2060 ::Info("Alias:\tOK",aliases->At(i)->GetName());
2061 }
2062 }
2063}