Line data Source code
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 : /// \file TStatToolkit.cxx
19 : /// \class TStatToolkit
20 : /// \brief Summary of statistics functions
21 : /// Subset of matheamtical functions not included in the TMath
22 : //
23 : //
24 : /////////////////////////////////////////////////////////////////////////
25 : #include "TStopwatch.h"
26 : #include "TStatToolkit.h"
27 : #include "TTreeFormula.h"
28 :
29 : using std::cout;
30 : using std::cerr;
31 : using std::endl;
32 :
33 : //_____________________________________________________________________________
34 : void TStatToolkit::EvaluateUni(Int_t nvectors, Double_t *data, Double_t &mean
35 : , Double_t &sigma, Int_t hh)
36 : {
37 : //
38 : // Robust estimator in 1D case MI version - (faster than ROOT version)
39 : //
40 : // For the univariate case
41 : // estimates of location and scatter are returned in mean and sigma parameters
42 : // the algorithm works on the same principle as in multivariate case -
43 : // it finds a subset of size hh with smallest sigma, and then returns mean and
44 : // sigma of this subset
45 : //
46 :
47 0 : if (hh==0)
48 0 : hh=(nvectors+2)/2;
49 0 : 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};
50 0 : Int_t *index=new Int_t[nvectors];
51 0 : TMath::Sort(nvectors, data, index, kFALSE);
52 :
53 0 : Int_t nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11);
54 0 : Double_t factor = faclts[TMath::Max(0,nquant-1)];
55 :
56 : Double_t sumx =0;
57 : Double_t sumx2 =0;
58 : Int_t bestindex = -1;
59 : Double_t bestmean = 0;
60 0 : Double_t bestsigma = (data[index[nvectors-1]]-data[index[0]]+1.); // maximal possible sigma
61 0 : bestsigma *=bestsigma;
62 :
63 0 : for (Int_t i=0; i<hh; i++){
64 0 : sumx += data[index[i]];
65 0 : sumx2 += data[index[i]]*data[index[i]];
66 : }
67 :
68 0 : Double_t norm = 1./Double_t(hh);
69 0 : Double_t norm2 = (hh-1)>0 ? 1./Double_t(hh-1):1;
70 0 : for (Int_t i=hh; i<nvectors; i++){
71 0 : Double_t cmean = sumx*norm;
72 0 : Double_t csigma = (sumx2 - hh*cmean*cmean)*norm2;
73 0 : if (csigma<bestsigma){
74 : bestmean = cmean;
75 : bestsigma = csigma;
76 : bestindex = i-hh;
77 0 : }
78 :
79 0 : sumx += data[index[i]]-data[index[i-hh]];
80 0 : sumx2 += data[index[i]]*data[index[i]]-data[index[i-hh]]*data[index[i-hh]];
81 : }
82 :
83 0 : Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma));
84 0 : mean = bestmean;
85 0 : sigma = bstd;
86 0 : delete [] index;
87 :
88 0 : }
89 :
90 :
91 :
92 : void TStatToolkit::EvaluateUniExternal(Int_t nvectors, Double_t *data, Double_t &mean, Double_t &sigma, Int_t hh, Float_t externalfactor)
93 : {
94 : // Modified version of ROOT robust EvaluateUni
95 : // robust estimator in 1D case MI version
96 : // added external factor to include precision of external measurement
97 : //
98 :
99 0 : if (hh==0)
100 0 : hh=(nvectors+2)/2;
101 0 : 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};
102 0 : Int_t *index=new Int_t[nvectors];
103 0 : TMath::Sort(nvectors, data, index, kFALSE);
104 : //
105 0 : Int_t nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11);
106 0 : Double_t factor = faclts[0];
107 0 : if (nquant>0){
108 : // fix proper normalization - Anja
109 0 : factor = faclts[nquant-1];
110 0 : }
111 :
112 : //
113 : //
114 : Double_t sumx =0;
115 : Double_t sumx2 =0;
116 : Int_t bestindex = -1;
117 : Double_t bestmean = 0;
118 : Double_t bestsigma = -1;
119 0 : for (Int_t i=0; i<hh; i++){
120 0 : sumx += data[index[i]];
121 0 : sumx2 += data[index[i]]*data[index[i]];
122 : }
123 : //
124 0 : Double_t kfactor = 2.*externalfactor - externalfactor*externalfactor;
125 0 : Double_t norm = 1./Double_t(hh);
126 0 : for (Int_t i=hh; i<nvectors; i++){
127 0 : Double_t cmean = sumx*norm;
128 0 : Double_t csigma = (sumx2*norm - cmean*cmean*kfactor);
129 0 : if (csigma<bestsigma || bestsigma<0){
130 : bestmean = cmean;
131 : bestsigma = csigma;
132 : bestindex = i-hh;
133 0 : }
134 : //
135 : //
136 0 : sumx += data[index[i]]-data[index[i-hh]];
137 0 : sumx2 += data[index[i]]*data[index[i]]-data[index[i-hh]]*data[index[i-hh]];
138 : }
139 :
140 0 : Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma));
141 0 : mean = bestmean;
142 0 : sigma = bstd;
143 0 : delete [] index;
144 0 : }
145 :
146 :
147 : //_____________________________________________________________________________
148 : Int_t TStatToolkit::Freq(Int_t n, const Int_t *inlist
149 : , Int_t *outlist, Bool_t down)
150 : {
151 : //
152 : // Sort eleements according occurancy
153 : // The size of output array has is 2*n
154 : //
155 :
156 0 : Int_t * sindexS = new Int_t[n]; // temp array for sorting
157 0 : Int_t * sindexF = new Int_t[2*n];
158 0 : for (Int_t i=0;i<n;i++) sindexS[i]=0;
159 0 : for (Int_t i=0;i<2*n;i++) sindexF[i]=0;
160 : //
161 0 : TMath::Sort(n,inlist, sindexS, down);
162 0 : Int_t last = inlist[sindexS[0]];
163 : Int_t val = last;
164 0 : sindexF[0] = 1;
165 0 : sindexF[0+n] = last;
166 : Int_t countPos = 0;
167 : //
168 : // find frequency
169 0 : for(Int_t i=1;i<n; i++){
170 0 : val = inlist[sindexS[i]];
171 0 : if (last == val) sindexF[countPos]++;
172 : else{
173 0 : countPos++;
174 0 : sindexF[countPos+n] = val;
175 0 : sindexF[countPos]++;
176 : last =val;
177 : }
178 : }
179 0 : if (last==val) countPos++;
180 : // sort according frequency
181 0 : TMath::Sort(countPos, sindexF, sindexS, kTRUE);
182 0 : for (Int_t i=0;i<countPos;i++){
183 0 : outlist[2*i ] = sindexF[sindexS[i]+n];
184 0 : outlist[2*i+1] = sindexF[sindexS[i]];
185 : }
186 0 : delete [] sindexS;
187 0 : delete [] sindexF;
188 :
189 0 : return countPos;
190 :
191 : }
192 :
193 :
194 :
195 : void TStatToolkit::MedianFilter(TH1 * his1D, Int_t nmedian){
196 : //
197 : // Algorithm to filter histogram
198 : // author: marian.ivanov@cern.ch
199 : // Details of algorithm:
200 : // http://en.wikipedia.org/w/index.php?title=Median_filter&oldid=582191524
201 : // Input parameters:
202 : // his1D - input histogam - to be modiefied by Medianfilter
203 : // nmendian - number of bins in median filter
204 : //
205 0 : Int_t nbins = his1D->GetNbinsX();
206 0 : TVectorD vectorH(nbins);
207 0 : for (Int_t ibin=0; ibin<nbins; ibin++) vectorH[ibin]=his1D->GetBinContent(ibin+1);
208 0 : for (Int_t ibin=0; ibin<nbins; ibin++) {
209 0 : Int_t index0=ibin-nmedian;
210 0 : Int_t index1=ibin+nmedian;
211 0 : if (index0<0) {index1+=-index0; index0=0;}
212 0 : if (index1>=nbins) {index0-=index1-nbins+1; index1=nbins-1;}
213 0 : Double_t value= TMath::Median(index1-index0,&(vectorH.GetMatrixArray()[index0]));
214 0 : his1D->SetBinContent(ibin+1, value);
215 : }
216 0 : }
217 :
218 :
219 :
220 : Float_t TStatToolkit::GetCOG(const Short_t *arr, Int_t nBins, Float_t xMin, Float_t xMax, Float_t *rms, Float_t *sum)
221 : {
222 : //
223 : // calculate center of gravity rms and sum for array 'arr' with nBins an a x range xMin to xMax
224 : // return COG; in case of failure return xMin
225 : //
226 : Float_t meanCOG = 0;
227 : Float_t rms2COG = 0;
228 : Float_t sumCOG = 0;
229 : Int_t npoints = 0;
230 :
231 0 : Float_t binWidth = (xMax-xMin)/(Float_t)nBins;
232 :
233 0 : for (Int_t ibin=0; ibin<nBins; ibin++){
234 0 : Float_t entriesI = (Float_t)arr[ibin];
235 0 : Double_t xcenter = xMin+(ibin+0.5)*binWidth;
236 0 : if ( entriesI>0 ){
237 0 : meanCOG += xcenter*entriesI;
238 0 : rms2COG += xcenter*entriesI*xcenter;
239 0 : sumCOG += entriesI;
240 0 : npoints++;
241 0 : }
242 : }
243 0 : if ( sumCOG == 0 ) return xMin;
244 0 : meanCOG/=sumCOG;
245 :
246 0 : if ( rms ){
247 0 : rms2COG /=sumCOG;
248 0 : (*rms) = TMath::Sqrt(TMath::Abs(meanCOG*meanCOG-rms2COG));
249 0 : if ( npoints == 1 ) (*rms) = binWidth/TMath::Sqrt(12);
250 : }
251 :
252 0 : if ( sum )
253 0 : (*sum) = sumCOG;
254 :
255 0 : return meanCOG;
256 0 : }
257 :
258 :
259 :
260 : ///////////////////////////////////////////////////////////////
261 : ////////////// TEST functions /////////////////////////
262 : ///////////////////////////////////////////////////////////////
263 :
264 :
265 :
266 :
267 :
268 : void TStatToolkit::TestGausFit(Int_t nhistos){
269 : //
270 : // Test performance of the parabolic - gaussian fit - compare it with
271 : // ROOT gauss fit
272 : // nhistos - number of histograms to be used for test
273 : //
274 0 : TTreeSRedirector *pcstream = new TTreeSRedirector("fitdebug.root","recreate");
275 :
276 0 : Float_t *xTrue = new Float_t[nhistos];
277 0 : Float_t *sTrue = new Float_t[nhistos];
278 0 : TVectorD **par1 = new TVectorD*[nhistos];
279 0 : TVectorD **par2 = new TVectorD*[nhistos];
280 0 : TMatrixD dummy(3,3);
281 :
282 :
283 0 : TH1F **h1f = new TH1F*[nhistos];
284 0 : TF1 *myg = new TF1("myg","gaus");
285 0 : TF1 *fit = new TF1("fit","gaus");
286 0 : gRandom->SetSeed(0);
287 :
288 : //init
289 0 : for (Int_t i=0;i<nhistos; i++){
290 0 : par1[i] = new TVectorD(3);
291 0 : par2[i] = new TVectorD(3);
292 0 : h1f[i] = new TH1F(Form("h1f%d",i),Form("h1f%d",i),20,-10,10);
293 0 : xTrue[i]= gRandom->Rndm();
294 0 : gSystem->Sleep(2);
295 0 : sTrue[i]= .75+gRandom->Rndm()*.5;
296 0 : myg->SetParameters(1,xTrue[i],sTrue[i]);
297 0 : h1f[i]->FillRandom("myg");
298 : }
299 :
300 0 : TStopwatch s;
301 0 : s.Start();
302 : //standard gaus fit
303 0 : for (Int_t i=0; i<nhistos; i++){
304 0 : h1f[i]->Fit(fit,"0q");
305 0 : (*par1[i])(0) = fit->GetParameter(0);
306 0 : (*par1[i])(1) = fit->GetParameter(1);
307 0 : (*par1[i])(2) = fit->GetParameter(2);
308 : }
309 0 : s.Stop();
310 0 : printf("Gaussian fit\t");
311 0 : s.Print();
312 :
313 0 : s.Start();
314 : //TStatToolkit gaus fit
315 0 : for (Int_t i=0; i<nhistos; i++){
316 0 : TStatToolkit::FitGaus(h1f[i]->GetArray()+1,h1f[i]->GetNbinsX(),h1f[i]->GetXaxis()->GetXmin(),h1f[i]->GetXaxis()->GetXmax(),par2[i],&dummy);
317 : }
318 :
319 0 : s.Stop();
320 0 : printf("Parabolic fit\t");
321 0 : s.Print();
322 : //write stream
323 0 : for (Int_t i=0;i<nhistos; i++){
324 0 : Float_t xt = xTrue[i];
325 0 : Float_t st = sTrue[i];
326 0 : (*pcstream)<<"data"
327 0 : <<"xTrue="<<xt
328 0 : <<"sTrue="<<st
329 0 : <<"pg.="<<(par1[i])
330 0 : <<"pa.="<<(par2[i])
331 0 : <<"\n";
332 0 : }
333 : //delete pointers
334 0 : for (Int_t i=0;i<nhistos; i++){
335 0 : delete par1[i];
336 0 : delete par2[i];
337 0 : delete h1f[i];
338 : }
339 0 : delete pcstream;
340 0 : delete []h1f;
341 0 : delete []xTrue;
342 0 : delete []sTrue;
343 : //
344 0 : delete []par1;
345 0 : delete []par2;
346 :
347 0 : }
348 :
349 :
350 :
351 : TGraph2D * TStatToolkit::MakeStat2D(TH3 * his, Int_t delta0, Int_t delta1, Int_t type){
352 : //
353 : //
354 : //
355 : // delta - number of bins to integrate
356 : // type - 0 - mean value
357 :
358 0 : TAxis * xaxis = his->GetXaxis();
359 0 : TAxis * yaxis = his->GetYaxis();
360 : // TAxis * zaxis = his->GetZaxis();
361 0 : Int_t nbinx = xaxis->GetNbins();
362 0 : Int_t nbiny = yaxis->GetNbins();
363 0 : char name[1000];
364 : Int_t icount=0;
365 0 : TGraph2D *graph = new TGraph2D(nbinx*nbiny);
366 0 : TF1 f1("f1","gaus");
367 0 : for (Int_t ix=0; ix<nbinx;ix++)
368 0 : for (Int_t iy=0; iy<nbiny;iy++){
369 0 : Float_t xcenter = xaxis->GetBinCenter(ix);
370 0 : Float_t ycenter = yaxis->GetBinCenter(iy);
371 0 : snprintf(name,1000,"%s_%d_%d",his->GetName(), ix,iy);
372 0 : TH1 *projection = his->ProjectionZ(name,ix-delta0,ix+delta0,iy-delta1,iy+delta1);
373 : Float_t stat= 0;
374 0 : if (type==0) stat = projection->GetMean();
375 0 : if (type==1) stat = projection->GetRMS();
376 0 : if (type==2 || type==3){
377 0 : TVectorD vec(10);
378 0 : TStatToolkit::LTM((TH1F*)projection,&vec,0.7);
379 0 : if (type==2) stat= vec[1];
380 0 : if (type==3) stat= vec[0];
381 0 : }
382 0 : if (type==4|| type==5){
383 0 : projection->Fit(&f1);
384 0 : if (type==4) stat= f1.GetParameter(1);
385 0 : if (type==5) stat= f1.GetParameter(2);
386 : }
387 : //printf("%d\t%f\t%f\t%f\n", icount,xcenter, ycenter, stat);
388 0 : graph->SetPoint(icount,xcenter, ycenter, stat);
389 0 : icount++;
390 : }
391 : return graph;
392 0 : }
393 :
394 : TGraphErrors * TStatToolkit::MakeStat1D(TH2 * his, Int_t deltaBin, Double_t fraction, Int_t returnType, Int_t markerStyle, Int_t markerColor){
395 : //
396 : // function to retrieve the "mean and RMS estimate" of 2D histograms
397 : //
398 : // Robust statistic to estimate properties of the distribution
399 : // See http://en.wikipedia.org/wiki/Trimmed_estimator
400 : //
401 : // deltaBin - number of bins to integrate (bin+-deltaBin)
402 : // fraction - fraction of values for the LTM and for the gauss fit
403 : // returnType -
404 : // 0 - mean value
405 : // 1 - RMS
406 : // 2 - LTM mean
407 : // 3 - LTM sigma
408 : // 4 - Gaus fit mean - on LTM range
409 : // 5 - Gaus fit sigma - on LTM range
410 : // 6 - Robust bin median
411 : //
412 0 : TAxis * xaxis = his->GetXaxis();
413 0 : Int_t nbinx = xaxis->GetNbins();
414 0 : char name[1000];
415 : Int_t icount=0;
416 : //
417 0 : TVectorD vecX(nbinx);
418 0 : TVectorD vecXErr(nbinx);
419 0 : TVectorD vecY(nbinx);
420 0 : TVectorD vecYErr(nbinx);
421 : //
422 0 : TF1 f1("f1","gaus");
423 0 : TVectorD vecLTM(10);
424 :
425 0 : for (Int_t jx=1; jx<=nbinx;jx++){
426 0 : Int_t ix=jx-1;
427 0 : Float_t xcenter = xaxis->GetBinCenter(jx);
428 0 : snprintf(name,1000,"%s_%d",his->GetName(), ix);
429 0 : TH1 *projection = his->ProjectionY(name,TMath::Max(jx-deltaBin,1),TMath::Min(jx+deltaBin,nbinx));
430 : Double_t stat= 0;
431 : Double_t err =0;
432 0 : TStatToolkit::LTMHisto((TH1F*)projection,vecLTM,fraction);
433 : //
434 0 : if (returnType==0) {
435 0 : stat = projection->GetMean();
436 0 : err = projection->GetMeanError();
437 0 : }
438 0 : else if (returnType==1) {
439 0 : stat = projection->GetRMS();
440 0 : err = projection->GetRMSError();
441 0 : }
442 0 : else if (returnType==2 || returnType==3){
443 0 : if (returnType==2) {stat= vecLTM[1]; err =projection->GetRMSError();}
444 0 : if (returnType==3) {stat= vecLTM[2]; err =projection->GetRMSError();}
445 : }
446 0 : else if (returnType==4|| returnType==5){
447 0 : f1.SetParameters(vecLTM[0], vecLTM[1], vecLTM[2]+0.05);
448 0 : projection->Fit(&f1,"QN","QN", vecLTM[7]-vecLTM[2], vecLTM[8]+vecLTM[2]);
449 0 : if (returnType==4) {
450 0 : stat= f1.GetParameter(1);
451 0 : err=f1.GetParError(1);
452 0 : }
453 0 : if (returnType==5) {
454 0 : stat= f1.GetParameter(2);
455 0 : err=f1.GetParError(2);
456 0 : }
457 : }
458 0 : else if (returnType==6) {
459 0 : stat=RobustBinMedian(projection,fraction);
460 0 : }
461 :
462 0 : vecX[icount]=xcenter;
463 0 : vecY[icount]=stat;
464 0 : vecYErr[icount]=err;
465 0 : icount++;
466 0 : delete projection;
467 : }
468 0 : TGraphErrors *graph = new TGraphErrors(icount,vecX.GetMatrixArray(), vecY.GetMatrixArray(),0, vecYErr.GetMatrixArray());
469 0 : graph->SetMarkerStyle(markerStyle);
470 0 : graph->SetMarkerColor(markerColor);
471 : return graph;
472 0 : }
473 :
474 : Double_t TStatToolkit::RobustBinMedian(TH1* hist, Double_t fractionCut/*=1.*/)
475 : {
476 : // Robust median with average from neighbouring bins
477 0 : const Int_t nbins1D=hist->GetNbinsX();
478 : Double_t binMedian=0;
479 0 : Double_t limits[2]={hist->GetBinCenter(1), hist->GetBinCenter(nbins1D)};
480 :
481 0 : Double_t* integral=hist->GetIntegral();
482 0 : for (Int_t i=1; i<nbins1D-1; i++){
483 0 : if (integral[i-1]<0.5 && integral[i]>=0.5){
484 0 : if (hist->GetBinContent(i-1)+hist->GetBinContent(i)>0){
485 0 : binMedian=hist->GetBinCenter(i);
486 0 : Double_t dIdx=-(integral[i-1]-integral[i]);
487 0 : Double_t dx=(0.5+(0.5-integral[i])/dIdx)*hist->GetBinWidth(i);
488 0 : binMedian+=dx;
489 0 : }
490 : }
491 0 : if (integral[i-1]<fractionCut && integral[i]>=fractionCut){
492 0 : limits[0]=hist->GetBinCenter(i-1)-hist->GetBinWidth(i);
493 0 : }
494 0 : if (integral[i]<1-fractionCut && integral[i+1]>=1-fractionCut){
495 0 : limits[1]=hist->GetBinCenter(i+1)+hist->GetBinWidth(i);
496 0 : }
497 : }
498 :
499 0 : return binMedian;
500 : }
501 :
502 :
503 : TString* TStatToolkit::FitPlane(TTree *tree, const char* drawCommand, const char* formula, const char* cuts, Double_t & chi2, Int_t &npoints, TVectorD &fitParam, TMatrixD &covMatrix, Float_t frac, Int_t start, Int_t stop,Bool_t fix0){
504 : //
505 : // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts
506 : // returns chi2, fitParam and covMatrix
507 : // returns TString with fitted formula
508 : //
509 :
510 0 : TString formulaStr(formula);
511 0 : TString drawStr(drawCommand);
512 0 : TString cutStr(cuts);
513 0 : TString ferr("1");
514 :
515 0 : TString strVal(drawCommand);
516 0 : if (strVal.Contains(":")){
517 0 : TObjArray* valTokens = strVal.Tokenize(":");
518 0 : drawStr = valTokens->At(0)->GetName();
519 0 : ferr = valTokens->At(1)->GetName();
520 0 : delete valTokens;
521 0 : }
522 :
523 :
524 0 : formulaStr.ReplaceAll("++", "~");
525 0 : TObjArray* formulaTokens = formulaStr.Tokenize("~");
526 0 : Int_t dim = formulaTokens->GetEntriesFast();
527 :
528 0 : fitParam.ResizeTo(dim);
529 0 : covMatrix.ResizeTo(dim,dim);
530 :
531 0 : TLinearFitter* fitter = new TLinearFitter(dim+1, Form("hyp%d",dim));
532 0 : fitter->StoreData(kTRUE);
533 0 : fitter->ClearPoints();
534 :
535 0 : Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start, start);
536 0 : if (entries == -1) {
537 0 : delete formulaTokens;
538 0 : return new TString(TString::Format("ERROR expr: %s\t%s\tEntries==0",drawStr.Data(),cutStr.Data()));
539 : }
540 0 : Double_t **values = new Double_t*[dim+1] ;
541 0 : for (Int_t i=0; i<dim+1; i++) values[i]=NULL;
542 : //
543 0 : entries = tree->Draw(ferr.Data(), cutStr.Data(), "goff", stop-start, start);
544 0 : if (entries == -1) {
545 0 : delete formulaTokens;
546 0 : delete []values;
547 0 : return new TString(TString::Format("ERROR error part: %s\t%s\tEntries==0",ferr.Data(),cutStr.Data()));
548 : }
549 0 : Double_t *errors = new Double_t[entries];
550 0 : memcpy(errors, tree->GetV1(), entries*sizeof(Double_t));
551 :
552 0 : for (Int_t i = 0; i < dim + 1; i++){
553 : Int_t centries = 0;
554 0 : if (i < dim) centries = tree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start);
555 0 : else centries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start);
556 :
557 0 : if (entries != centries) {
558 0 : delete []errors;
559 0 : delete []values;
560 0 : return new TString(TString::Format("ERROR: %s\t%s\tEntries==%d\tEntries2=%d\n",drawStr.Data(),cutStr.Data(),entries,centries));
561 : }
562 0 : values[i] = new Double_t[entries];
563 0 : memcpy(values[i], tree->GetV1(), entries*sizeof(Double_t));
564 0 : }
565 :
566 : // add points to the fitter
567 0 : for (Int_t i = 0; i < entries; i++){
568 0 : Double_t x[1000];
569 0 : for (Int_t j=0; j<dim;j++) x[j]=values[j][i];
570 0 : fitter->AddPoint(x, values[dim][i], errors[i]);
571 0 : }
572 :
573 0 : fitter->Eval();
574 0 : if (frac>0.5 && frac<1){
575 0 : fitter->EvalRobust(frac);
576 : }else{
577 0 : if (fix0) {
578 0 : fitter->FixParameter(0,0);
579 0 : fitter->Eval();
580 : }
581 : }
582 0 : fitter->GetParameters(fitParam);
583 0 : fitter->GetCovarianceMatrix(covMatrix);
584 0 : chi2 = fitter->GetChisquare();
585 0 : npoints = entries;
586 0 : TString *preturnFormula = new TString(Form("( %f+",fitParam[0])), &returnFormula = *preturnFormula;
587 :
588 0 : for (Int_t iparam = 0; iparam < dim; iparam++) {
589 0 : returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam+1]));
590 0 : if (iparam < dim-1) returnFormula.Append("+");
591 : }
592 0 : returnFormula.Append(" )");
593 :
594 :
595 0 : for (Int_t j=0; j<dim+1;j++) delete [] values[j];
596 :
597 :
598 0 : delete formulaTokens;
599 0 : delete fitter;
600 0 : delete[] values;
601 0 : delete[] errors;
602 : return preturnFormula;
603 0 : }
604 :
605 : TString* TStatToolkit::FitPlaneConstrain(TTree *tree, const char* drawCommand, const char* formula, const char* cuts, Double_t & chi2, Int_t &npoints, TVectorD &fitParam, TMatrixD &covMatrix, Float_t frac, Int_t start, Int_t stop,Double_t constrain){
606 : //
607 : // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts
608 : // returns chi2, fitParam and covMatrix
609 : // returns TString with fitted formula
610 : //
611 :
612 0 : TString formulaStr(formula);
613 0 : TString drawStr(drawCommand);
614 0 : TString cutStr(cuts);
615 0 : TString ferr("1");
616 :
617 0 : TString strVal(drawCommand);
618 0 : if (strVal.Contains(":")){
619 0 : TObjArray* valTokens = strVal.Tokenize(":");
620 0 : drawStr = valTokens->At(0)->GetName();
621 0 : ferr = valTokens->At(1)->GetName();
622 0 : delete valTokens;
623 0 : }
624 :
625 :
626 0 : formulaStr.ReplaceAll("++", "~");
627 0 : TObjArray* formulaTokens = formulaStr.Tokenize("~");
628 0 : Int_t dim = formulaTokens->GetEntriesFast();
629 :
630 0 : fitParam.ResizeTo(dim);
631 0 : covMatrix.ResizeTo(dim,dim);
632 :
633 0 : TLinearFitter* fitter = new TLinearFitter(dim+1, Form("hyp%d",dim));
634 0 : fitter->StoreData(kTRUE);
635 0 : fitter->ClearPoints();
636 :
637 0 : Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start, start);
638 0 : if (entries == -1) {
639 0 : delete formulaTokens;
640 0 : return new TString("An ERROR has occured during fitting!");
641 : }
642 0 : Double_t **values = new Double_t*[dim+1] ;
643 0 : for (Int_t i=0; i<dim+1; i++) values[i]=NULL;
644 : //
645 0 : entries = tree->Draw(ferr.Data(), cutStr.Data(), "goff", stop-start, start);
646 0 : if (entries == -1) {
647 0 : delete formulaTokens;
648 0 : delete [] values;
649 0 : return new TString("An ERROR has occured during fitting!");
650 : }
651 0 : Double_t *errors = new Double_t[entries];
652 0 : memcpy(errors, tree->GetV1(), entries*sizeof(Double_t));
653 :
654 0 : for (Int_t i = 0; i < dim + 1; i++){
655 : Int_t centries = 0;
656 0 : if (i < dim) centries = tree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start);
657 0 : else centries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start);
658 :
659 0 : if (entries != centries) {
660 0 : delete []errors;
661 0 : delete []values;
662 0 : delete formulaTokens;
663 0 : return new TString("An ERROR has occured during fitting!");
664 : }
665 0 : values[i] = new Double_t[entries];
666 0 : memcpy(values[i], tree->GetV1(), entries*sizeof(Double_t));
667 0 : }
668 :
669 : // add points to the fitter
670 0 : for (Int_t i = 0; i < entries; i++){
671 0 : Double_t x[1000];
672 0 : for (Int_t j=0; j<dim;j++) x[j]=values[j][i];
673 0 : fitter->AddPoint(x, values[dim][i], errors[i]);
674 0 : }
675 0 : if (constrain>0){
676 0 : for (Int_t i = 0; i < dim; i++){
677 0 : Double_t x[1000];
678 0 : for (Int_t j=0; j<dim;j++) if (i!=j) x[j]=0;
679 0 : x[i]=1.;
680 0 : fitter->AddPoint(x, 0, constrain);
681 0 : }
682 0 : }
683 :
684 :
685 0 : fitter->Eval();
686 0 : if (frac>0.5 && frac<1){
687 0 : fitter->EvalRobust(frac);
688 : }
689 0 : fitter->GetParameters(fitParam);
690 0 : fitter->GetCovarianceMatrix(covMatrix);
691 0 : chi2 = fitter->GetChisquare();
692 0 : npoints = entries;
693 :
694 0 : TString *preturnFormula = new TString(Form("( %f+",fitParam[0])), &returnFormula = *preturnFormula;
695 :
696 0 : for (Int_t iparam = 0; iparam < dim; iparam++) {
697 0 : returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam+1]));
698 0 : if (iparam < dim-1) returnFormula.Append("+");
699 : }
700 0 : returnFormula.Append(" )");
701 :
702 0 : for (Int_t j=0; j<dim+1;j++) delete [] values[j];
703 :
704 :
705 :
706 0 : delete formulaTokens;
707 0 : delete fitter;
708 0 : delete[] values;
709 0 : delete[] errors;
710 : return preturnFormula;
711 0 : }
712 :
713 :
714 :
715 : TString* TStatToolkit::FitPlaneFixed(TTree *tree, const char* drawCommand, const char* formula, const char* cuts, Double_t & chi2, Int_t &npoints, TVectorD &fitParam, TMatrixD &covMatrix, Float_t frac, Int_t start, Int_t stop){
716 : //
717 : // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts
718 : // returns chi2, fitParam and covMatrix
719 : // returns TString with fitted formula
720 : //
721 :
722 0 : TString formulaStr(formula);
723 0 : TString drawStr(drawCommand);
724 0 : TString cutStr(cuts);
725 0 : TString ferr("1");
726 :
727 0 : TString strVal(drawCommand);
728 0 : if (strVal.Contains(":")){
729 0 : TObjArray* valTokens = strVal.Tokenize(":");
730 0 : drawStr = valTokens->At(0)->GetName();
731 0 : ferr = valTokens->At(1)->GetName();
732 0 : delete valTokens;
733 0 : }
734 :
735 :
736 0 : formulaStr.ReplaceAll("++", "~");
737 0 : TObjArray* formulaTokens = formulaStr.Tokenize("~");
738 0 : Int_t dim = formulaTokens->GetEntriesFast();
739 :
740 0 : fitParam.ResizeTo(dim);
741 0 : covMatrix.ResizeTo(dim,dim);
742 0 : TString fitString="x0";
743 0 : for (Int_t i=1; i<dim; i++) fitString+=Form("++x%d",i);
744 0 : TLinearFitter* fitter = new TLinearFitter(dim, fitString.Data());
745 0 : fitter->StoreData(kTRUE);
746 0 : fitter->ClearPoints();
747 :
748 0 : Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start, start);
749 0 : if (entries == -1) {
750 0 : delete formulaTokens;
751 0 : return new TString("An ERROR has occured during fitting!");
752 : }
753 0 : Double_t **values = new Double_t*[dim+1] ;
754 0 : for (Int_t i=0; i<dim+1; i++) values[i]=NULL;
755 : //
756 0 : entries = tree->Draw(ferr.Data(), cutStr.Data(), "goff", stop-start, start);
757 0 : if (entries == -1) {
758 0 : delete []values;
759 0 : delete formulaTokens;
760 0 : return new TString("An ERROR has occured during fitting!");
761 : }
762 0 : Double_t *errors = new Double_t[entries];
763 0 : memcpy(errors, tree->GetV1(), entries*sizeof(Double_t));
764 :
765 0 : for (Int_t i = 0; i < dim + 1; i++){
766 : Int_t centries = 0;
767 0 : if (i < dim) centries = tree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start);
768 0 : else centries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start);
769 :
770 0 : if (entries != centries) {
771 0 : delete []errors;
772 0 : delete []values;
773 0 : delete formulaTokens;
774 0 : return new TString("An ERROR has occured during fitting!");
775 : }
776 0 : values[i] = new Double_t[entries];
777 0 : memcpy(values[i], tree->GetV1(), entries*sizeof(Double_t));
778 0 : }
779 :
780 : // add points to the fitter
781 0 : for (Int_t i = 0; i < entries; i++){
782 0 : Double_t x[1000];
783 0 : for (Int_t j=0; j<dim;j++) x[j]=values[j][i];
784 0 : fitter->AddPoint(x, values[dim][i], errors[i]);
785 0 : }
786 :
787 0 : fitter->Eval();
788 0 : if (frac>0.5 && frac<1){
789 0 : fitter->EvalRobust(frac);
790 : }
791 0 : fitter->GetParameters(fitParam);
792 0 : fitter->GetCovarianceMatrix(covMatrix);
793 0 : chi2 = fitter->GetChisquare();
794 0 : npoints = entries;
795 :
796 0 : TString *preturnFormula = new TString("("), &returnFormula = *preturnFormula;
797 :
798 0 : for (Int_t iparam = 0; iparam < dim; iparam++) {
799 0 : returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam]));
800 0 : if (iparam < dim-1) returnFormula.Append("+");
801 : }
802 0 : returnFormula.Append(" )");
803 :
804 :
805 0 : for (Int_t j=0; j<dim+1;j++) delete [] values[j];
806 :
807 0 : delete formulaTokens;
808 0 : delete fitter;
809 0 : delete[] values;
810 0 : delete[] errors;
811 : return preturnFormula;
812 0 : }
813 :
814 :
815 :
816 :
817 :
818 : Int_t TStatToolkit::GetFitIndex(const TString fString, const TString subString){
819 : //
820 : // fitString - ++ separated list of fits
821 : // substring - ++ separated list of the requiered substrings
822 : //
823 : // return the last occurance of substring in fit string
824 : //
825 0 : TObjArray *arrFit = fString.Tokenize("++");
826 0 : TObjArray *arrSub = subString.Tokenize("++");
827 : Int_t index=-1;
828 0 : for (Int_t i=0; i<arrFit->GetEntries(); i++){
829 : Bool_t isOK=kTRUE;
830 0 : TString str =arrFit->At(i)->GetName();
831 0 : for (Int_t isub=0; isub<arrSub->GetEntries(); isub++){
832 0 : if (str.Contains(arrSub->At(isub)->GetName())==0) isOK=kFALSE;
833 : }
834 0 : if (isOK) index=i;
835 0 : }
836 0 : delete arrFit;
837 0 : delete arrSub;
838 0 : return index;
839 0 : }
840 :
841 :
842 : TString TStatToolkit::FilterFit(const TString &input, const TString filter, TVectorD ¶m, TMatrixD & covar){
843 : //
844 : // Filter fit expression make sub-fit
845 : //
846 0 : TObjArray *array0= input.Tokenize("++");
847 0 : TObjArray *array1= filter.Tokenize("++");
848 : //TString *presult=new TString("(0");
849 0 : TString result="(0.0";
850 0 : for (Int_t i=0; i<array0->GetEntries(); i++){
851 : Bool_t isOK=kTRUE;
852 0 : TString str(array0->At(i)->GetName());
853 0 : for (Int_t j=0; j<array1->GetEntries(); j++){
854 0 : if (str.Contains(array1->At(j)->GetName())==0) isOK=kFALSE;
855 : }
856 0 : if (isOK) {
857 0 : result+="+"+str;
858 0 : result+=Form("*(%f)",param[i+1]);
859 0 : printf("%f\t%f\t%s\n",param[i+1], TMath::Sqrt(covar(i+1,i+1)),str.Data());
860 : }
861 0 : }
862 0 : result+="-0.)";
863 0 : delete array0;
864 0 : delete array1;
865 : return result;
866 0 : }
867 :
868 : void TStatToolkit::Update1D(Double_t delta, Double_t sigma, Int_t s1, TMatrixD &vecXk, TMatrixD &covXk){
869 : //
870 : // Update parameters and covariance - with one measurement
871 : // Input:
872 : // vecXk - input vector - Updated in function
873 : // covXk - covariance matrix - Updated in function
874 : // delta, sigma, s1 - new measurement, rms of new measurement and the index of measurement
875 : const Int_t knMeas=1;
876 0 : Int_t knElem=vecXk.GetNrows();
877 :
878 0 : TMatrixD mat1(knElem,knElem); // update covariance matrix
879 0 : TMatrixD matHk(1,knElem); // vector to mesurement
880 0 : TMatrixD vecYk(knMeas,1); // Innovation or measurement residual
881 0 : TMatrixD matHkT(knElem,knMeas); // helper matrix Hk transpose
882 0 : TMatrixD matSk(knMeas,knMeas); // Innovation (or residual) covariance
883 0 : TMatrixD matKk(knElem,knMeas); // Optimal Kalman gain
884 0 : TMatrixD covXk2(knElem,knElem); // helper matrix
885 0 : TMatrixD covXk3(knElem,knElem); // helper matrix
886 0 : TMatrixD vecZk(1,1);
887 0 : TMatrixD measR(1,1);
888 0 : vecZk(0,0)=delta;
889 0 : measR(0,0)=sigma*sigma;
890 : //
891 : // reset matHk
892 0 : for (Int_t iel=0;iel<knElem;iel++)
893 0 : for (Int_t ip=0;ip<knMeas;ip++) matHk(ip,iel)=0;
894 : //mat1
895 0 : for (Int_t iel=0;iel<knElem;iel++) {
896 0 : for (Int_t jel=0;jel<knElem;jel++) mat1(iel,jel)=0;
897 0 : mat1(iel,iel)=1;
898 : }
899 : //
900 0 : matHk(0, s1)=1;
901 0 : vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual
902 0 : matHkT=matHk.T(); matHk.T();
903 0 : matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance
904 0 : matSk.Invert();
905 0 : matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
906 0 : vecXk += matKk*vecYk; // updated vector
907 0 : covXk2= (mat1-(matKk*matHk));
908 0 : covXk3 = covXk2*covXk;
909 0 : covXk = covXk3;
910 0 : Int_t nrows=covXk3.GetNrows();
911 :
912 0 : for (Int_t irow=0; irow<nrows; irow++)
913 0 : for (Int_t icol=0; icol<nrows; icol++){
914 : // rounding problems - make matrix again symteric
915 0 : covXk(irow,icol)=(covXk3(irow,icol)+covXk3(icol,irow))*0.5;
916 : }
917 0 : }
918 :
919 :
920 :
921 : void TStatToolkit::Constrain1D(const TString &input, const TString filter, TVectorD ¶m, TMatrixD & covar, Double_t mean, Double_t sigma){
922 : //
923 : // constrain linear fit
924 : // input - string description of fit function
925 : // filter - string filter to select sub fits
926 : // param,covar - parameters and covariance matrix of the fit
927 : // mean,sigma - new measurement uning which the fit is updated
928 : //
929 :
930 0 : TObjArray *array0= input.Tokenize("++");
931 0 : TObjArray *array1= filter.Tokenize("++");
932 0 : TMatrixD paramM(param.GetNrows(),1);
933 0 : for (Int_t i=0; i<=array0->GetEntries(); i++){paramM(i,0)=param(i);}
934 :
935 0 : if (filter.Length()==0){
936 0 : TStatToolkit::Update1D(mean, sigma, 0, paramM, covar);//
937 : }else{
938 0 : for (Int_t i=0; i<array0->GetEntries(); i++){
939 : Bool_t isOK=kTRUE;
940 0 : TString str(array0->At(i)->GetName());
941 0 : for (Int_t j=0; j<array1->GetEntries(); j++){
942 0 : if (str.Contains(array1->At(j)->GetName())==0) isOK=kFALSE;
943 : }
944 0 : if (isOK) {
945 0 : TStatToolkit::Update1D(mean, sigma, i+1, paramM, covar);//
946 : }
947 0 : }
948 : }
949 0 : for (Int_t i=0; i<=array0->GetEntries(); i++){
950 0 : param(i)=paramM(i,0);
951 : }
952 0 : delete array0;
953 0 : delete array1;
954 0 : }
955 :
956 : TString TStatToolkit::MakeFitString(const TString &input, const TVectorD ¶m, const TMatrixD & covar, Bool_t verbose){
957 : //
958 : //
959 : //
960 0 : TObjArray *array0= input.Tokenize("++");
961 0 : TString result=Form("(%f",param[0]);
962 0 : printf("%f\t%f\t\n", param[0], TMath::Sqrt(covar(0,0)));
963 0 : for (Int_t i=0; i<array0->GetEntries(); i++){
964 0 : TString str(array0->At(i)->GetName());
965 0 : result+="+"+str;
966 0 : result+=Form("*(%f)",param[i+1]);
967 0 : if (verbose) printf("%f\t%f\t%s\n", param[i+1], TMath::Sqrt(covar(i+1,i+1)),str.Data());
968 0 : }
969 0 : result+="-0.)";
970 0 : delete array0;
971 : return result;
972 0 : }
973 :
974 : TGraphErrors * TStatToolkit::MakeGraphErrors(TTree * tree, const char * expr, const char * cut, Int_t mstyle, Int_t mcolor, Float_t msize, Float_t offset, Int_t drawEntries, Int_t firstEntry){
975 : //
976 : // Query a graph errors
977 : // return TGraphErrors specified by expr and cut
978 : // Example usage TStatToolkit::MakeGraphError(tree,"Y:X:ErrY","X>0", 25,2,0.4)
979 : // tree - tree with variable
980 : // expr - examp
981 0 : const Int_t entries = tree->Draw(expr,cut,"goff",drawEntries,firstEntry);
982 :
983 0 : if (entries<=0) {
984 0 : ::Error("TStatToolkit::MakeGraphError","Empty or Not valid expression (%s) or cut *%s)", expr,cut);
985 0 : return 0;
986 : }
987 0 : if ( tree->GetV2()==0){
988 0 : ::Error("TStatToolkit::MakeGraphError","Not valid expression (%s) ", expr);
989 0 : return 0;
990 : }
991 : TGraphErrors * graph=0;
992 0 : if ( tree->GetV3()!=0){
993 0 : graph = new TGraphErrors (entries, tree->GetV2(),tree->GetV1(),0,tree->GetV3());
994 0 : }else{
995 0 : graph = new TGraphErrors (entries, tree->GetV2(),tree->GetV1(),0,0);
996 : }
997 :
998 0 : graph->SetMarkerStyle(mstyle);
999 0 : graph->SetMarkerColor(mcolor);
1000 0 : graph->SetLineColor(mcolor);
1001 0 : graph->SetTitle(expr);
1002 0 : TString chstring(expr);
1003 0 : TObjArray *charray = chstring.Tokenize(":");
1004 0 : graph->GetXaxis()->SetTitle(charray->At(1)->GetName());
1005 0 : graph->GetYaxis()->SetTitle(charray->At(0)->GetName());
1006 0 : THashList * metaData = (THashList*) tree->GetUserInfo()->FindObject("metaTable");
1007 0 : if (!metaData == 0){
1008 0 : TNamed *nmdTitle0 = TStatToolkit::GetMetadata(tree,Form("%s.Title",charray->At(0)->GetName()));
1009 0 : TNamed *nmdTitle1 = TStatToolkit::GetMetadata(tree,Form("%s.Title",charray->At(1)->GetName()));
1010 0 : TNamed *nmdYAxis = TStatToolkit::GetMetadata(tree,Form("%s.AxisTitle",charray->At(0)->GetName()));
1011 0 : TNamed *nmdXAxis = TStatToolkit::GetMetadata(tree,Form("%s.AxisTitle",charray->At(1)->GetName()));
1012 : //
1013 0 : TString grTitle=charray->At(0)->GetName();
1014 0 : if (nmdTitle0) grTitle=nmdTitle0->GetTitle();
1015 0 : if (nmdTitle1) {
1016 0 : grTitle+=":";
1017 0 : grTitle+=nmdTitle1->GetTitle();
1018 : }else{
1019 0 : grTitle+=":";
1020 0 : grTitle+=charray->At(1)->GetName();
1021 : }
1022 0 : if (nmdYAxis) {graph->GetYaxis()->SetTitle(nmdYAxis->GetTitle());}
1023 0 : if (nmdXAxis) {graph->GetXaxis()->SetTitle(nmdXAxis->GetTitle());}
1024 0 : graph->SetTitle(grTitle.Data());
1025 0 : }
1026 0 : delete charray;
1027 0 : if (msize>0) graph->SetMarkerSize(msize);
1028 0 : for(Int_t i=0;i<graph->GetN();i++) graph->GetX()[i]+=offset;
1029 : //
1030 0 : if (tree->GetVar(1)->IsInteger()){
1031 0 : TAxis * axis = tree->GetHistogram()->GetXaxis();
1032 0 : axis->Copy(*(graph->GetXaxis()));
1033 0 : }
1034 0 : if (tree->GetVar(0)->IsInteger()){
1035 0 : TAxis * axis = tree->GetHistogram()->GetYaxis();
1036 0 : axis->Copy(*(graph->GetYaxis()));
1037 0 : }
1038 0 : graph->Sort();
1039 : return graph;
1040 :
1041 0 : }
1042 :
1043 : THashList* TStatToolkit::AddMetadata(TTree* tree, const char *varTagName,const char *varTagValue){
1044 : //
1045 : // Add metadata infromation as user info to the tree - see https://alice.its.cern.ch/jira/browse/ATO-290
1046 : // TTree metdata are used for the Drawing methods in the folling drawing functions
1047 : /*
1048 : Supported metadata:
1049 : - <varName>.AxisTitle
1050 : - <varName>.Legend
1051 : - <varname>.Color
1052 : - <varname>.MarkerStyle
1053 : This metadata than can be used by the TStatToolkit
1054 : - TStatToolkit::MakeGraphSparse
1055 : - TStatToolkit::MakeGraphErrors
1056 : */
1057 : //
1058 0 : if (!tree) return NULL;
1059 0 : THashList * metaData = (THashList*) tree->GetUserInfo()->FindObject("metaTable");
1060 0 : if (metaData == NULL){
1061 0 : metaData=new THashList;
1062 0 : metaData->SetName("metaTable");
1063 0 : tree->GetUserInfo()->AddLast(metaData);
1064 0 : }
1065 0 : if (varTagName!=NULL && varTagValue!=NULL){
1066 0 : TNamed * named = TStatToolkit::GetMetadata(tree, varTagName);
1067 0 : if (named==NULL){
1068 0 : metaData->AddLast(new TNamed(varTagName,varTagValue));
1069 0 : }else{
1070 0 : named->SetTitle(varTagValue);
1071 : }
1072 0 : }
1073 : return metaData;
1074 0 : }
1075 :
1076 : TNamed* TStatToolkit::GetMetadata(TTree* tree, const char *vartagName){
1077 : //
1078 : // Get metadata description
1079 : //
1080 0 : if (!tree) return 0;
1081 0 : THashList * metaData = (THashList*) tree->GetUserInfo()->FindObject("metaTable");
1082 0 : if (metaData == NULL){
1083 0 : metaData=new THashList;
1084 0 : metaData->SetName("metaTable");
1085 0 : tree->GetUserInfo()->AddLast(metaData);
1086 0 : return 0;
1087 : }
1088 0 : TNamed * named = (TNamed*)metaData->FindObject(vartagName);
1089 : return named;
1090 0 : }
1091 :
1092 :
1093 :
1094 : TGraph * TStatToolkit::MakeGraphSparse(TTree * tree, const char * expr, const char * cut, Int_t mstyle, Int_t mcolor, Float_t msize, Float_t offset){
1095 : //
1096 : // Make a sparse draw of the variables
1097 : // Format of expr : Var:Run or Var:Run:ErrorY or Var:Run:ErrorY:ErrorX
1098 : // offset : points can slightly be shifted in x for better visibility with more graphs
1099 : //
1100 : // Patrick Reichelt and Marian Ivanov
1101 : // maintained and updated by Marian Ivanov
1102 0 : const Int_t entries = tree->Draw(expr,cut,"goff");
1103 0 : if (entries<=0) {
1104 0 : ::Error("TStatToolkit::MakeGraphSparse","Empty or Not valid expression (%s) or cut (%s)", expr, cut);
1105 0 : return 0;
1106 : }
1107 : // TGraph * graph = (TGraph*)gPad->GetPrimitive("Graph"); // 2D
1108 :
1109 : Double_t *graphY, *graphX;
1110 0 : graphY = tree->GetV1();
1111 0 : graphX = tree->GetV2();
1112 :
1113 : // sort according to run number
1114 0 : Int_t *index = new Int_t[entries*4];
1115 0 : TMath::Sort(entries,graphX,index,kFALSE);
1116 :
1117 : // define arrays for the new graph
1118 0 : Double_t *unsortedX = new Double_t[entries];
1119 0 : Int_t *runNumber = new Int_t[entries];
1120 : Double_t count = 0.5;
1121 :
1122 : // evaluate arrays for the new graph according to the run-number
1123 : Int_t icount=0;
1124 : //first entry
1125 0 : unsortedX[index[0]] = count;
1126 0 : runNumber[0] = graphX[index[0]];
1127 : // loop the rest of entries
1128 0 : for(Int_t i=1;i<entries;i++)
1129 : {
1130 0 : if(graphX[index[i]]==graphX[index[i-1]])
1131 0 : unsortedX[index[i]] = count;
1132 0 : else if(graphX[index[i]]!=graphX[index[i-1]]){
1133 0 : count++;
1134 0 : icount++;
1135 0 : unsortedX[index[i]] = count;
1136 0 : runNumber[icount]=graphX[index[i]];
1137 0 : }
1138 : }
1139 :
1140 : // count the number of xbins (run-wise) for the new graph
1141 0 : const Int_t newNbins = int(count+0.5);
1142 0 : Double_t *newBins = new Double_t[newNbins+1];
1143 0 : for(Int_t i=0; i<=count+1;i++){
1144 0 : newBins[i] = i;
1145 : }
1146 :
1147 : // define and fill the new graph
1148 : TGraph *graphNew = 0;
1149 0 : if (tree->GetV3()) {
1150 0 : if (tree->GetV4()) {
1151 0 : graphNew = new TGraphErrors(entries,unsortedX,graphY,tree->GetV4(),tree->GetV3());
1152 0 : }
1153 0 : else { graphNew = new TGraphErrors(entries,unsortedX,graphY,0,tree->GetV3()); }
1154 : }
1155 0 : else { graphNew = new TGraphErrors(entries,unsortedX,graphY,0,0); }
1156 : // with "Set(...)", the x-axis is being sorted
1157 0 : graphNew->GetXaxis()->Set(newNbins,newBins);
1158 :
1159 : // set the bins for the x-axis, apply shifting of points
1160 0 : Char_t xName[50];
1161 0 : for(Int_t i=0;i<count;i++){
1162 0 : snprintf(xName,50,"%d",runNumber[i]);
1163 0 : graphNew->GetXaxis()->SetBinLabel(i+1,xName);
1164 0 : graphNew->GetX()[i]+=offset;
1165 : }
1166 0 : if (tree->GetVar(1)->IsInteger() && strlen(tree->GetHistogram()->GetXaxis()->GetBinLabel(1))>0){
1167 0 : for(Int_t i=0;i<count;i++){
1168 0 : graphNew->GetXaxis()->SetBinLabel(i+1,tree->GetHistogram()->GetXaxis()->GetBinLabel(i+1));
1169 : }
1170 0 : }
1171 0 : if (tree->GetVar(0)->IsInteger() && strlen(tree->GetHistogram()->GetXaxis()->GetBinLabel(1))>0 ){
1172 0 : for(Int_t i=0;i<count;i++){
1173 0 : graphNew->GetYaxis()->SetBinLabel(i+1,tree->GetHistogram()->GetYaxis()->GetBinLabel(i+1));
1174 : }
1175 0 : }
1176 :
1177 :
1178 :
1179 :
1180 0 : graphNew->GetHistogram()->SetTitle("");
1181 0 : graphNew->SetMarkerStyle(mstyle);
1182 0 : graphNew->SetMarkerColor(mcolor); graphNew->SetLineColor(mcolor);
1183 0 : if (msize>0) { graphNew->SetMarkerSize(msize); graphNew->SetLineWidth(msize); }
1184 0 : delete [] unsortedX;
1185 0 : delete [] runNumber;
1186 0 : delete [] index;
1187 0 : delete [] newBins;
1188 : //
1189 0 : TString chstring(expr);
1190 0 : if (cut) chstring+=TString::Format(" ( %s )", cut);
1191 0 : graphNew->SetTitle(chstring);
1192 :
1193 0 : THashList * metaData = (THashList*) tree->GetUserInfo()->FindObject("metaTable");
1194 0 : if (!metaData == 0){
1195 0 : chstring=expr;
1196 0 : TObjArray *charray = chstring.Tokenize(":");
1197 0 : graphNew->GetXaxis()->SetTitle(charray->At(1)->GetName());
1198 0 : graphNew->GetYaxis()->SetTitle(charray->At(0)->GetName());
1199 0 : TNamed *nmdTitle0 = TStatToolkit::GetMetadata(tree,Form("%s.Title",charray->At(0)->GetName()));
1200 0 : TNamed *nmdTitle1 = TStatToolkit::GetMetadata(tree,Form("%s.Title",charray->At(1)->GetName()));
1201 0 : TNamed *nmdYAxis = TStatToolkit::GetMetadata(tree,Form("%s.AxisTitle",charray->At(0)->GetName()));
1202 0 : TNamed *nmdXAxis = TStatToolkit::GetMetadata(tree,Form("%s.AxisTitle",charray->At(1)->GetName()));
1203 : //
1204 0 : TString grTitle=charray->At(0)->GetName();
1205 0 : if (nmdTitle0) grTitle=nmdTitle0->GetTitle();
1206 0 : if (nmdTitle1) {
1207 0 : grTitle+=":";
1208 0 : grTitle+=nmdTitle1->GetTitle();
1209 : }else{
1210 0 : grTitle+=":";
1211 0 : grTitle+=charray->At(1)->GetName();
1212 : }
1213 0 : if (cut) grTitle+=TString::Format(" ( %s )", cut);
1214 0 : graphNew->SetTitle(grTitle);
1215 0 : if (nmdYAxis) {graphNew->GetYaxis()->SetTitle(nmdYAxis->GetTitle());}
1216 0 : if (nmdXAxis) {graphNew->GetXaxis()->SetTitle(nmdXAxis->GetTitle());}
1217 0 : delete charray;
1218 0 : }
1219 : return graphNew;
1220 0 : }
1221 :
1222 :
1223 :
1224 : //
1225 : // functions used for the trending
1226 : //
1227 :
1228 :
1229 : Int_t TStatToolkit::MakeStatAlias(TTree * tree, const char * expr, const char * cut, const char * alias)
1230 : {
1231 : //
1232 : // Add alias using statistical values of a given variable.
1233 : // (by MI, Patrick Reichelt)
1234 : //
1235 : // tree - input tree
1236 : // expr - variable expression
1237 : // cut - selection criteria
1238 : // Output - return number of entries used to define variable
1239 : // In addition mean, rms, median, and robust mean and rms (choosing fraction of data with smallest RMS)
1240 : //
1241 : /* Example usage:
1242 : 1.) create the robust estimators for variable expr="QA.TPC.CPass1.meanTPCncl" and create a corresponding
1243 : aliases with the prefix alias[0]="ncl", calculated using fraction alias[1]="0.90"
1244 :
1245 : TStatToolkit::MakeStatAlias(tree,"QA.TPC.CPass1.meanTPCncl","QA.TPC.CPass1.status>0","ncl:0.9");
1246 : root [4] tree->GetListOfAliases().Print()
1247 : OBJ: TNamed ncl_Median (130.964333+0)
1248 : OBJ: TNamed ncl_Mean (122.120387+0)
1249 : OBJ: TNamed ncl_RMS (33.509623+0)
1250 : OBJ: TNamed ncl_Mean90 (131.503862+0)
1251 : OBJ: TNamed ncl_RMS90 (3.738260+0)
1252 : */
1253 : //
1254 0 : Int_t entries = tree->Draw(expr,cut,"goff");
1255 0 : if (entries<=1){
1256 0 : printf("Expression or cut not valid:\t%s\t%s\n", expr, cut);
1257 0 : return 0;
1258 : }
1259 : //
1260 0 : TObjArray* oaAlias = TString(alias).Tokenize(":");
1261 0 : if (oaAlias->GetEntries()<2) {
1262 0 : printf("Alias must have 2 arguments:\t%s\n", alias);
1263 0 : return 0;
1264 : }
1265 0 : Float_t entryFraction = atof( oaAlias->At(1)->GetName() );
1266 : //
1267 0 : Double_t median = TMath::Median(entries,tree->GetV1());
1268 0 : Double_t mean = TMath::Mean(entries,tree->GetV1());
1269 0 : Double_t rms = TMath::RMS(entries,tree->GetV1());
1270 0 : Double_t meanEF=0, rmsEF=0;
1271 0 : TStatToolkit::EvaluateUni(entries, tree->GetV1(), meanEF, rmsEF, entries*entryFraction);
1272 : //
1273 0 : tree->SetAlias(Form("%s_Median",oaAlias->At(0)->GetName()), Form("(%f+0)",median));
1274 0 : tree->SetAlias(Form("%s_Mean",oaAlias->At(0)->GetName()), Form("(%f+0)",mean));
1275 0 : tree->SetAlias(Form("%s_RMS",oaAlias->At(0)->GetName()), Form("(%f+0)",rms));
1276 0 : tree->SetAlias(Form("%s_Mean%d",oaAlias->At(0)->GetName(),Int_t(entryFraction*100)), Form("(%f+0)",meanEF));
1277 0 : tree->SetAlias(Form("%s_RMS%d",oaAlias->At(0)->GetName(),Int_t(entryFraction*100)), Form("(%f+0)",rmsEF));
1278 0 : delete oaAlias;
1279 : return entries;
1280 0 : }
1281 :
1282 : Int_t TStatToolkit::SetStatusAlias(TTree * tree, const char * expr, const char * cut, const char * alias)
1283 : {
1284 : //
1285 : // Add alias to trending tree using statistical values of a given variable.
1286 : // (by MI, Patrick Reichelt)
1287 : //
1288 : // format of expr : varname (e.g. meanTPCncl)
1289 : // format of cut : char like in TCut
1290 : // format of alias: alias:query:entryFraction(EF) (fraction of entries used for uniformity evaluation)
1291 : // e.g.: varname_Out:(abs(varname-meanEF)>6.*rmsEF):0.8
1292 : // available internal variables are: 'varname, Median, Mean, MeanEF, RMS, RMSEF'
1293 : // in the alias, 'varname' will be replaced by its content, and 'EF' by the percentage (e.g. MeanEF -> Mean80)
1294 : //
1295 : /* Example usage:
1296 : 1.) Define robust mean (possible, but easier done with TStatToolkit::MakeStatAlias(...))
1297 : TStatToolkit::SetStatusAlias(tree, "meanTPCnclF", "meanTPCnclF>0", "meanTPCnclF_MeanEF:MeanEF:0.80") ;
1298 : root [10] tree->GetListOfAliases()->Print()
1299 : Collection name='TList', class='TList', size=1
1300 : OBJ: TNamed meanTPCnclF_Mean80 0.899308
1301 : 2.) create alias outlyers - 6 sigma cut
1302 : TStatToolkit::SetStatusAlias(tree, "meanTPCnclF", "meanTPCnclF>0", "meanTPCnclF_Out:(abs(meanTPCnclF-MeanEF)>6.*RMSEF):0.8")
1303 : meanTPCnclF_Out ==> (abs(meanTPCnclF-0.899308)>6.*0.016590)
1304 : 3.) the same functionality as in 2.)
1305 : TStatToolkit::SetStatusAlias(tree, "meanTPCnclF", "meanTPCnclF>0", "varname_Out2:(abs(varname-MeanEF)>6.*RMSEF):0.8")
1306 : meanTPCnclF_Out2 ==> (abs(meanTPCnclF-0.899308)>6.*0.016590)
1307 : */
1308 : //
1309 0 : Int_t entries = tree->Draw(expr,cut,"goff");
1310 0 : if (entries<1){
1311 0 : printf("Expression or cut not valid:\t%s\t%s\n", expr, cut);
1312 0 : return 0;
1313 : }
1314 : //
1315 0 : TObjArray* oaVar = TString(expr).Tokenize(":");
1316 0 : char varname[50];
1317 0 : snprintf(varname,50,"%s", oaVar->At(0)->GetName());
1318 : Float_t entryFraction = 0.8;
1319 : //
1320 0 : TObjArray* oaAlias = TString(alias).Tokenize(":");
1321 0 : if (oaAlias->GetEntries()<2) {
1322 0 : printf("Alias must have at least 2 arguments:\t%s\n", alias);
1323 0 : return 0;
1324 : }
1325 0 : else if (oaAlias->GetEntries()<3) {
1326 : //printf("Using default entryFraction if needed:\t%f\n", entryFraction);
1327 : }
1328 0 : else entryFraction = atof( oaAlias->At(2)->GetName() );
1329 : //
1330 0 : Double_t median = TMath::Median(entries,tree->GetV1());
1331 0 : Double_t mean = TMath::Mean(entries,tree->GetV1());
1332 0 : Double_t rms = TMath::RMS(entries,tree->GetV1());
1333 0 : Double_t meanEF=0, rmsEF=0;
1334 0 : TStatToolkit::EvaluateUni(entries, tree->GetV1(), meanEF, rmsEF, entries*entryFraction);
1335 : //
1336 0 : TString sAlias( oaAlias->At(0)->GetName() );
1337 0 : sAlias.ReplaceAll("varname",varname);
1338 0 : sAlias.ReplaceAll("MeanEF", Form("Mean%1.0f",entryFraction*100) );
1339 0 : sAlias.ReplaceAll("RMSEF", Form("RMS%1.0f",entryFraction*100) );
1340 0 : TString sQuery( oaAlias->At(1)->GetName() );
1341 0 : sQuery.ReplaceAll("varname",varname);
1342 0 : sQuery.ReplaceAll("MeanEF", Form("%f",meanEF) );
1343 0 : sQuery.ReplaceAll("RMSEF", Form("%f",rmsEF) ); //make sure to replace 'RMSEF' before 'RMS'...
1344 0 : sQuery.ReplaceAll("Median", Form("%f",median) );
1345 0 : sQuery.ReplaceAll("Mean", Form("%f",mean) );
1346 0 : sQuery.ReplaceAll("RMS", Form("%f",rms) );
1347 0 : printf("define alias:\t%s = %s\n", sAlias.Data(), sQuery.Data());
1348 : //
1349 0 : char query[200];
1350 0 : char aname[200];
1351 0 : snprintf(query,200,"%s", sQuery.Data());
1352 0 : snprintf(aname,200,"%s", sAlias.Data());
1353 0 : tree->SetAlias(aname, query);
1354 0 : delete oaVar;
1355 0 : delete oaAlias;
1356 : return entries;
1357 0 : }
1358 :
1359 : TMultiGraph* TStatToolkit::MakeStatusMultGr(TTree * tree, const char * expr, const char * cut, const char * alias, Int_t igr)
1360 : {
1361 : //
1362 : // Compute a trending multigraph that shows for which runs a variable has outliers.
1363 : // (by MI, Patrick Reichelt)
1364 : //
1365 : // format of expr : varname:xaxis (e.g. meanTPCncl:run, but 'varname' can be any string that you need for seach-and-replace)
1366 : // format of cut : char like in TCut
1367 : // format of alias: (1):(statisticOK):(varname_Warning):(varname_Out)[:(varname_PhysAcc):(varname_Extra)]
1368 : //
1369 : // function MakeGraphSparse() is called for each alias argument, which will be used as tree expression.
1370 : // each alias argument is supposed to be a Boolean statement which can be evaluated as tree expression.
1371 : // the order of these criteria should be kept, as the marker styles and colors are chosen to be meaningful this way!
1372 : // 'statisticOK' could e.g. be an alias for '(meanTPCncl>0)'.
1373 : // if you dont need e.g. a 'warning' condition, then just replace it by (0).
1374 : // in the alias, 'varname' will be replaced by its content (e.g. varname_Out -> meanTPCncl_Out)
1375 : // note: the aliases 'varname_Out' etc have to be defined by function TStatToolkit::SetStatusAlias(...)
1376 : // counter igr is used to shift the multigraph in y when filling a TObjArray.
1377 : //
1378 : //
1379 : // To create the Status Bar, the following is done in principle.
1380 : // ( example current usage in $ALICE_ROOT/PWGPP/TPC/macros/drawPerformanceTPCQAMatchTrends.C and ./qaConfig.C. )
1381 : //
1382 : // TStatToolkit::SetStatusAlias(tree, "meanTPCncl", "", "varname_Out:(abs(varname-MeanEF)>6.*RMSEF):0.8");
1383 : // TStatToolkit::SetStatusAlias(tree, "tpcItsMatchA", "", "varname_Out:(abs(varname-MeanEF)>6.*RMSEF):0.8");
1384 : // TStatToolkit::SetStatusAlias(tree, "meanTPCncl", "", "varname_Warning:(abs(varname-MeanEF)>3.*RMSEF):0.8");
1385 : // TStatToolkit::SetStatusAlias(tree, "tpcItsMatchA", "", "varname_Warning:(abs(varname-MeanEF)>3.*RMSEF):0.8");
1386 : // TObjArray* oaMultGr = new TObjArray(); int igr=0;
1387 : // oaMultGr->Add( TStatToolkit::MakeStatusMultGr(tree, "tpcItsMatchA:run", "", "(1):(meanTPCncl>0):(varname_Warning):(varname_Outlier):", igr) ); igr++;
1388 : // oaMultGr->Add( TStatToolkit::MakeStatusMultGr(tree, "meanTPCncl:run", "", "(1):(meanTPCncl>0):(varname_Warning):(varname_Outlier):", igr) ); igr++;
1389 : // TCanvas *c1 = new TCanvas("c1","c1");
1390 : // TStatToolkit::AddStatusPad(c1, 0.30, 0.40);
1391 : // TStatToolkit::DrawStatusGraphs(oaMultGr);
1392 :
1393 :
1394 0 : TObjArray* oaVar = TString(expr).Tokenize(":");
1395 0 : if (oaVar->GetEntries()<2) {
1396 0 : printf("Expression has to be of type 'varname:xaxis':\t%s\n", expr);
1397 0 : return 0;
1398 : }
1399 0 : char varname[50];
1400 0 : char var_x[50];
1401 0 : snprintf(varname,50,"%s", oaVar->At(0)->GetName());
1402 0 : snprintf(var_x ,50,"%s", oaVar->At(1)->GetName());
1403 : //
1404 0 : TString sAlias(alias);
1405 0 : sAlias.ReplaceAll("varname",varname);
1406 0 : TObjArray* oaAlias = TString(sAlias.Data()).Tokenize(":");
1407 0 : if (oaAlias->GetEntries()<2) {
1408 0 : printf("Alias must have 2-6 arguments:\t%s\n", alias);
1409 0 : return 0;
1410 : }
1411 0 : char query[200];
1412 0 : TMultiGraph* multGr = new TMultiGraph();
1413 0 : Int_t marArr[6] = {24+igr%2, 20+igr%2, 20+igr%2, 20+igr%2, 20+igr%2, 20+igr%2};
1414 0 : Int_t colArr[6] = {kBlack, kBlack, kOrange, kRed, kGreen+1, kBlue};
1415 0 : Double_t sizeArr[6] = {1.4, 1.1, 1.5, 1.1, 1.4, 0.8};
1416 0 : Double_t shiftArr[6] = {0., 0., 0.25, 0.25, -0.25, -0.25};
1417 0 : const Int_t ngr = oaAlias->GetEntriesFast();
1418 0 : for (Int_t i=0; i<ngr; i++){
1419 0 : snprintf(query,200, "%f*(%s-0.5):%s", 1.+igr, oaAlias->At(i)->GetName(), var_x);
1420 0 : TGraphErrors * gr = (TGraphErrors*) TStatToolkit::MakeGraphSparse(tree,query,cut,marArr[i],colArr[i],sizeArr[i],shiftArr[i]);
1421 0 : if (gr) multGr->Add(gr);
1422 : }
1423 : //
1424 0 : multGr->SetName(varname);
1425 0 : multGr->SetTitle(varname); // used for y-axis labels of status bar, can be modified by calling function.
1426 0 : delete oaVar;
1427 0 : delete oaAlias;
1428 : return multGr;
1429 0 : }
1430 :
1431 :
1432 : void TStatToolkit::AddStatusPad(TCanvas* c1, Float_t padratio, Float_t bottommargin)
1433 : {
1434 : //
1435 : // add pad to bottom of canvas for Status graphs (by Patrick Reichelt)
1436 : // call function "DrawStatusGraphs(...)" afterwards
1437 : //
1438 0 : TCanvas* c1_clone = (TCanvas*) c1->Clone("c1_clone");
1439 0 : c1->Clear();
1440 : // produce new pads
1441 0 : c1->cd();
1442 0 : TPad* pad1 = new TPad("pad1", "pad1", 0., padratio, 1., 1.);
1443 0 : pad1->Draw();
1444 0 : pad1->SetNumber(1); // so it can be called via "c1->cd(1);"
1445 0 : c1->cd();
1446 0 : TPad* pad2 = new TPad("pad2", "pad2", 0., 0., 1., padratio);
1447 0 : pad2->Draw();
1448 0 : pad2->SetNumber(2);
1449 : // draw original canvas into first pad
1450 0 : c1->cd(1);
1451 0 : c1_clone->DrawClonePad();
1452 0 : pad1->SetBottomMargin(0.001);
1453 0 : pad1->SetRightMargin(0.01);
1454 : // set up second pad
1455 0 : c1->cd(2);
1456 0 : pad2->SetGrid(3);
1457 0 : pad2->SetTopMargin(0);
1458 0 : pad2->SetBottomMargin(bottommargin); // for the long x-axis labels (runnumbers)
1459 0 : pad2->SetRightMargin(0.01);
1460 0 : }
1461 :
1462 :
1463 : void TStatToolkit::DrawStatusGraphs(TObjArray* oaMultGr)
1464 : {
1465 : //
1466 : // draw Status graphs into active pad of canvas (by MI, Patrick Reichelt)
1467 : // ...into bottom pad, if called after "AddStatusPad(...)"
1468 : //
1469 0 : const Int_t nvars = oaMultGr->GetEntriesFast();
1470 0 : TGraph* grAxis = (TGraph*) ((TMultiGraph*) oaMultGr->At(0))->GetListOfGraphs()->At(0);
1471 0 : grAxis->SetMaximum(0.5*nvars+0.5);
1472 0 : grAxis->SetMinimum(0);
1473 0 : grAxis->GetYaxis()->SetLabelSize(0);
1474 0 : grAxis->GetYaxis()->SetTitle("");
1475 0 : grAxis->SetTitle("");
1476 0 : Int_t entries = grAxis->GetN();
1477 0 : grAxis->GetXaxis()->SetLabelSize(5.7*TMath::Min(TMath::Max(5./entries,0.01),0.03));
1478 0 : grAxis->GetXaxis()->LabelsOption("v");
1479 0 : grAxis->Draw("ap");
1480 : //
1481 : // draw multigraphs & names of status variables on the y axis
1482 0 : for (Int_t i=0; i<nvars; i++){
1483 0 : ((TMultiGraph*) oaMultGr->At(i))->Draw("p");
1484 0 : TLatex* ylabel = new TLatex(-0.1, 0.5*i+0.5, ((TMultiGraph*) oaMultGr->At(i))->GetTitle());
1485 0 : ylabel->SetTextAlign(32); //hor:right & vert:centered
1486 0 : ylabel->SetTextSize(0.025/gPad->GetHNDC());
1487 0 : ylabel->Draw();
1488 : }
1489 0 : }
1490 :
1491 :
1492 : TTree* TStatToolkit::WriteStatusToTree(TObject* oStatusGr)
1493 : {
1494 : //
1495 : // Create Tree with Integers for each status variable flag (warning, outlier, physacc).
1496 : // (by Patrick Reichelt)
1497 : //
1498 : // input: either a TMultiGraph with status of single variable, which
1499 : // was computed by TStatToolkit::MakeStatusMultGr(),
1500 : // or a TObjArray which contains up to 10 of such variables.
1501 : // example: TTree* statusTree = WriteStatusToTree( TStatToolkit::MakeStatusMultGr(tree, "tpcItsMatch:run", "", sCriteria.Data(), 0) );
1502 : // or : TTree* statusTree = TStatToolkit::WriteStatusToTree(oaMultGr);
1503 : //
1504 : // output tree: 1=flag is true, 0=flag is false, -1=flag was not computed.
1505 : // To be rewritten to the pcstream
1506 :
1507 : TObjArray* oaMultGr = NULL;
1508 : Bool_t needDeletion=kFALSE;
1509 0 : if (oStatusGr->IsA() == TObjArray::Class()) {
1510 0 : oaMultGr = (TObjArray*) oStatusGr;
1511 0 : }
1512 0 : else if (oStatusGr->IsA() == TMultiGraph::Class()) {
1513 0 : oaMultGr = new TObjArray(); needDeletion=kTRUE;
1514 0 : oaMultGr->Add((TMultiGraph*) oStatusGr);
1515 : }
1516 : else {
1517 0 : Printf("WriteStatusToTree(): Error! 'oStatusGr' must be a TMultiGraph or a TObjArray of them!");
1518 0 : return 0;
1519 : }
1520 : // variables for output tree
1521 : const int nvarsMax=10;
1522 : const int ncritMax=5;
1523 0 : Int_t currentRun;
1524 0 : Int_t treevars[nvarsMax*ncritMax];
1525 0 : TString varnames[nvarsMax*ncritMax];
1526 0 : for (int i=0; i<nvarsMax*ncritMax; i++) treevars[i]=-1;
1527 :
1528 0 : Printf("WriteStatusToTree(): writing following variables to TTree (maybe only subset of listed criteria filled)");
1529 0 : for (Int_t vari=0; vari<nvarsMax; vari++)
1530 : {
1531 0 : if (vari < oaMultGr->GetEntriesFast()) {
1532 0 : varnames[vari*ncritMax+0] = Form("%s_statisticOK", ((TMultiGraph*) oaMultGr->At(vari))->GetName());
1533 0 : varnames[vari*ncritMax+1] = Form("%s_Warning", ((TMultiGraph*) oaMultGr->At(vari))->GetName());
1534 0 : varnames[vari*ncritMax+2] = Form("%s_Outlier", ((TMultiGraph*) oaMultGr->At(vari))->GetName());
1535 0 : varnames[vari*ncritMax+3] = Form("%s_PhysAcc", ((TMultiGraph*) oaMultGr->At(vari))->GetName());
1536 0 : varnames[vari*ncritMax+4] = Form("%s_Extra", ((TMultiGraph*) oaMultGr->At(vari))->GetName());
1537 : }
1538 : else {
1539 0 : varnames[vari*ncritMax+0] = Form("dummy");
1540 0 : varnames[vari*ncritMax+1] = Form("dummy");
1541 0 : varnames[vari*ncritMax+2] = Form("dummy");
1542 0 : varnames[vari*ncritMax+3] = Form("dummy");
1543 0 : varnames[vari*ncritMax+4] = Form("dummy");
1544 : }
1545 0 : 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;
1546 : }
1547 :
1548 0 : TTree* statusTree = new TTree("statusTree","statusTree");
1549 0 : statusTree->Branch("run", ¤tRun );
1550 0 : statusTree->Branch(varnames[ 0].Data(), &treevars[ 0]);
1551 0 : statusTree->Branch(varnames[ 1].Data(), &treevars[ 1]);
1552 0 : statusTree->Branch(varnames[ 2].Data(), &treevars[ 2]);
1553 0 : statusTree->Branch(varnames[ 3].Data(), &treevars[ 3]);
1554 0 : statusTree->Branch(varnames[ 4].Data(), &treevars[ 4]);
1555 0 : statusTree->Branch(varnames[ 5].Data(), &treevars[ 5]);
1556 0 : statusTree->Branch(varnames[ 6].Data(), &treevars[ 6]);
1557 0 : statusTree->Branch(varnames[ 7].Data(), &treevars[ 7]);
1558 0 : statusTree->Branch(varnames[ 8].Data(), &treevars[ 8]);
1559 0 : statusTree->Branch(varnames[ 9].Data(), &treevars[ 9]);
1560 0 : statusTree->Branch(varnames[10].Data(), &treevars[10]);
1561 0 : statusTree->Branch(varnames[11].Data(), &treevars[11]);
1562 0 : statusTree->Branch(varnames[12].Data(), &treevars[12]);
1563 0 : statusTree->Branch(varnames[13].Data(), &treevars[13]);
1564 0 : statusTree->Branch(varnames[14].Data(), &treevars[14]);
1565 0 : statusTree->Branch(varnames[15].Data(), &treevars[15]);
1566 0 : statusTree->Branch(varnames[16].Data(), &treevars[16]);
1567 0 : statusTree->Branch(varnames[17].Data(), &treevars[17]);
1568 0 : statusTree->Branch(varnames[18].Data(), &treevars[18]);
1569 0 : statusTree->Branch(varnames[19].Data(), &treevars[19]);
1570 0 : statusTree->Branch(varnames[20].Data(), &treevars[20]);
1571 0 : statusTree->Branch(varnames[21].Data(), &treevars[21]);
1572 0 : statusTree->Branch(varnames[22].Data(), &treevars[22]);
1573 0 : statusTree->Branch(varnames[23].Data(), &treevars[23]);
1574 0 : statusTree->Branch(varnames[24].Data(), &treevars[24]);
1575 0 : statusTree->Branch(varnames[25].Data(), &treevars[25]);
1576 0 : statusTree->Branch(varnames[26].Data(), &treevars[26]);
1577 0 : statusTree->Branch(varnames[27].Data(), &treevars[27]);
1578 0 : statusTree->Branch(varnames[28].Data(), &treevars[28]);
1579 0 : statusTree->Branch(varnames[29].Data(), &treevars[29]);
1580 0 : statusTree->Branch(varnames[30].Data(), &treevars[30]);
1581 0 : statusTree->Branch(varnames[31].Data(), &treevars[31]);
1582 0 : statusTree->Branch(varnames[32].Data(), &treevars[32]);
1583 0 : statusTree->Branch(varnames[33].Data(), &treevars[33]);
1584 0 : statusTree->Branch(varnames[34].Data(), &treevars[34]);
1585 0 : statusTree->Branch(varnames[35].Data(), &treevars[35]);
1586 0 : statusTree->Branch(varnames[36].Data(), &treevars[36]);
1587 0 : statusTree->Branch(varnames[37].Data(), &treevars[37]);
1588 0 : statusTree->Branch(varnames[38].Data(), &treevars[38]);
1589 0 : statusTree->Branch(varnames[39].Data(), &treevars[39]);
1590 0 : statusTree->Branch(varnames[40].Data(), &treevars[40]);
1591 0 : statusTree->Branch(varnames[41].Data(), &treevars[41]);
1592 0 : statusTree->Branch(varnames[42].Data(), &treevars[42]);
1593 0 : statusTree->Branch(varnames[43].Data(), &treevars[43]);
1594 0 : statusTree->Branch(varnames[44].Data(), &treevars[44]);
1595 0 : statusTree->Branch(varnames[45].Data(), &treevars[45]);
1596 0 : statusTree->Branch(varnames[46].Data(), &treevars[46]);
1597 0 : statusTree->Branch(varnames[47].Data(), &treevars[47]);
1598 0 : statusTree->Branch(varnames[48].Data(), &treevars[48]);
1599 0 : statusTree->Branch(varnames[49].Data(), &treevars[49]);
1600 :
1601 : // run loop
1602 0 : Double_t graphX; // x-position of marker (0.5, 1.5, ...)
1603 0 : Double_t graphY; // if >0 -> warning/outlier/physacc! if =-0.5 -> no warning/outlier/physacc
1604 0 : TList* arrRuns = (TList*) ((TGraph*) ((TMultiGraph*) oaMultGr->At(0))->GetListOfGraphs()->At(0))->GetXaxis()->GetLabels();
1605 : //'TAxis->GetLabels()' returns THashList of TObjString, but using THashList gives compilation error "... incomplete type 'struct THashList' "
1606 0 : for (Int_t runi=0; runi<arrRuns->GetSize(); runi++)
1607 : {
1608 0 : currentRun = atoi( arrRuns->At(runi)->GetName() );
1609 : //Printf(" runi=%2i, name: %s \t run number: %i", runi, arrRuns->At(runi)->GetName(), currentRun);
1610 :
1611 : // status variable loop
1612 0 : for (Int_t vari=0; vari<oaMultGr->GetEntriesFast(); vari++)
1613 : {
1614 0 : TMultiGraph* multGr = (TMultiGraph*) oaMultGr->At(vari);
1615 :
1616 : // criteria loop
1617 : // the order is given by TStatToolkit::MakeStatusMultGr().
1618 : // criterion #1 is 'statisticOK' and mandatory, the rest is optional. (#0 is always True, thus skipped)
1619 0 : for (Int_t criti=1; criti<multGr->GetListOfGraphs()->GetEntries(); criti++)
1620 : {
1621 0 : TGraph* grCriterion = (TGraph*) multGr->GetListOfGraphs()->At(criti);
1622 0 : graphX = -1, graphY = -1;
1623 0 : grCriterion->GetPoint(runi, graphX, graphY);
1624 0 : treevars[(vari)*ncritMax+(criti-1)] = (graphY>0)?1:0;
1625 : }
1626 : }
1627 0 : statusTree->Fill();
1628 : }
1629 :
1630 0 : if (needDeletion) delete oaMultGr;
1631 :
1632 : return statusTree;
1633 0 : }
1634 :
1635 :
1636 : void TStatToolkit::MakeSummaryTree(TTree* treeIn, TTreeSRedirector *pcstream, TObjString & sumID, TCut &selection){
1637 : //
1638 : // Make a summary tree for the input tree
1639 : // For the moment statistic works only for the primitive branches (Float/Double/Int)
1640 : // Extension recursive version planned for graphs a and histograms
1641 : //
1642 : // Following statistics are exctracted:
1643 : // - Standard: mean, meadian, rms
1644 : // - LTM robust statistic: mean60, rms60, mean90, rms90
1645 : // Parameters:
1646 : // treeIn - input tree
1647 : // pctream - Output redirector
1648 : // sumID - ID as will be used in output tree
1649 : // selection - selection criteria define the set of entries used to evaluat statistic
1650 : //
1651 : // Curently only predefined statistic used to fill summary information
1652 : // Future plans an option user defined statistic descriptor instead of the defualt (if exist)
1653 : //
1654 : // e.g
1655 : // default.Branches=median:mean90:rms90:mean60:rms60
1656 : // interactionRate.Branches mean90:median:rms90:mean95:rms95:mean60:rms60
1657 : //
1658 0 : TObjArray * brArray = treeIn->GetListOfBranches();
1659 0 : Int_t tEntries= treeIn->GetEntries();
1660 0 : Int_t nBranches=brArray->GetEntries();
1661 0 : TString treeName = treeIn->GetName();
1662 :
1663 0 : (*pcstream)<<treeName.Data()<<"entries="<<tEntries;
1664 0 : (*pcstream)<<treeName.Data()<<"ID.="<<&sumID;
1665 :
1666 0 : TMatrixD valBranch(nBranches,7);
1667 0 : for (Int_t iBr=0; iBr<nBranches; iBr++){
1668 0 : TString brName= brArray->At(iBr)->GetName();
1669 0 : Int_t entries=treeIn->Draw(TString::Format("%s>>dummy(10,0,1)",brArray->At(iBr)->GetName()).Data(),selection,"goff");
1670 0 : if (entries==0) continue;
1671 0 : Double_t median, mean, rms, mean60,rms60, mean90, rms90;
1672 0 : mean = TMath::Mean(entries,treeIn->GetV1());
1673 0 : median= TMath::Median(entries,treeIn->GetV1());
1674 0 : rms = TMath::RMS(entries,treeIn->GetV1());
1675 0 : TStatToolkit::EvaluateUni(entries, treeIn->GetV1(), mean60,rms60,TMath::Min(TMath::Max(2., 0.60*entries),Double_t(entries)));
1676 0 : TStatToolkit::EvaluateUni(entries, treeIn->GetV1(), mean90,rms90,TMath::Min(TMath::Max(2., 0.90*entries),Double_t(entries)));
1677 0 : valBranch(iBr,0)=mean;
1678 0 : valBranch(iBr,1)=median;
1679 0 : valBranch(iBr,2)=rms;
1680 0 : valBranch(iBr,3)=mean60;
1681 0 : valBranch(iBr,4)=rms60;
1682 0 : valBranch(iBr,5)=mean90;
1683 0 : valBranch(iBr,6)=rms90;
1684 0 : (*pcstream)<<treeName.Data()<<
1685 0 : brName+"="<<valBranch(iBr,1)<< // use as an default median estimator
1686 0 : brName+"_Mean="<<valBranch(iBr,0)<<
1687 0 : brName+"_Median="<<valBranch(iBr,1)<<
1688 0 : brName+"_RMS="<<valBranch(iBr,2)<<
1689 0 : brName+"_Mean60="<<valBranch(iBr,3)<<
1690 0 : brName+"_RMS60="<<valBranch(iBr,4)<<
1691 0 : brName+"_Mean90="<<valBranch(iBr,5)<<
1692 0 : brName+"_RMS90="<<valBranch(iBr,6);
1693 0 : }
1694 0 : (*pcstream)<<treeName.Data()<<"\n";
1695 0 : }
1696 :
1697 :
1698 :
1699 : TMultiGraph* TStatToolkit::MakeStatusLines(TTree * tree, const char * expr, const char * cut, const char * alias)
1700 : {
1701 : //
1702 : // Create status lines for trending using MakeGraphSparse(), very similar to MakeStatusMultGr().
1703 : // (by Patrick Reichelt)
1704 : //
1705 : // format of expr : varname:xaxis (e.g. meanTPCncl:run, but 'varname' can be any string that you need for seach-and-replace)
1706 : // format of cut : char like in TCut
1707 : // format of alias: varname_OutlierMin:varname_OutlierMax:varname_WarningMin:varname_WarningMax:varname_PhysAccMin:varname_PhysAccMax:varname_RobustMean
1708 : //
1709 0 : TObjArray* oaVar = TString(expr).Tokenize(":");
1710 0 : if (oaVar->GetEntries()<2) {
1711 0 : printf("Expression has to be of type 'varname:xaxis':\t%s\n", expr);
1712 0 : return 0;
1713 : }
1714 0 : char varname[50];
1715 0 : char var_x[50];
1716 0 : snprintf(varname,50,"%s", oaVar->At(0)->GetName());
1717 0 : snprintf(var_x ,50,"%s", oaVar->At(1)->GetName());
1718 : //
1719 0 : TString sAlias(alias);
1720 0 : if (sAlias.IsNull()) { // alias for default usage set here:
1721 0 : sAlias = "varname_OutlierMin:varname_OutlierMax:varname_WarningMin:varname_WarningMax:varname_PhysAccMin:varname_PhysAccMax:varname_RobustMean";
1722 : }
1723 0 : sAlias.ReplaceAll("varname",varname);
1724 0 : TObjArray* oaAlias = TString(sAlias.Data()).Tokenize(":");
1725 0 : if (oaAlias->GetEntries()<2) {
1726 0 : printf("Alias must have 2-7 arguments:\t%s\n", alias);
1727 0 : return 0;
1728 : }
1729 0 : char query[200];
1730 0 : TMultiGraph* multGr = new TMultiGraph();
1731 0 : Int_t colArr[7] = {kRed, kRed, kOrange, kOrange, kGreen+1, kGreen+1, kGray+2};
1732 0 : const Int_t ngr = oaAlias->GetEntriesFast();
1733 0 : for (Int_t i=0; i<ngr; i++){
1734 0 : snprintf(query,200, "%s:%s", oaAlias->At(i)->GetName(), var_x);
1735 0 : multGr->Add( (TGraphErrors*) TStatToolkit::MakeGraphSparse(tree,query,cut,29,colArr[i],1.5) );
1736 : }
1737 : //
1738 0 : multGr->SetName(varname);
1739 0 : multGr->SetTitle(varname);
1740 0 : delete oaVar;
1741 0 : delete oaAlias;
1742 : return multGr;
1743 0 : }
1744 :
1745 :
1746 : TH1* TStatToolkit::DrawHistogram(TTree * tree, const char* drawCommand, const char* cuts, const char* histoname, const char* histotitle, Int_t nsigma, Float_t fraction, TObjArray *description )
1747 : {
1748 : //
1749 : // Draw histogram from TTree with robust range
1750 : // Only for 1D so far!
1751 : //
1752 : // Parameters:
1753 : // - histoname: name of histogram
1754 : // - histotitle: title of histgram
1755 : // - fraction: fraction of data to define the robust mean
1756 : // - nsigma: nsigma value for range
1757 : //
1758 : // To add:
1759 : // automatic ranges - separatelly for X, Y and Z nbins - as string
1760 : // names for the variables
1761 : // option, entries, first entry like in tree draw
1762 : //
1763 0 : TString drawStr(drawCommand);
1764 0 : TString cutStr(cuts);
1765 : Int_t dim = 1;
1766 :
1767 0 : if(!tree) {
1768 0 : ::Error("TStatToolkit::DrawHistogram","Tree pointer is NULL!");
1769 0 : return 0;
1770 : }
1771 :
1772 : // get entries
1773 0 : Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff");
1774 0 : if (entries == -1) {
1775 0 : ::Error("TStatToolkit::DrawHistogram","Tree draw returns -!");
1776 0 : return 0;
1777 : }
1778 0 : TObjArray *charray = drawStr.Tokenize(":");
1779 :
1780 : // get dimension
1781 0 : if(tree->GetV1()) dim = 1;
1782 0 : if(tree->GetV2()) dim = 2;
1783 0 : if(tree->GetV3()) dim = 3;
1784 0 : if(dim > 2){
1785 0 : cerr<<"TTree has more than 2 dimensions (not yet supported)"<<endl;
1786 0 : return 0;
1787 : }
1788 :
1789 : // draw robust
1790 : // Get estimators
1791 0 : Double_t mean1=0, rms1=0, min1=0, max1=0;
1792 0 : Double_t mean2=0, rms2=0, min2=0, max2=0;
1793 0 : Double_t mean3=0, rms3=0, min3=0, max3=0;
1794 :
1795 0 : TStatToolkit::GetMinMaxMean( tree->GetV1(),entries, min1,max1, mean1);
1796 0 : TStatToolkit::EvaluateUni(entries, tree->GetV1(),mean1,rms1, fraction*entries);
1797 0 : if(dim>1){
1798 0 : TStatToolkit::GetMinMaxMean( tree->GetV2(),entries, min2,max2, mean2);
1799 0 : TStatToolkit::EvaluateUni(entries, tree->GetV1(),mean2,rms2, fraction*entries);
1800 : }
1801 0 : if(dim>2){
1802 0 : TStatToolkit::GetMinMaxMean( tree->GetV3(),entries, min3,max3, mean3);
1803 0 : TStatToolkit::EvaluateUni(entries, tree->GetV3(),mean3,rms3, fraction*entries);
1804 : }
1805 :
1806 : TH1* hOut=NULL;
1807 0 : if(dim==1){
1808 0 : hOut = new TH1F(histoname, histotitle, 200, mean1-nsigma*rms1, mean1+nsigma*rms1);
1809 0 : for (Int_t i=0; i<entries; i++) hOut->Fill(tree->GetV1()[i]);
1810 0 : hOut->GetXaxis()->SetTitle(tree->GetHistogram()->GetXaxis()->GetTitle());
1811 : }
1812 0 : else if(dim==2){
1813 0 : hOut = new TH2F(histoname, histotitle, 200, min2, max2,200, mean1-nsigma*rms1, mean1+nsigma*rms1);
1814 0 : for (Int_t i=0; i<entries; i++) hOut->Fill(tree->GetV2()[i],tree->GetV1()[i]);
1815 0 : hOut->GetXaxis()->SetTitle(tree->GetHistogram()->GetXaxis()->GetTitle());
1816 0 : hOut->GetYaxis()->SetTitle(tree->GetHistogram()->GetYaxis()->GetTitle());
1817 : }
1818 0 : THashList * metaData = (THashList*) tree->GetUserInfo()->FindObject("metaTable");
1819 :
1820 0 : if (!metaData == 0){
1821 0 : TNamed *nmdTitle0 = TStatToolkit::GetMetadata(tree,Form("%s.Title",charray->At(0)->GetName()));
1822 0 : TNamed *nmdXAxis = TStatToolkit::GetMetadata(tree,Form("%s.AxisTitle",charray->At(1)->GetName()));
1823 0 : TNamed *nmdTitle1 = TStatToolkit::GetMetadata(tree,Form("%s.Title",charray->At(1)->GetName()));
1824 0 : TNamed *nmdYAxis = TStatToolkit::GetMetadata(tree,Form("%s.AxisTitle",charray->At(0)->GetName()));
1825 : //
1826 0 : TString hisTitle=charray->At(0)->GetName();
1827 0 : if (nmdTitle0) hisTitle=nmdTitle0->GetTitle();
1828 0 : if (nmdTitle1) {
1829 0 : hisTitle+=":";
1830 0 : hisTitle+=nmdTitle1->GetTitle();
1831 : }else{
1832 0 : hisTitle+=":";
1833 0 : hisTitle+=charray->At(1)->GetName();
1834 : }
1835 0 : if (nmdYAxis) {hOut->GetYaxis()->SetTitle(nmdYAxis->GetTitle());}
1836 0 : if (nmdXAxis) {hOut->GetXaxis()->SetTitle(nmdXAxis->GetTitle());}
1837 0 : hOut->SetTitle(hisTitle);
1838 0 : }
1839 0 : delete charray;
1840 : // if (option) hOut->Draw(option);
1841 : return hOut;
1842 0 : }
1843 :
1844 : void TStatToolkit::CheckTreeAliases(TTree * tree, Int_t ncheck){
1845 : //
1846 : // Check consistency of tree aliases
1847 : //
1848 : Int_t nCheck=100;
1849 0 : TList * aliases = (TList*)tree->GetListOfAliases();
1850 0 : Int_t entries = aliases->GetEntries();
1851 0 : for (Int_t i=0; i<entries; i++){
1852 0 : TObject * object= aliases->At(i);
1853 0 : if (!object) continue;
1854 0 : Int_t ndraw=tree->Draw(aliases->At(i)->GetName(),"1","goff",nCheck);
1855 0 : if (ndraw==0){
1856 0 : ::Error("Alias:\tProblem","%s",aliases->At(i)->GetName());
1857 0 : }else{
1858 0 : ::Info("Alias:\tOK","%s",aliases->At(i)->GetName());
1859 : }
1860 0 : }
1861 0 : }
1862 :
1863 :
1864 :
1865 :
1866 : Double_t TStatToolkit::GetDefaultStat(TTree * tree, const char * var, const char * selection, TStatType statType){
1867 : //
1868 : //
1869 : //
1870 0 : Int_t entries = tree->Draw(var,selection,"goff");
1871 0 : if (entries==0) return 0;
1872 0 : switch(statType){
1873 : case kEntries:
1874 0 : return entries;
1875 : case kSum:
1876 0 : return entries*TMath::Mean(entries, tree->GetV1());
1877 : case kMean:
1878 0 : return TMath::Mean(entries, tree->GetV1());
1879 : case kRMS:
1880 0 : return TMath::RMS(entries, tree->GetV1());
1881 : case kMedian:
1882 0 : return TMath::Median(entries, tree->GetV1());
1883 : }
1884 0 : return 0;
1885 0 : }
1886 :
1887 : //_____________________________________________________________________________
1888 : void TStatToolkit::CombineArray(TTree *tree, TVectorD &values)
1889 : {
1890 : /// Collect all variables from the last draw in one array.
1891 : ///
1892 : /// It is assumed that the Draw function of the TTree was called before
1893 : /// if e.g. Draw("v1:v2:v3") had been called, then values will contain
1894 : /// the concatenated array of the values from v1,v2 and v3.
1895 : /// E.g. if the v1[0..n], v2[0..n], v3[0..n] then
1896 : /// values[0..3n] = [v1, v2, v3]
1897 : /// \param[in] tree input tree
1898 : /// \param[out] values array in which to summarise all 'drawn' values
1899 0 : const Int_t numberOfDimensions = tree->GetPlayer()->GetDimension();
1900 0 : if (numberOfDimensions==1) {
1901 0 : values.Use(tree->GetSelectedRows(), tree->GetVal(0));
1902 0 : return;
1903 : }
1904 :
1905 0 : const Int_t numberOfSelectedRows = tree->GetSelectedRows();
1906 0 : values.ResizeTo(numberOfDimensions * numberOfSelectedRows);
1907 :
1908 : Int_t nfill=0;
1909 0 : for (Int_t idim=0; idim<numberOfDimensions; ++idim) {
1910 0 : const Double_t *arr = tree->GetVal(idim);
1911 0 : if (!arr) continue;
1912 :
1913 0 : for (Int_t ival=0; ival<numberOfSelectedRows; ++ival) {
1914 0 : values.GetMatrixArray()[nfill++] = arr[ival];
1915 : }
1916 0 : }
1917 :
1918 0 : }
1919 :
1920 : //_____________________________________________________________________________
1921 : Double_t TStatToolkit::GetDistance(const TVectorD &values, const ENormType normType,
1922 : const Bool_t normaliseToEntries/*=kFALSE*/, const Double_t pvalue/*=1*/)
1923 : {
1924 : /// Calculate the distance of the elements in values using a certain norm
1925 : /// \param[in] values array with input values
1926 : /// \param[in] normType normalisation to use
1927 : /// \param[in] normaliseToEntries divide the norm by the number of eleements ('average norm')
1928 : /// \param[in] pvalue the p value for the p-type norm, ignored for all other norms
1929 : /// \return calculated distance
1930 :
1931 : Double_t norm=0.;
1932 :
1933 0 : switch (normType) {
1934 : case kL1:
1935 0 : norm=values.Norm1();
1936 0 : break;
1937 : case kL2:
1938 0 : norm=TMath::Sqrt(values.Norm2Sqr());
1939 0 : break;
1940 : case kLp:
1941 : {
1942 0 : if (pvalue<1.) {
1943 0 : ::Error("TStatToolkit::GetDistance","Lp norm: p-value=%5.3g not valid. Only p-value>=1 is allowed", pvalue);
1944 0 : break;
1945 : }
1946 : Double_t sum=0.;
1947 0 : for (Int_t ival=0; ival<values.GetNrows(); ++ival) {
1948 0 : sum+=TMath::Power(TMath::Abs(values.GetMatrixArray()[ival]), pvalue);
1949 : }
1950 0 : norm=TMath::Power(sum, 1./pvalue);
1951 : }
1952 0 : break;
1953 : case kMax:
1954 0 : norm=values.NormInf();
1955 0 : break;
1956 : case kHamming:
1957 : {
1958 : Double_t sum=0.;
1959 0 : for (Int_t ival=0; ival<values.GetNrows(); ++ival) {
1960 0 : if (TMath::Abs(values.GetMatrixArray()[ival])>1e-30) ++sum;
1961 : }
1962 : norm=sum;
1963 : }
1964 0 : break;
1965 : }
1966 0 : if (normaliseToEntries && values.GetNrows()>0) {
1967 0 : norm/=values.GetNrows();
1968 0 : }
1969 0 : return norm;
1970 : }
1971 :
1972 : //_____________________________________________________________________________
1973 : Double_t TStatToolkit::GetDistance(const Int_t size, const Double_t *values, const ENormType normType,
1974 : const Bool_t normaliseToEntries/*=kFALSE*/, const Double_t pvalue/*=1*/)
1975 : {
1976 : /// Calculate the distance of the elements in values using a certain norm
1977 : /// \sa GetDistance()
1978 0 : TVectorD vecvalues;
1979 0 : vecvalues.Use(size, values);
1980 0 : return GetDistance(vecvalues, normType, normaliseToEntries, pvalue);
1981 0 : }
1982 :
1983 : //_____________________________________________________________________________
1984 : Double_t TStatToolkit::GetDistance(TTree * tree, const char* var, const char * selection,
1985 : const ENormType normType, const Bool_t normaliseToEntries/*=kFALSE*/, const Double_t pvalue/*=1*/)
1986 : {
1987 : /// Calculate the distance of the values selecte in tree->Draw(var, selection)
1988 : ///
1989 : /// If var contains more than one variable (separated by ':' as usual) the arrays
1990 : /// are concatenated:<BR>
1991 : /// E.g. if var="v1:v2:v3", then the norm of the
1992 : /// the concatenated array of the values from v1,v2 and v3 will be calculated:<BR>
1993 : /// This means if the internal tree arrays for each variable are v1[0..n], v2[0..n], v3[0..n] then
1994 : /// the norm of vx[0..3n] = [v1, v2, v3] is calculated.
1995 : /// \param[in] tree input tree
1996 : /// \param[in] var variable expression for the tree->Draw()
1997 : /// \param[in] selection selection for the tree->Draw()
1998 : /// \param[in] normType norm to use for calculating the point distances
1999 : /// \param[in] normaliseToEntries divide the norm by the number of eleements ('average norm')
2000 : /// \param[in] pvalue p-value for the p-norm (ignored for other norm types
2001 : /// \return calculated distnace
2002 0 : Int_t entries = tree->Draw(var,selection,"goff");
2003 0 : if (entries==0) return 0.;
2004 :
2005 0 : TVectorD values;
2006 0 : CombineArray(tree, values);
2007 0 : return GetDistance(values, normType, normaliseToEntries, pvalue);
2008 0 : }
2009 :
2010 :
2011 :
2012 : void TStatToolkit::MakeDistortionMap(Int_t iter, THnBase * histo, TTreeSRedirector *pcstream, TMatrixD &projectionInfo,Int_t dumpHisto, Int_t verbose){
2013 : //
2014 : // Recursive function to calculate Distortion maps from the residual histograms
2015 : // Input:
2016 : // iter - ndim..0
2017 : // histo - THn histogram
2018 : // pcstream -
2019 : // projectionInfo - TMatrix speicifiing distortion map cration setup
2020 : // user specify columns:
2021 : // 0.) sequence of dimensions
2022 : // 1.) grouping in dimensions (how many bins will be groupd in specific dimension - 0 means onl specified bin 1, curren +-1 bin ...)
2023 : // 2.) step in dimension ( in case >1 some n(projectionInfo(<dim>,2) bins will be not exported
2024 : // internally used collumns (needed to pass current bin index and bin center to the recursive function)
2025 : // 3.) current bin value
2026 : // 4.) current bin center
2027 : //
2028 : // Output:
2029 : // pcstream - file with output distortion tree
2030 : // 1.) distortion characteristic: mean, rms, gaussian fit parameters, meang, rmsG chi2 ... at speciefied bin
2031 : // 2.) specidfied bins (tree branches) are defined by the name of the histogram axis in input histograms
2032 : //
2033 : //
2034 : // Example projection info
2035 : /*
2036 : TFile *f = TFile::Open("/hera/alice/hellbaer/alice-tpc-notes/SpaceChargeDistortion/data/ATO-108/fullMerge/SCcalibMergeLHC12d.root");
2037 : THnF* histof= (THnF*) f->Get("deltaY_ClTPC_ITSTOF");
2038 : histof->SetName("deltaRPhiTPCTISTOF");
2039 : histof->GetAxis(4)->SetName("qpt");
2040 : TH1::SetDirectory(0);
2041 : TTreeSRedirector * pcstream = new TTreeSRedirector("distortion.root","recreate");
2042 : TMatrixD projectionInfo(5,5);
2043 : projectionInfo(0,0)=0; projectionInfo(0,1)=0; projectionInfo(0,2)=0;
2044 : projectionInfo(1,0)=4; projectionInfo(1,1)=0; projectionInfo(1,2)=1;
2045 : projectionInfo(2,0)=3; projectionInfo(2,1)=3; projectionInfo(2,2)=2;
2046 : projectionInfo(3,0)=2; projectionInfo(3,1)=0; projectionInfo(3,2)=5;
2047 : projectionInfo(4,0)=1; projectionInfo(4,1)=5; projectionInfo(4,2)=20;
2048 : MakeDistortionMap(4, histof, pcstream, projectionInfo);
2049 : delete pcstream;
2050 : */
2051 : //
2052 0 : static TF1 fgaus("fgaus","gaus",-10,10);
2053 : const Double_t kMinEntries=50;
2054 0 : Int_t ndim=histo->GetNdimensions();
2055 0 : Int_t axis[ndim];
2056 0 : Double_t meanVector[ndim];
2057 0 : Int_t binVector[ndim];
2058 0 : Double_t centerVector[ndim];
2059 0 : for (Int_t idim=0; idim<ndim; idim++) axis[idim]=idim;
2060 : //
2061 0 : if (iter==0){
2062 0 : char tname[100];
2063 0 : char aname[100];
2064 0 : char bname[100];
2065 0 : char cname[100];
2066 :
2067 0 : snprintf(tname, 100, "%sDist",histo->GetName());
2068 : //
2069 : //
2070 : // 1.) Calculate properties - mean, rms, gaus fit, chi2, entries
2071 : // 2.) Dump properties to tree 1D properties - plus dimension descriptor f
2072 : Int_t axis1D[1]={0};
2073 0 : Int_t dimProject = TMath::Nint(projectionInfo(iter,0));
2074 : axis1D[0]=dimProject;
2075 0 : TH1 *his1DFull = histo->Projection(dimProject);
2076 0 : Double_t mean= his1DFull->GetMean();
2077 0 : Double_t rms= his1DFull->GetRMS();
2078 0 : Int_t entries= his1DFull->GetEntries();
2079 0 : TString hname="his_";
2080 0 : for (Int_t idim=0; idim<ndim; idim++) {hname+="_"; hname+=TMath::Nint(projectionInfo(idim,3));}
2081 0 : Double_t meanG=0, rmsG=0, chi2G=0;
2082 0 : if (entries>kMinEntries){
2083 0 : fgaus.SetParameters(entries,mean,rms);
2084 0 : his1DFull->Fit(&fgaus,"qnr","qnr");
2085 0 : meanG = fgaus.GetParameter(1);
2086 0 : rmsG = fgaus.GetParameter(2);
2087 0 : chi2G = fgaus.GetChisquare()/fgaus.GetNumberFreeParameters();
2088 0 : }
2089 0 : if (dumpHisto>=0) {
2090 : static Int_t histoCounter=0;
2091 0 : if ((histoCounter%dumpHisto)==0) his1DFull->Write(hname.Data());
2092 0 : histoCounter++;
2093 0 : }
2094 0 : delete his1DFull;
2095 0 : (*pcstream)<<tname<<
2096 0 : "entries="<<entries<< // number of entries
2097 0 : "mean="<<mean<< // mean value of the last dimension
2098 0 : "rms="<<rms<< // rms value of the last dimension
2099 0 : "meanG="<<meanG<< // mean of the gaus fit
2100 0 : "rmsG="<<rmsG<< // rms of the gaus fit
2101 0 : "chi2G="<<chi2G; // chi2 of the gaus fit
2102 :
2103 0 : for (Int_t idim=0; idim<ndim; idim++){
2104 : axis1D[0]=idim;
2105 0 : TH1 *his1DAxis = histo->Projection(idim);
2106 0 : meanVector[idim] = his1DAxis->GetMean();
2107 0 : snprintf(aname, 100, "%sMean=",histo->GetAxis(idim)->GetName());
2108 0 : (*pcstream)<<tname<<
2109 0 : aname<<meanVector[idim]; // current bin means
2110 0 : delete his1DAxis;
2111 : }
2112 0 : for (Int_t iIter=0; iIter<ndim; iIter++){
2113 0 : Int_t idim = TMath::Nint(projectionInfo(iIter,0));
2114 0 : binVector[idim] = TMath::Nint(projectionInfo(iIter,3));
2115 0 : centerVector[idim] = projectionInfo(iIter,4);
2116 0 : snprintf(bname, 100, "%sBin=",histo->GetAxis(idim)->GetName());
2117 0 : snprintf(cname, 100, "%sCenter=",histo->GetAxis(idim)->GetName());
2118 0 : (*pcstream)<<tname<<
2119 0 : bname<<binVector[idim]<< // current bin values
2120 0 : cname<<centerVector[idim]; // current bin centers
2121 : }
2122 0 : (*pcstream)<<tname<<"\n";
2123 0 : }else{
2124 : // loop over the diminsion of interest
2125 : // project selecting bin+-deltabin histoProj
2126 : // MakeDistortionMap(histoProj ...)
2127 : //
2128 0 : Int_t dimProject = TMath::Nint(projectionInfo(iter,0));
2129 0 : Int_t groupProject = TMath::Nint(projectionInfo(iter,1));
2130 0 : Int_t stepProject = TMath::Nint(projectionInfo(iter,2));
2131 0 : if (stepProject<1) stepProject=1;
2132 0 : Int_t nbins = histo->GetAxis(dimProject)->GetNbins();
2133 :
2134 0 : for (Int_t ibin=1; ibin<=nbins; ibin+=stepProject){
2135 0 : if (iter>1 && verbose){
2136 0 : for (Int_t idim=0; idim<ndim; idim++){
2137 0 : printf("\t%d(%d,%d)",TMath::Nint(projectionInfo(idim,3)),TMath::Nint(projectionInfo(idim,0)),TMath::Nint(projectionInfo(idim,1) ));
2138 : }
2139 0 : printf("\n");
2140 0 : AliSysInfo::AddStamp("xxx",iter, dimProject);
2141 0 : }
2142 0 : Int_t bin0=TMath::Max(ibin-groupProject,1);
2143 0 : Int_t bin1=TMath::Min(ibin+groupProject,nbins);
2144 0 : histo->GetAxis(dimProject)->SetRange(bin0,bin1);
2145 0 : projectionInfo(iter,3)=ibin;
2146 0 : projectionInfo(iter,4)=histo->GetAxis(dimProject)->GetBinCenter(ibin);
2147 0 : Int_t iterProject=iter-1;
2148 0 : MakeDistortionMap(iterProject, histo, pcstream, projectionInfo);
2149 : }
2150 : }
2151 : //
2152 0 : }
2153 :
2154 : void TStatToolkit::MakeDistortionMapFast(THnBase * histo, TTreeSRedirector *pcstream, TMatrixD &projectionInfo,Int_t verbose, Double_t fractionCut, const char * estimators)
2155 : {
2156 : //
2157 : // Function to calculate Distortion maps from the residual histograms
2158 : // Input:
2159 : // histo - THn histogram
2160 : // pcstream -
2161 : // projectionInfo - TMatrix speicifiing distortion map cration setup
2162 : // user specify columns:
2163 : // 0.) sequence of dimensions
2164 : // 1.) grouping in dimensions (how many bins will be groupd in specific dimension - 0 means onl specified bin 1, curren +-1 bin ...)
2165 : // 2.) step in dimension ( in case >1 some n(projectionInfo(<dim>,2) bins will be not exported
2166 : //
2167 : // Output:
2168 : // pcstream - file with output distortion tree
2169 : // 1.) distortion characteristic: mean, rms, gaussian fit parameters, meang, rmsG chi2 ... at speciefied bin
2170 : // 2.) specidfied bins (tree branches) are defined by the name of the histogram axis in input histograms
2171 : // 3.) in debug mode - controlled by env variable "gDumpHistoFraction" fractio of histogram + fits dumped to the file
2172 : //
2173 : // Example projection info
2174 : /*
2175 : TFile *f = TFile::Open("/hera/alice/hellbaer/alice-tpc-notes/SpaceChargeDistortion/data/ATO-108/fullMerge/SCcalibMergeLHC12d.root");
2176 : THnF* histof= (THnF*) f->Get("deltaY_ClTPC_ITSTOF");
2177 : histof->SetName("deltaRPhiTPCTISTOF");
2178 : histof->GetAxis(4)->SetName("qpt");
2179 : TH1::SetDirectory(0);
2180 : TTreeSRedirector * pcstream = new TTreeSRedirector("distortion.root","recreate");
2181 : TMatrixD projectionInfo(5,3);
2182 : projectionInfo(0,0)=0; projectionInfo(0,1)=0; projectionInfo(0,2)=0;
2183 : projectionInfo(1,0)=4; projectionInfo(1,1)=0; projectionInfo(1,2)=1;
2184 : projectionInfo(2,0)=3; projectionInfo(2,1)=3; projectionInfo(2,2)=2;
2185 : projectionInfo(3,0)=2; projectionInfo(3,1)=0; projectionInfo(3,2)=5;
2186 : projectionInfo(4,0)=1; projectionInfo(4,1)=5; projectionInfo(4,2)=20;
2187 : MakeDistortionMap(histof, pcstream, projectionInfo);
2188 : delete pcstream;
2189 : */
2190 : //
2191 : const Double_t kMinEntries=30, kUseLLFrom=20;
2192 0 : const Float_t kDumpHistoFraction = TString(gSystem->Getenv("gDumpHistoFraction")).Atof(); // in debug mode - controlled by env variable "gDumpHistoFraction" fractio of histogram + fits dumped to the file
2193 0 : char tname[100];
2194 0 : char aname[100];
2195 0 : char bname[100];
2196 0 : char cname[100];
2197 0 : Float_t fractionLTM[100]={0.8};
2198 0 : TVectorF *vecLTM[100]={0};
2199 : Int_t nestimators=1;
2200 0 : if (estimators!=NULL){
2201 0 : TObjArray * array=TString(estimators).Tokenize(":");
2202 0 : nestimators=array->GetEntries();
2203 0 : for (Int_t iest=0; iest<nestimators; iest++){
2204 0 : fractionLTM[iest]=TString(array->At(iest)->GetName()).Atof();
2205 : }
2206 0 : }
2207 0 : for (Int_t iest=0; iest<nestimators; iest++) {
2208 0 : vecLTM[iest]=new TVectorF(10);
2209 0 : (*(vecLTM[iest]))[9]= fractionLTM[iest];
2210 : }
2211 :
2212 : //
2213 0 : int ndim = histo->GetNdimensions();
2214 0 : int nbins[ndim],idx[ndim],idxmin[ndim],idxmax[ndim],idxSav[ndim];
2215 0 : for (int id=0;id<ndim;id++) nbins[id] = histo->GetAxis(id)->GetNbins();
2216 : //
2217 0 : int axOrd[ndim],binSt[ndim],binGr[ndim];
2218 0 : for (int i=0;i<ndim;i++) {
2219 0 : axOrd[i] = TMath::Nint(projectionInfo(i,0));
2220 0 : binGr[i] = TMath::Nint(projectionInfo(i,1));
2221 0 : binSt[i] = TMath::Max(1,TMath::Nint(projectionInfo(i,2)));
2222 : }
2223 0 : int tgtDim = axOrd[0],tgtStep=binSt[0],tgtNb=nbins[tgtDim],tgtNb1=tgtNb+1;
2224 0 : double binY[tgtNb],binX[tgtNb],meanVector[ndim],centerVector[ndim];
2225 0 : Int_t binVector[ndim];
2226 : // prepare X axis
2227 0 : TAxis* xax = histo->GetAxis(tgtDim);
2228 0 : for (int i=tgtNb;i--;) binX[i] = xax->GetBinCenter(i+1);
2229 0 : for (int i=ndim;i--;) idx[i]=1;
2230 : Bool_t grpOn = kFALSE;
2231 0 : for (int i=1;i<ndim;i++) if (binGr[i]) grpOn = kTRUE;
2232 : //
2233 : // estimate number of output fits
2234 0 : histo->GetListOfAxes()->Print();
2235 : ULong64_t nfits = 1, fitCount=0;
2236 0 : printf("index\tdim\t|\tnbins\tgrouping\tstep\tnfits\n");
2237 0 : for (int i=1;i<ndim;i++) {
2238 0 : int idim = axOrd[i];
2239 0 : nfits *= TMath::Max(1,nbins[idim]/binSt[idim]);
2240 0 : printf("%d %d | %d %d %d %lld\n",i,idim,nbins[idim],binGr[idim], binSt[idim],nfits);
2241 : }
2242 0 : printf("Expect %lld nfits\n",nfits);
2243 0 : ULong64_t fitProgress = nfits/100;
2244 : //
2245 : // setup fit function, at the moment full root fit
2246 0 : static TF1 fgaus("fgaus","gaus",-10,10);
2247 0 : fgaus.SetRange(xax->GetXmin(),xax->GetXmax());
2248 : // TGraph grafFit(tgtNb);
2249 0 : TH1F* hfit = new TH1F("hfit","hfit",tgtNb,xax->GetXmin(),xax->GetXmax());
2250 : //
2251 0 : snprintf(tname, 100, "%sDist",histo->GetName());
2252 0 : TStopwatch sw;
2253 0 : sw.Start();
2254 0 : int dimVar=1, dimVarID = axOrd[dimVar];
2255 : //
2256 : // TVectorF vecLTM(9);
2257 0 : while(1) {
2258 : //
2259 0 : double dimVarCen = histo->GetAxis(dimVarID)->GetBinCenter(idx[dimVarID]); // center of currently varied bin
2260 : //
2261 0 : if (grpOn) { //>> grouping requested?
2262 0 : memset(binY,0,tgtNb*sizeof(double)); // need to accumulate
2263 : //
2264 0 : for (int idim=1;idim<ndim;idim++) {
2265 0 : int grp = binGr[idim];
2266 0 : int idimR = axOrd[idim]; // real axis id
2267 0 : idxSav[idimR]=idx[idimR]; // save central bins
2268 0 : idxmax[idimR] = TMath::Min(idx[idimR]+grp,nbins[idimR]);
2269 0 : idx[idimR] = idxmin[idimR] = TMath::Max(1,idx[idimR]-grp);
2270 : //
2271 : // effective bin center
2272 0 : meanVector[idimR] = 0;
2273 0 : TAxis* ax = histo->GetAxis(idimR);
2274 0 : if (grp>0) {
2275 0 : for (int k=idxmin[idimR];k<=idxmax[idimR];k++) meanVector[idimR] += ax->GetBinCenter(k);
2276 0 : meanVector[idimR] /= (1+(grp<<1));
2277 0 : }
2278 0 : else meanVector[idimR] = ax->GetBinCenter(idxSav[idimR]);
2279 : } // set limits for grouping
2280 0 : if (verbose>0) {
2281 0 : printf("output bin: ");
2282 0 : for (int i=0;i<ndim;i++) if (i!=tgtDim) printf("[D:%d]:%d ",i,idxSav[i]); printf("\n");
2283 0 : printf("integrates: ");
2284 0 : for (int i=0;i<ndim;i++) if (i!=tgtDim) printf("[D:%d]:%d-%d ",i,idxmin[i],idxmax[i]); printf("\n");
2285 : }
2286 : //
2287 : while(1) {
2288 : // loop over target dimension: accumulation
2289 0 : int &it = idx[tgtDim];
2290 0 : for (it=1;it<tgtNb1;it+=tgtStep) {
2291 0 : binY[it-1] += histo->GetBinContent(idx);
2292 0 : if (verbose>1) {for (int i=0;i<ndim;i++) printf("%d ",idx[i]); printf(" | accumulation\n");}
2293 : }
2294 : //
2295 : int idim;
2296 0 : for (idim=1;idim<ndim;idim++) { // dimension being groupped
2297 0 : int idimR = axOrd[idim]; // real axis id in the histo
2298 0 : if ( (++idx[idimR]) > idxmax[idimR] ) idx[idimR]=idxmin[idimR];
2299 0 : else break;
2300 0 : }
2301 0 : if (idim==ndim) break;
2302 0 : }
2303 : } // <<grouping requested
2304 : else {
2305 0 : int &it = idx[tgtDim];
2306 0 : for (it=1;it<tgtNb1;it+=tgtStep) {
2307 0 : binY[it-1] = histo->GetBinContent(idx);
2308 : //
2309 : //for (int i=0;i<ndim;i++) printf("%d ",idx[i]); printf(" | \n");
2310 : }
2311 0 : for (int idim=1;idim<ndim;idim++) {
2312 0 : int idimR = axOrd[idim]; // real axis id
2313 0 : meanVector[idimR] = histo->GetAxis(idimR)->GetBinCenter(idx[idimR]);
2314 : }
2315 : }
2316 0 : if (grpOn) for (int i=ndim;i--;) idx[i]=idxSav[i]; // restore central bins
2317 0 : idx[tgtDim] = 0;
2318 0 : if (verbose>0) {for (int i=0;i<ndim;i++) printf("%d ",idx[i]); printf(" | central bin fit\n");}
2319 : //
2320 : // >> ------------- do fit
2321 0 : float mean=0,mom2=0,rms=0,m3=0, m4=0, nrm=0,meanG=0,rmsG=0,chi2G=0,maxVal=0,entriesG=0,mean0=0, rms0=0, curt0=0;
2322 0 : hfit->Reset();
2323 0 : for (int ip=tgtNb;ip--;) {
2324 : //grafFit.SetPoint(ip,binX[ip],binY[ip]);
2325 0 : hfit->SetBinContent(ip+1,binY[ip]);
2326 0 : nrm += binY[ip];
2327 0 : mean += binX[ip]*binY[ip];
2328 0 : mom2 += binX[ip]*binX[ip]*binY[ip];
2329 0 : if (maxVal<binY[ip]) maxVal = binY[ip];
2330 : }
2331 0 : if (nrm>0) {
2332 0 : mean /= nrm;
2333 0 : mom2 /= nrm;
2334 0 : rms = mom2 - mean*mean;
2335 0 : rms = rms>0 ? TMath::Sqrt(rms):0;
2336 0 : }
2337 0 : mean0=mean;
2338 0 : rms0=rms;
2339 :
2340 :
2341 0 : Int_t nbins1D=hfit->GetNbinsX();
2342 0 : Float_t binMedian=0;
2343 0 : Double_t limits[2]={hfit->GetBinCenter(1), hfit->GetBinCenter(nbins1D)};
2344 0 : if (nrm>5) {
2345 0 : for (Int_t iest=0; iest<nestimators; iest++){
2346 0 : TStatToolkit::LTMHisto(hfit, *(vecLTM[iest]), fractionLTM[iest]);
2347 : }
2348 0 : Double_t* integral=hfit->GetIntegral();
2349 0 : for (Int_t i=1; i<nbins1D-1; i++){
2350 0 : if (integral[i-1]<0.5 && integral[i]>=0.5){
2351 0 : if (hfit->GetBinContent(i-1)+hfit->GetBinContent(i)>0){
2352 0 : binMedian=hfit->GetBinCenter(i);
2353 0 : Double_t dIdx=-(integral[i-1]-integral[i]);
2354 0 : Double_t dx=(0.5+(0.5-integral[i])/dIdx)*hfit->GetBinWidth(i);
2355 0 : binMedian+=dx;
2356 0 : }
2357 : }
2358 0 : if (integral[i-1]<fractionCut && integral[i]>=fractionCut){
2359 0 : limits[0]=hfit->GetBinCenter(i-1)-hfit->GetBinWidth(i);
2360 0 : }
2361 0 : if (integral[i]<1-fractionCut && integral[i+1]>=1-fractionCut){
2362 0 : limits[1]=hfit->GetBinCenter(i+1)+hfit->GetBinWidth(i);
2363 0 : }
2364 : }
2365 0 : }
2366 0 : if (nrm>5&&fractionCut>0 &&rms>0) {
2367 0 : hfit->GetXaxis()->SetRangeUser(limits[0], limits[1]);
2368 0 : mean=hfit->GetMean();
2369 0 : rms=hfit->GetRMS();
2370 0 : if (nrm>0 && rms>0) {
2371 0 : m3=hfit->GetSkewness();
2372 0 : m4=hfit->GetKurtosis();
2373 0 : }
2374 0 : fgaus.SetRange(limits[0]-rms, limits[1]+rms);
2375 : }else{
2376 0 : fgaus.SetRange(xax->GetXmin(),xax->GetXmax());
2377 : }
2378 :
2379 :
2380 0 : Bool_t isFitValid=kFALSE;
2381 0 : if (nrm>=kMinEntries && rms>0) {
2382 0 : fgaus.SetParameters(nrm/(rms/hfit->GetBinWidth(nbins1D)),mean,rms);
2383 : //grafFit.Fit(&fgaus,/*maxVal<kUseLLFrom ? "qnrl":*/"qnr");
2384 0 : TFitResultPtr fitPtr= hfit->Fit(&fgaus,maxVal<kUseLLFrom ? "qnrlS":"qnrS");
2385 0 : entriesG = fgaus.GetParameter(0);
2386 0 : meanG = fgaus.GetParameter(1);
2387 0 : rmsG = fgaus.GetParameter(2);
2388 0 : chi2G = fgaus.GetChisquare()/fgaus.GetNumberFreeParameters();
2389 0 : TFitResult * result = fitPtr.Get();
2390 0 : if (result!=NULL){
2391 0 : isFitValid = result->IsValid();
2392 0 : }
2393 : //
2394 0 : }
2395 : TH1 * hDump=0;
2396 0 : if (nrm>=kMinEntries&& kDumpHistoFraction>0 && (gRandom->Rndm()<kDumpHistoFraction || isFitValid!=kTRUE)){
2397 : hDump=hfit;
2398 0 : }
2399 0 : if (hDump){
2400 0 : (*pcstream)<<TString::Format("%sDump", tname).Data()<<
2401 0 : "entries="<<nrm<< // number of entries
2402 0 : "isFitValid="<<isFitValid<< // true if the gaus fit converged
2403 0 : "hDump.="<<hDump<< // histogram - by default not filled
2404 0 : "mean0="<<mean0<< // mean value of the last dimension - without fraction cut
2405 0 : "rms0="<<rms0<< // rms value of the last dimension - without fraction cut
2406 0 : "mean="<<mean<< // mean value of the last dimension
2407 0 : "rms="<<rms<< // rms value of the last dimension
2408 0 : "m3="<<m3<< // m3 (skewnes) of the last dimension
2409 0 : "m4="<<m4<< // m4 (kurtosis) of the last dimension
2410 0 : "binMedian="<<binMedian<< //binned median value of 1D histogram
2411 0 : "entriesG="<<entriesG<<
2412 0 : "meanG="<<meanG<< // mean of the gaus fit
2413 0 : "rmsG="<<rmsG<< // rms of the gaus fit
2414 0 : "vecLTM.="<<vecLTM[0]<< // LTM frac% statistic
2415 0 : "chi2G="<<chi2G; // chi2 of the gaus fit
2416 0 : for (Int_t iest=1; iest<nestimators; iest++)
2417 0 : (*pcstream)<<TString::Format("%sDump", tname).Data()<<TString::Format("vecLTM%d.=",iest)<<vecLTM[iest]; // LTM frac% statistic
2418 :
2419 0 : }
2420 :
2421 0 : (*pcstream)<<tname<<
2422 0 : "entries="<<nrm<< // number of entries
2423 0 : "isFitValid="<<isFitValid<< // true if the gaus fit converged
2424 0 : "mean0="<<mean0<< // mean value of the last dimension - without fraction cut
2425 0 : "rms0="<<rms0<< // rms value of the last dimension - without fraction cut
2426 0 : "mean="<<mean<< // mean value of the last dimension
2427 0 : "rms="<<rms<< // rms value of the last dimension
2428 0 : "m3="<<m3<< // m3 (skewnes) of the last dimension
2429 0 : "m4="<<m4<< // m4 (kurtosis) of the last dimension
2430 0 : "binMedian="<<binMedian<< //binned median value of 1D histogram
2431 0 : "entriesG="<<entriesG<< //
2432 0 : "meanG="<<meanG<< // mean of the gaus fit
2433 0 : "rmsG="<<rmsG<< // rms of the gaus fit
2434 0 : "vecLTM.="<<vecLTM[0]<< // LTM frac% statistic
2435 0 : "chi2G="<<chi2G; // chi2 of the gaus fit
2436 0 : for (Int_t iest=1; iest<nestimators; iest++)
2437 0 : (*pcstream)<<tname<<TString::Format("vecLTM%d.=",iest)<<vecLTM[iest]; // LTM frac% statistic
2438 :
2439 :
2440 : //
2441 0 : meanVector[tgtDim] = mean; // what's a point of this?
2442 0 : for (Int_t idim=0; idim<ndim; idim++){
2443 0 : snprintf(aname, 100, "%sMean=",histo->GetAxis(idim)->GetName());
2444 0 : (*pcstream)<<tname<<
2445 0 : aname<<meanVector[idim]; // current bin means
2446 : }
2447 : //
2448 0 : for (Int_t iIter=0; iIter<ndim; iIter++){
2449 0 : Int_t idim = axOrd[iIter];
2450 0 : binVector[idim] = idx[idim];
2451 0 : centerVector[idim] = histo->GetAxis(idim)->GetBinCenter(idx[idim]);
2452 0 : snprintf(bname, 100, "%sBin=",histo->GetAxis(idim)->GetName());
2453 0 : snprintf(cname, 100, "%sCenter=",histo->GetAxis(idim)->GetName());
2454 0 : (*pcstream)<<tname<<
2455 0 : bname<<binVector[idim]<< // current bin values
2456 0 : cname<<centerVector[idim]; // current bin centers
2457 0 : if (hDump){
2458 0 : (*pcstream)<<TString::Format("%sDump", tname).Data()<<
2459 0 : bname<<binVector[idim]<< // current bin values
2460 0 : cname<<centerVector[idim]; // current bin centers
2461 0 : }
2462 : }
2463 0 : (*pcstream)<<tname<<"\n";
2464 0 : if (hDump) (*pcstream)<<TString::Format("%sDump", tname).Data()<<"\n";
2465 : // << ------------- do fit
2466 : //
2467 0 : if (((++fitCount)%fitProgress)==0) {
2468 0 : printf("fit %lld %4.1f%% done\n",fitCount,100*double(fitCount)/nfits);
2469 0 : AliSysInfo::AddStamp("fitCout", 1,fitCount,100*double(fitCount)/nfits);
2470 : }
2471 : //
2472 : //next global bin in which target dimention will be looped
2473 0 : for (dimVar=1;dimVar<ndim;dimVar++) { // varying dimension
2474 0 : dimVarID = axOrd[dimVar]; // real axis id in the histo
2475 0 : if ( (idx[dimVarID]+=binSt[dimVar]) > nbins[dimVarID] ) idx[dimVarID]=1;
2476 : else break;
2477 : }
2478 0 : if (dimVar==ndim) break;
2479 0 : }
2480 0 : delete hfit;
2481 0 : sw.Stop();
2482 0 : sw.Print();
2483 : /*
2484 : int nb = histo->GetNbins();
2485 : int prc = nb/100;
2486 : for (int i=0;i<nb;i++) {
2487 : histo->GetBinContent(i);
2488 : if (i && (i%prc)==0) printf("Done %d%%\n",int(float(100*i)/prc));
2489 : }
2490 : */
2491 0 : }
|
