AliRoot Core  3dc7879 (3dc7879)
TStatToolkit.cxx
Go to the documentation of this file.
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 //-----------------------------------------------------------------------------
20 //-----------------------------------------------------------------------------
21 
22 //
24 #include "TStopwatch.h"
25 #include "TStatToolkit.h"
26 #include "TTreeFormula.h"
27 #include "TLegend.h"
28 #include "TPRegexp.h"
29 #include "AliDrawStyle.h"
30 #include <stdlib.h>
31 #include <stdexcept> // std::invalid_argument
32 using std::cout;
33 using std::cerr;
34 using std::endl;
35 using namespace std;
36 
37 //_____________________________________________________________________________
38 void TStatToolkit::EvaluateUni(Int_t nvectors, Double_t *data, Double_t &mean, Double_t &sigma, Int_t hh) {
39  //
40  // Robust estimator in 1D case MI version - (faster than ROOT version)
41  //
42  // For the univariate case
43  // estimates of location and scatter are returned in mean and sigma parameters
44  // the algorithm works on the same principle as in multivariate case -
45  // it finds a subset of size hh with smallest sigma, and then returns mean and
46  // sigma of this subset
47  //
48 
49  if (hh==0)
50  hh=(nvectors+2)/2;
51  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};
52  Int_t *index=new Int_t[nvectors];
53  TMath::Sort(nvectors, data, index, kFALSE);
54 
55  Int_t nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11);
56  Double_t factor = faclts[TMath::Max(0,nquant-1)];
57 
58  Double_t sumx =0;
59  Double_t sumx2 =0;
60  Double_t bestmean = 0;
61  Double_t bestsigma = (data[index[nvectors-1]]-data[index[0]]+1.); // maximal possible sigma
62  bestsigma *=bestsigma;
63 
64  for (Int_t i=0; i<hh; i++){
65  sumx += data[index[i]];
66  sumx2 += data[index[i]]*data[index[i]];
67  }
68 
69  Double_t norm = 1./Double_t(hh);
70  Double_t norm2 = (hh-1)>0 ? 1./Double_t(hh-1):1;
71  for (Int_t i=hh; i<nvectors; i++){
72  Double_t cmean = sumx*norm;
73  Double_t csigma = (sumx2 - hh*cmean*cmean)*norm2;
74  if (csigma<bestsigma){
75  bestmean = cmean;
76  bestsigma = csigma;
77  }
78 
79  sumx += data[index[i]]-data[index[i-hh]];
80  sumx2 += data[index[i]]*data[index[i]]-data[index[i-hh]]*data[index[i-hh]];
81  }
82 
83  Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma));
84  mean = bestmean;
85  sigma = bstd;
86  delete [] index;
87 
88 }
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  // Modified version of ROOT robust EvaluateUni
94  // robust estimator in 1D case MI version
95  // added external factor to include precision of external measurement
96  //
97 
98  if (hh==0)
99  hh=(nvectors+2)/2;
100  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};
101  Int_t *index=new Int_t[nvectors];
102  TMath::Sort(nvectors, data, index, kFALSE);
103  //
104  Int_t nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11);
105  Double_t factor = faclts[0];
106  if (nquant>0){
107  // fix proper normalization - Anja
108  factor = faclts[nquant-1];
109  }
110 
111  //
112  //
113  Double_t sumx =0;
114  Double_t sumx2 =0;
115  Double_t bestmean = 0;
116  Double_t bestsigma = -1;
117  for (Int_t i=0; i<hh; i++){
118  sumx += data[index[i]];
119  sumx2 += data[index[i]]*data[index[i]];
120  }
121  //
122  Double_t kfactor = 2.*externalfactor - externalfactor*externalfactor;
123  Double_t norm = 1./Double_t(hh);
124  for (Int_t i=hh; i<nvectors; i++){
125  Double_t cmean = sumx*norm;
126  Double_t csigma = (sumx2*norm - cmean*cmean*kfactor);
127  if (csigma<bestsigma || bestsigma<0){
128  bestmean = cmean;
129  bestsigma = csigma;
130  }
131  //
132  //
133  sumx += data[index[i]]-data[index[i-hh]];
134  sumx2 += data[index[i]]*data[index[i]]-data[index[i-hh]]*data[index[i-hh]];
135  }
136 
137  Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma));
138  mean = bestmean;
139  sigma = bstd;
140  delete [] index;
141 }
142 
143 
144 //_____________________________________________________________________________
145 Int_t TStatToolkit::Freq(Int_t n, const Int_t *inlist
146  , Int_t *outlist, Bool_t down)
147 {
148  //
149  // Sort eleements according occurancy
150  // The size of output array has is 2*n
151  //
152 
153  Int_t * sindexS = new Int_t[n]; // temp array for sorting
154  Int_t * sindexF = new Int_t[2*n];
155  for (Int_t i=0;i<n;i++) sindexS[i]=0;
156  for (Int_t i=0;i<2*n;i++) sindexF[i]=0;
157  //
158  TMath::Sort(n,inlist, sindexS, down);
159  Int_t last = inlist[sindexS[0]];
160  Int_t val = last;
161  sindexF[0] = 1;
162  sindexF[0+n] = last;
163  Int_t countPos = 0;
164  //
165  // find frequency
166  for(Int_t i=1;i<n; i++){
167  val = inlist[sindexS[i]];
168  if (last == val) sindexF[countPos]++;
169  else{
170  countPos++;
171  sindexF[countPos+n] = val;
172  sindexF[countPos]++;
173  last =val;
174  }
175  }
176  if (last==val) countPos++;
177  // sort according frequency
178  TMath::Sort(countPos, sindexF, sindexS, kTRUE);
179  for (Int_t i=0;i<countPos;i++){
180  outlist[2*i ] = sindexF[sindexS[i]+n];
181  outlist[2*i+1] = sindexF[sindexS[i]];
182  }
183  delete [] sindexS;
184  delete [] sindexF;
185 
186  return countPos;
187 
188 }
189 
190 
191 
192 void TStatToolkit::MedianFilter(TH1 * his1D, Int_t nmedian){
193  //
194  // Algorithm to filter histogram
195  // author: marian.ivanov@cern.ch
196  // Details of algorithm:
197  // http://en.wikipedia.org/w/index.php?title=Median_filter&oldid=582191524
198  // Input parameters:
199  // his1D - input histogam - to be modiefied by Medianfilter
200  // nmendian - number of bins in median filter
201  //
202  Int_t nbins = his1D->GetNbinsX();
203  TVectorD vectorH(nbins);
204  for (Int_t ibin=0; ibin<nbins; ibin++) vectorH[ibin]=his1D->GetBinContent(ibin+1);
205  for (Int_t ibin=0; ibin<nbins; ibin++) {
206  Int_t index0=ibin-nmedian;
207  Int_t index1=ibin+nmedian;
208  if (index0<0) {index1+=-index0; index0=0;}
209  if (index1>=nbins) {index0-=index1-nbins+1; index1=nbins-1;}
210  Double_t value= TMath::Median(index1-index0,&(vectorH.GetMatrixArray()[index0]));
211  his1D->SetBinContent(ibin+1, value);
212  }
213 }
214 
215 
216 
217 Float_t TStatToolkit::GetCOG(const Short_t *arr, Int_t nBins, Float_t xMin, Float_t xMax, Float_t *rms, Float_t *sum)
218 {
219  //
220  // calculate center of gravity rms and sum for array 'arr' with nBins an a x range xMin to xMax
221  // return COG; in case of failure return xMin
222  //
223  Float_t meanCOG = 0;
224  Float_t rms2COG = 0;
225  Float_t sumCOG = 0;
226  Int_t npoints = 0;
227 
228  Float_t binWidth = (xMax-xMin)/(Float_t)nBins;
229 
230  for (Int_t ibin=0; ibin<nBins; ibin++){
231  Float_t entriesI = (Float_t)arr[ibin];
232  Double_t xcenter = xMin+(ibin+0.5)*binWidth;
233  if ( entriesI>0 ){
234  meanCOG += xcenter*entriesI;
235  rms2COG += xcenter*entriesI*xcenter;
236  sumCOG += entriesI;
237  npoints++;
238  }
239  }
240  if ( sumCOG == 0 ) return xMin;
241  meanCOG/=sumCOG;
242 
243  if ( rms ){
244  rms2COG /=sumCOG;
245  (*rms) = TMath::Sqrt(TMath::Abs(meanCOG*meanCOG-rms2COG));
246  if ( npoints == 1 ) (*rms) = binWidth/TMath::Sqrt(12);
247  }
248 
249  if ( sum )
250  (*sum) = sumCOG;
251 
252  return meanCOG;
253 }
254 
255 
256 
260 
261 
262 
263 
264 
265 void TStatToolkit::TestGausFit(Int_t nhistos){
266  //
267  // Test performance of the parabolic - gaussian fit - compare it with
268  // ROOT gauss fit
269  // nhistos - number of histograms to be used for test
270  //
271  TTreeSRedirector *pcstream = new TTreeSRedirector("fitdebug.root","recreate");
272 
273  Float_t *xTrue = new Float_t[nhistos];
274  Float_t *sTrue = new Float_t[nhistos];
275  TVectorD **par1 = new TVectorD*[nhistos];
276  TVectorD **par2 = new TVectorD*[nhistos];
277  TMatrixD dummy(3,3);
278 
279 
280  TH1F **h1f = new TH1F*[nhistos];
281  TF1 *myg = new TF1("myg","gaus");
282  TF1 *fit = new TF1("fit","gaus");
283  gRandom->SetSeed(0);
284 
285  //init
286  for (Int_t i=0;i<nhistos; i++){
287  par1[i] = new TVectorD(3);
288  par2[i] = new TVectorD(3);
289  h1f[i] = new TH1F(Form("h1f%d",i),Form("h1f%d",i),20,-10,10);
290  xTrue[i]= gRandom->Rndm();
291  gSystem->Sleep(2);
292  sTrue[i]= .75+gRandom->Rndm()*.5;
293  myg->SetParameters(1,xTrue[i],sTrue[i]);
294  h1f[i]->FillRandom("myg");
295  }
296 
297  TStopwatch s;
298  s.Start();
299  //standard gaus fit
300  for (Int_t i=0; i<nhistos; i++){
301  h1f[i]->Fit(fit,"0q");
302  (*par1[i])(0) = fit->GetParameter(0);
303  (*par1[i])(1) = fit->GetParameter(1);
304  (*par1[i])(2) = fit->GetParameter(2);
305  }
306  s.Stop();
307  printf("Gaussian fit\t");
308  s.Print();
309 
310  s.Start();
311  //TStatToolkit gaus fit
312  for (Int_t i=0; i<nhistos; i++){
313  TStatToolkit::FitGaus(h1f[i]->GetArray()+1,h1f[i]->GetNbinsX(),h1f[i]->GetXaxis()->GetXmin(),h1f[i]->GetXaxis()->GetXmax(),par2[i],&dummy);
314  }
315 
316  s.Stop();
317  printf("Parabolic fit\t");
318  s.Print();
319  //write stream
320  for (Int_t i=0;i<nhistos; i++){
321  Float_t xt = xTrue[i];
322  Float_t st = sTrue[i];
323  (*pcstream)<<"data"
324  <<"xTrue="<<xt
325  <<"sTrue="<<st
326  <<"pg.="<<(par1[i])
327  <<"pa.="<<(par2[i])
328  <<"\n";
329  }
330  //delete pointers
331  for (Int_t i=0;i<nhistos; i++){
332  delete par1[i];
333  delete par2[i];
334  delete h1f[i];
335  }
336  delete pcstream;
337  delete []h1f;
338  delete []xTrue;
339  delete []sTrue;
340  //
341  delete []par1;
342  delete []par2;
343 
344 }
345 
346 
347 
348 TGraph2D * TStatToolkit::MakeStat2D(TH3 * his, Int_t delta0, Int_t delta1, Int_t type){
349  //
350  //
351  //
352  // delta - number of bins to integrate
353  // type - 0 - mean value
354 
355  TAxis * xaxis = his->GetXaxis();
356  TAxis * yaxis = his->GetYaxis();
357  // TAxis * zaxis = his->GetZaxis();
358  Int_t nbinx = xaxis->GetNbins();
359  Int_t nbiny = yaxis->GetNbins();
360  char name[1000];
361  Int_t icount=0;
362  TGraph2D *graph = new TGraph2D(nbinx*nbiny);
363  TF1 f1("f1","gaus");
364  for (Int_t ix=0; ix<nbinx;ix++)
365  for (Int_t iy=0; iy<nbiny;iy++){
366  Float_t xcenter = xaxis->GetBinCenter(ix);
367  Float_t ycenter = yaxis->GetBinCenter(iy);
368  snprintf(name,1000,"%s_%d_%d",his->GetName(), ix,iy);
369  TH1 *projection = his->ProjectionZ(name,ix-delta0,ix+delta0,iy-delta1,iy+delta1);
370  Float_t stat= 0;
371  if (type==0) stat = projection->GetMean();
372  if (type==1) stat = projection->GetRMS();
373  if (type==2 || type==3){
374  TVectorD vec(10);
375  TStatToolkit::LTM((TH1F*)projection,&vec,0.7);
376  if (type==2) stat= vec[1];
377  if (type==3) stat= vec[0];
378  }
379  if (type==4|| type==5){
380  projection->Fit(&f1);
381  if (type==4) stat= f1.GetParameter(1);
382  if (type==5) stat= f1.GetParameter(2);
383  }
384  //printf("%d\t%f\t%f\t%f\n", icount,xcenter, ycenter, stat);
385  graph->SetPoint(icount,xcenter, ycenter, stat);
386  icount++;
387  }
388  return graph;
389 }
390 
391 TGraphErrors * TStatToolkit::MakeStat1D(TH2 * his, Int_t deltaBin, Double_t fraction, Int_t returnType, Int_t markerStyle, Int_t markerColor){
392  //
393  // function to retrieve the "mean and RMS estimate" of 2D histograms
394  //
395  // Robust statistic to estimate properties of the distribution
396  // See http://en.wikipedia.org/wiki/Trimmed_estimator
397  //
398  // deltaBin - number of bins to integrate (bin+-deltaBin)
399  // fraction - fraction of values for the LTM and for the gauss fit
400  // returnType -
401  // 0 - mean value
402  // 1 - RMS
403  // 2 - LTM mean
404  // 3 - LTM sigma
405  // 4 - Gaus fit mean - on LTM range
406  // 5 - Gaus fit sigma - on LTM range
407  // 6 - Robust bin median
408  // 7 - Gaus fit around maxium of the distribution, fraction is +- relative range of full y scale
409  //
410  TAxis * xaxis = his->GetXaxis();
411  Int_t nbinx = xaxis->GetNbins();
412  char name[1000];
413  Int_t icount=0;
414  //
415  TVectorD vecX(nbinx);
416  TVectorD vecXErr(nbinx);
417  TVectorD vecY(nbinx);
418  TVectorD vecYErr(nbinx);
419  //
420  TF1 f1("f1","gaus");
421  TVectorD vecLTM(10);
422 
423  for (Int_t jx=1; jx<=nbinx;jx++){
424  Int_t ix=jx-1;
425  Float_t xcenter = xaxis->GetBinCenter(jx);
426  snprintf(name,1000,"%s_%d",his->GetName(), ix);
427  TH1 *projection = his->ProjectionY(name,TMath::Max(jx-deltaBin,1),TMath::Min(jx+deltaBin,nbinx));
428  Double_t stat= 0;
429  Double_t err =0;
430  if (projection->Integral()==0) {
431  vecX[icount] = xcenter;
432  vecY[icount] = stat;
433  vecYErr[icount] = err;
434  icount++;
435  delete projection;
436  continue;
437  }
438 
439  TStatToolkit::LTMHisto((TH1F*)projection,vecLTM,fraction);
440  //
441  if (returnType==0) {
442  stat = projection->GetMean();
443  err = projection->GetMeanError();
444  }
445  else if (returnType==1) {
446  stat = projection->GetRMS();
447  err = projection->GetRMSError();
448  }
449  else if (returnType==2 || returnType==3){
450  if (returnType==2) {stat= vecLTM[1]; err =projection->GetRMSError();}
451  if (returnType==3) {stat= vecLTM[2]; err =projection->GetRMSError();}
452  }
453  else if (returnType==4|| returnType==5){
454  f1.SetParameters(vecLTM[0], vecLTM[1], vecLTM[2]+0.05);
455  projection->Fit(&f1,"QN","QN", vecLTM[7]-vecLTM[2], vecLTM[8]+vecLTM[2]);
456  if (returnType==4) {
457  stat= f1.GetParameter(1);
458  err=f1.GetParError(1);
459  }
460  if (returnType==5) {
461  stat= f1.GetParameter(2);
462  err=f1.GetParError(2);
463  }
464  }
465  else if (returnType==6) {
466  stat=RobustBinMedian(projection,fraction);
467  }
468  else if (returnType==7) {
469  const Int_t maxBin = projection->GetMaximumBin();
470  const Double_t max = projection->GetXaxis()->GetBinCenter(maxBin);
471  const Double_t range = fraction*(projection->GetXaxis()->GetXmax()-projection->GetXaxis()->GetXmin());
472  f1.SetParameters(projection->GetMaximum(),
473  projection->GetMean(),
474  projection->GetRMS());
475  f1.SetRange(max-range, max+range);
476  projection->Fit(&f1,"QN","QN", max-range, max+range);
477  stat= f1.GetParameter(1);
478  err=f1.GetParError(1);
479  }
480 
481  vecX[icount]=xcenter;
482  vecY[icount]=stat;
483  vecYErr[icount]=err;
484  icount++;
485  delete projection;
486  }
487  TGraphErrors *graph = new TGraphErrors(icount,vecX.GetMatrixArray(), vecY.GetMatrixArray(),0, vecYErr.GetMatrixArray());
488  graph->SetMarkerStyle(markerStyle);
489  graph->SetMarkerColor(markerColor);
490  return graph;
491 }
492 
493 Double_t TStatToolkit::RobustBinMedian(TH1* hist, Double_t fractionCut/*=1.*/)
494 {
495  // Robust median with average from neighbouring bins
496  const Int_t nbins1D=hist->GetNbinsX();
497  Double_t binMedian=0;
498  Double_t limits[2]={hist->GetBinCenter(1), hist->GetBinCenter(nbins1D)};
499  (void) limits[0];
500 
501  Double_t* integral=hist->GetIntegral();
502  for (Int_t i=1; i<nbins1D-1; i++){
503  if (integral[i-1]<0.5 && integral[i]>=0.5){
504  if (hist->GetBinContent(i-1)+hist->GetBinContent(i)>0){
505  binMedian=hist->GetBinCenter(i);
506  Double_t dIdx=-(integral[i-1]-integral[i]);
507  Double_t dx=(0.5+(0.5-integral[i])/dIdx)*hist->GetBinWidth(i);
508  binMedian+=dx;
509  }
510  }
511  if (integral[i-1]<fractionCut && integral[i]>=fractionCut){
512  limits[0]=hist->GetBinCenter(i-1)-hist->GetBinWidth(i);
513  }
514  if (integral[i]<1-fractionCut && integral[i+1]>=1-fractionCut){
515  limits[1]=hist->GetBinCenter(i+1)+hist->GetBinWidth(i);
516  }
517  }
518 
519  return binMedian;
520 }
521 
522 
523 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){
524  //
525  // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts
526  // returns chi2, fitParam and covMatrix
527  // returns TString with fitted formula
528  //
529 
530  TString formulaStr(formula);
531  TString drawStr(drawCommand);
532  TString cutStr(cuts);
533  TString ferr("1");
534 
535  TString strVal(drawCommand);
536  if (strVal.Contains(":")){
537  TObjArray* valTokens = strVal.Tokenize(":");
538  drawStr = valTokens->At(0)->GetName();
539  ferr = valTokens->At(1)->GetName();
540  delete valTokens;
541  }
542 
543 
544  formulaStr.ReplaceAll("++", "~");
545  TObjArray* formulaTokens = formulaStr.Tokenize("~");
546  Int_t dim = formulaTokens->GetEntriesFast();
547 
548  fitParam.ResizeTo(dim);
549  covMatrix.ResizeTo(dim,dim);
550 
551  TLinearFitter* fitter = new TLinearFitter(dim+1, Form("hyp%d",dim));
552  fitter->StoreData(kTRUE);
553  fitter->ClearPoints();
554 
555  Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start, start);
556  if (entries == -1) {
557  delete formulaTokens;
558  return new TString(TString::Format("ERROR expr: %s\t%s\tEntries==0",drawStr.Data(),cutStr.Data()));
559  }
560  Double_t **values = new Double_t*[dim+1] ;
561  for (Int_t i=0; i<dim+1; i++) values[i]=NULL;
562  //
563  entries = tree->Draw(ferr.Data(), cutStr.Data(), "goff", stop-start, start);
564  if (entries == -1) {
565  delete formulaTokens;
566  delete []values;
567  return new TString(TString::Format("ERROR error part: %s\t%s\tEntries==0",ferr.Data(),cutStr.Data()));
568  }
569  Double_t *errors = new Double_t[entries];
570  memcpy(errors, tree->GetV1(), entries*sizeof(Double_t));
571 
572  for (Int_t i = 0; i < dim + 1; i++){
573  Int_t centries = 0;
574  if (i < dim) centries = tree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start);
575  else centries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start);
576 
577  if (entries != centries) {
578  delete []errors;
579  delete []values;
580  return new TString(TString::Format("ERROR: %s\t%s\tEntries==%d\tEntries2=%d\n",drawStr.Data(),cutStr.Data(),entries,centries));
581  }
582  values[i] = new Double_t[entries];
583  memcpy(values[i], tree->GetV1(), entries*sizeof(Double_t));
584  }
585 
586  // add points to the fitter
587  for (Int_t i = 0; i < entries; i++){
588  Double_t x[1000];
589  for (Int_t j=0; j<dim;j++) x[j]=values[j][i];
590  fitter->AddPoint(x, values[dim][i], errors[i]);
591  }
592 
593  fitter->Eval();
594  if (frac>0.5 && frac<1){
595  fitter->EvalRobust(frac);
596  }else{
597  if (fix0) {
598  fitter->FixParameter(0,0);
599  fitter->Eval();
600  }
601  }
602  fitter->GetParameters(fitParam);
603  fitter->GetCovarianceMatrix(covMatrix);
604  chi2 = fitter->GetChisquare();
605  npoints = entries;
606  TString *preturnFormula = new TString(Form("( %f+",fitParam[0])), &returnFormula = *preturnFormula;
607 
608  for (Int_t iparam = 0; iparam < dim; iparam++) {
609  returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam+1]));
610  if (iparam < dim-1) returnFormula.Append("+");
611  }
612  returnFormula.Append(" )");
613 
614 
615  for (Int_t j=0; j<dim+1;j++) delete [] values[j];
616 
617 
618  delete formulaTokens;
619  delete fitter;
620  delete[] values;
621  delete[] errors;
622  return preturnFormula;
623 }
624 
625 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){
626  //
627  // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts
628  // returns chi2, fitParam and covMatrix
629  // returns TString with fitted formula
630  //
631 
632  TString formulaStr(formula);
633  TString drawStr(drawCommand);
634  TString cutStr(cuts);
635  TString ferr("1");
636 
637  TString strVal(drawCommand);
638  if (strVal.Contains(":")){
639  TObjArray* valTokens = strVal.Tokenize(":");
640  drawStr = valTokens->At(0)->GetName();
641  ferr = valTokens->At(1)->GetName();
642  delete valTokens;
643  }
644 
645 
646  formulaStr.ReplaceAll("++", "~");
647  TObjArray* formulaTokens = formulaStr.Tokenize("~");
648  Int_t dim = formulaTokens->GetEntriesFast();
649 
650  fitParam.ResizeTo(dim);
651  covMatrix.ResizeTo(dim,dim);
652 
653  TLinearFitter* fitter = new TLinearFitter(dim+1, Form("hyp%d",dim));
654  fitter->StoreData(kTRUE);
655  fitter->ClearPoints();
656 
657  Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start, start);
658  if (entries == -1) {
659  delete formulaTokens;
660  return new TString("An ERROR has occured during fitting!");
661  }
662  Double_t **values = new Double_t*[dim+1] ;
663  for (Int_t i=0; i<dim+1; i++) values[i]=NULL;
664  //
665  entries = tree->Draw(ferr.Data(), cutStr.Data(), "goff", stop-start, start);
666  if (entries == -1) {
667  delete formulaTokens;
668  delete [] values;
669  return new TString("An ERROR has occured during fitting!");
670  }
671  Double_t *errors = new Double_t[entries];
672  memcpy(errors, tree->GetV1(), entries*sizeof(Double_t));
673 
674  for (Int_t i = 0; i < dim + 1; i++){
675  Int_t centries = 0;
676  if (i < dim) centries = tree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start);
677  else centries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start);
678 
679  if (entries != centries) {
680  delete []errors;
681  delete []values;
682  delete formulaTokens;
683  return new TString("An ERROR has occured during fitting!");
684  }
685  values[i] = new Double_t[entries];
686  memcpy(values[i], tree->GetV1(), entries*sizeof(Double_t));
687  }
688 
689  // add points to the fitter
690  for (Int_t i = 0; i < entries; i++){
691  Double_t x[1000];
692  for (Int_t j=0; j<dim;j++) x[j]=values[j][i];
693  fitter->AddPoint(x, values[dim][i], errors[i]);
694  }
695  if (constrain>0){
696  for (Int_t i = 0; i < dim; i++){
697  Double_t x[1000];
698  for (Int_t j=0; j<dim;j++) if (i!=j) x[j]=0;
699  x[i]=1.;
700  fitter->AddPoint(x, 0, constrain);
701  }
702  }
703 
704 
705  fitter->Eval();
706  if (frac>0.5 && frac<1){
707  fitter->EvalRobust(frac);
708  }
709  fitter->GetParameters(fitParam);
710  fitter->GetCovarianceMatrix(covMatrix);
711  chi2 = fitter->GetChisquare();
712  npoints = entries;
713 
714  TString *preturnFormula = new TString(Form("( %f+",fitParam[0])), &returnFormula = *preturnFormula;
715 
716  for (Int_t iparam = 0; iparam < dim; iparam++) {
717  returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam+1]));
718  if (iparam < dim-1) returnFormula.Append("+");
719  }
720  returnFormula.Append(" )");
721 
722  for (Int_t j=0; j<dim+1;j++) delete [] values[j];
723 
724 
725 
726  delete formulaTokens;
727  delete fitter;
728  delete[] values;
729  delete[] errors;
730  return preturnFormula;
731 }
732 
733 
734 
735 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){
736  //
737  // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts
738  // returns chi2, fitParam and covMatrix
739  // returns TString with fitted formula
740  //
741 
742  TString formulaStr(formula);
743  TString drawStr(drawCommand);
744  TString cutStr(cuts);
745  TString ferr("1");
746 
747  TString strVal(drawCommand);
748  if (strVal.Contains(":")){
749  TObjArray* valTokens = strVal.Tokenize(":");
750  drawStr = valTokens->At(0)->GetName();
751  ferr = valTokens->At(1)->GetName();
752  delete valTokens;
753  }
754 
755 
756  formulaStr.ReplaceAll("++", "~");
757  TObjArray* formulaTokens = formulaStr.Tokenize("~");
758  Int_t dim = formulaTokens->GetEntriesFast();
759 
760  fitParam.ResizeTo(dim);
761  covMatrix.ResizeTo(dim,dim);
762  TString fitString="x0";
763  for (Int_t i=1; i<dim; i++) fitString+=Form("++x%d",i);
764  TLinearFitter* fitter = new TLinearFitter(dim, fitString.Data());
765  fitter->StoreData(kTRUE);
766  fitter->ClearPoints();
767 
768  Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start, start);
769  if (entries == -1) {
770  delete formulaTokens;
771  return new TString("An ERROR has occured during fitting!");
772  }
773  Double_t **values = new Double_t*[dim+1] ;
774  for (Int_t i=0; i<dim+1; i++) values[i]=NULL;
775  //
776  entries = tree->Draw(ferr.Data(), cutStr.Data(), "goff", stop-start, start);
777  if (entries == -1) {
778  delete []values;
779  delete formulaTokens;
780  return new TString("An ERROR has occured during fitting!");
781  }
782  Double_t *errors = new Double_t[entries];
783  memcpy(errors, tree->GetV1(), entries*sizeof(Double_t));
784 
785  for (Int_t i = 0; i < dim + 1; i++){
786  Int_t centries = 0;
787  if (i < dim) centries = tree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start);
788  else centries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start);
789 
790  if (entries != centries) {
791  delete []errors;
792  delete []values;
793  delete formulaTokens;
794  return new TString("An ERROR has occured during fitting!");
795  }
796  values[i] = new Double_t[entries];
797  memcpy(values[i], tree->GetV1(), entries*sizeof(Double_t));
798  }
799 
800  // add points to the fitter
801  for (Int_t i = 0; i < entries; i++){
802  Double_t x[1000];
803  for (Int_t j=0; j<dim;j++) x[j]=values[j][i];
804  fitter->AddPoint(x, values[dim][i], errors[i]);
805  }
806 
807  fitter->Eval();
808  if (frac>0.5 && frac<1){
809  fitter->EvalRobust(frac);
810  }
811  fitter->GetParameters(fitParam);
812  fitter->GetCovarianceMatrix(covMatrix);
813  chi2 = fitter->GetChisquare();
814  npoints = entries;
815 
816  TString *preturnFormula = new TString("("), &returnFormula = *preturnFormula;
817 
818  for (Int_t iparam = 0; iparam < dim; iparam++) {
819  returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam]));
820  if (iparam < dim-1) returnFormula.Append("+");
821  }
822  returnFormula.Append(" )");
823 
824 
825  for (Int_t j=0; j<dim+1;j++) delete [] values[j];
826 
827  delete formulaTokens;
828  delete fitter;
829  delete[] values;
830  delete[] errors;
831  return preturnFormula;
832 }
833 
834 
835 
836 
837 
838 Int_t TStatToolkit::GetFitIndex(const TString fString, const TString subString){
839  //
840  // fitString - ++ separated list of fits
841  // substring - ++ separated list of the requiered substrings
842  //
843  // return the last occurance of substring in fit string
844  //
845  TObjArray *arrFit = fString.Tokenize("++");
846  TObjArray *arrSub = subString.Tokenize("++");
847  Int_t index=-1;
848  for (Int_t i=0; i<arrFit->GetEntries(); i++){
849  Bool_t isOK=kTRUE;
850  TString str =arrFit->At(i)->GetName();
851  for (Int_t isub=0; isub<arrSub->GetEntries(); isub++){
852  if (str.Contains(arrSub->At(isub)->GetName())==0) isOK=kFALSE;
853  }
854  if (isOK) index=i;
855  }
856  delete arrFit;
857  delete arrSub;
858  return index;
859 }
860 
861 
862 TString TStatToolkit::FilterFit(const TString &input, const TString filter, TVectorD &param, TMatrixD & covar){
863  //
864  // Filter fit expression make sub-fit
865  //
866  TObjArray *array0= input.Tokenize("++");
867  TObjArray *array1= filter.Tokenize("++");
868  //TString *presult=new TString("(0");
869  TString result="(0.0";
870  for (Int_t i=0; i<array0->GetEntries(); i++){
871  Bool_t isOK=kTRUE;
872  TString str(array0->At(i)->GetName());
873  for (Int_t j=0; j<array1->GetEntries(); j++){
874  if (str.Contains(array1->At(j)->GetName())==0) isOK=kFALSE;
875  }
876  if (isOK) {
877  result+="+"+str;
878  result+=Form("*(%f)",param[i+1]);
879  printf("%f\t%f\t%s\n",param[i+1], TMath::Sqrt(covar(i+1,i+1)),str.Data());
880  }
881  }
882  result+="-0.)";
883  delete array0;
884  delete array1;
885  return result;
886 }
887 
888 void TStatToolkit::Update1D(Double_t delta, Double_t sigma, Int_t s1, TMatrixD &vecXk, TMatrixD &covXk){
889  //
890  // Update parameters and covariance - with one measurement
891  // Input:
892  // vecXk - input vector - Updated in function
893  // covXk - covariance matrix - Updated in function
894  // delta, sigma, s1 - new measurement, rms of new measurement and the index of measurement
895  const Int_t knMeas=1;
896  Int_t knElem=vecXk.GetNrows();
897 
898  TMatrixD mat1(knElem,knElem); // update covariance matrix
899  TMatrixD matHk(1,knElem); // vector to mesurement
900  TMatrixD vecYk(knMeas,1); // Innovation or measurement residual
901  TMatrixD matHkT(knElem,knMeas); // helper matrix Hk transpose
902  TMatrixD matSk(knMeas,knMeas); // Innovation (or residual) covariance
903  TMatrixD matKk(knElem,knMeas); // Optimal Kalman gain
904  TMatrixD covXk2(knElem,knElem); // helper matrix
905  TMatrixD covXk3(knElem,knElem); // helper matrix
906  TMatrixD vecZk(1,1);
907  TMatrixD measR(1,1);
908  vecZk(0,0)=delta;
909  measR(0,0)=sigma*sigma;
910  //
911  // reset matHk
912  for (Int_t iel=0;iel<knElem;iel++)
913  for (Int_t ip=0;ip<knMeas;ip++) matHk(ip,iel)=0;
914  //mat1
915  for (Int_t iel=0;iel<knElem;iel++) {
916  for (Int_t jel=0;jel<knElem;jel++) mat1(iel,jel)=0;
917  mat1(iel,iel)=1;
918  }
919  //
920  matHk(0, s1)=1;
921  vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual
922  matHkT=matHk.T(); matHk.T();
923  matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance
924  matSk.Invert();
925  matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
926  vecXk += matKk*vecYk; // updated vector
927  covXk2= (mat1-(matKk*matHk));
928  covXk3 = covXk2*covXk;
929  covXk = covXk3;
930  Int_t nrows=covXk3.GetNrows();
931 
932  for (Int_t irow=0; irow<nrows; irow++)
933  for (Int_t icol=0; icol<nrows; icol++){
934  // rounding problems - make matrix again symteric
935  covXk(irow,icol)=(covXk3(irow,icol)+covXk3(icol,irow))*0.5;
936  }
937 }
938 
939 
940 
941 void TStatToolkit::Constrain1D(const TString &input, const TString filter, TVectorD &param, TMatrixD & covar, Double_t mean, Double_t sigma){
942  //
943  // constrain linear fit
944  // input - string description of fit function
945  // filter - string filter to select sub fits
946  // param,covar - parameters and covariance matrix of the fit
947  // mean,sigma - new measurement uning which the fit is updated
948  //
949 
950  TObjArray *array0= input.Tokenize("++");
951  TObjArray *array1= filter.Tokenize("++");
952  TMatrixD paramM(param.GetNrows(),1);
953  for (Int_t i=0; i<=array0->GetEntries(); i++){paramM(i,0)=param(i);}
954 
955  if (filter.Length()==0){
956  TStatToolkit::Update1D(mean, sigma, 0, paramM, covar);//
957  }else{
958  for (Int_t i=0; i<array0->GetEntries(); i++){
959  Bool_t isOK=kTRUE;
960  TString str(array0->At(i)->GetName());
961  for (Int_t j=0; j<array1->GetEntries(); j++){
962  if (str.Contains(array1->At(j)->GetName())==0) isOK=kFALSE;
963  }
964  if (isOK) {
965  TStatToolkit::Update1D(mean, sigma, i+1, paramM, covar);//
966  }
967  }
968  }
969  for (Int_t i=0; i<=array0->GetEntries(); i++){
970  param(i)=paramM(i,0);
971  }
972  delete array0;
973  delete array1;
974 }
975 
976 TString TStatToolkit::MakeFitString(const TString &input, const TVectorD &param, const TMatrixD & covar, Bool_t verbose){
977  //
978  //
979  //
980  TObjArray *array0= input.Tokenize("++");
981  TString result=Form("(%f",param[0]);
982  printf("%f\t%f\t\n", param[0], TMath::Sqrt(covar(0,0)));
983  for (Int_t i=0; i<array0->GetEntries(); i++){
984  TString str(array0->At(i)->GetName());
985  result+="+"+str;
986  result+=Form("*(%f)",param[i+1]);
987  if (verbose) printf("%f\t%f\t%s\n", param[i+1], TMath::Sqrt(covar(i+1,i+1)),str.Data());
988  }
989  result+="-0.)";
990  delete array0;
991  return result;
992 }
993 
994 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){
995  //
996  // Query a graph errors
997  // return TGraphErrors specified by expr and cut
998  // Example usage TStatToolkit::MakeGraphError(tree,"Y:X:ErrY","X>0", 25,2,0.4)
999  // tree - tree with variable
1000  // expr - examp
1001  const Int_t entries = tree->Draw(expr,cut,"goff",drawEntries,firstEntry);
1002 
1003  if (entries<=0) {
1004  ::Error("TStatToolkit::MakeGraphError","Empty or Not valid expression (%s) or cut *%s)", expr,cut);
1005  return 0;
1006  }
1007  if ( tree->GetV2()==0){
1008  ::Error("TStatToolkit::MakeGraphError","Not valid expression (%s) ", expr);
1009  return 0;
1010  }
1011  TGraphErrors * graph=0;
1012  if ( tree->GetV3()!=0){
1013  graph = new TGraphErrors (entries, tree->GetV2(),tree->GetV1(),0,tree->GetV3());
1014  }else{
1015  graph = new TGraphErrors (entries, tree->GetV2(),tree->GetV1(),0,0);
1016  }
1017 
1018  graph->SetMarkerStyle(mstyle);
1019  graph->SetMarkerColor(mcolor);
1020  graph->SetLineColor(mcolor);
1021  graph->SetTitle(expr);
1022  TString chstring(expr);
1023  TObjArray *charray = chstring.Tokenize(":");
1024  graph->GetXaxis()->SetTitle(charray->At(1)->GetName());
1025  graph->GetYaxis()->SetTitle(charray->At(0)->GetName());
1026  THashList * metaData = (THashList*) tree->GetUserInfo()->FindObject("metaTable");
1027  if ( metaData != NULL){
1028  TNamed *nmdTitle0 = TStatToolkit::GetMetadata(tree,Form("%s.Title",charray->At(0)->GetName()));
1029  TNamed *nmdTitle1 = TStatToolkit::GetMetadata(tree,Form("%s.Title",charray->At(1)->GetName()));
1030  TNamed *nmdYAxis = TStatToolkit::GetMetadata(tree,Form("%s.AxisTitle",charray->At(0)->GetName()));
1031  TNamed *nmdXAxis = TStatToolkit::GetMetadata(tree,Form("%s.AxisTitle",charray->At(1)->GetName()));
1032  //
1033  TString grTitle=charray->At(0)->GetName();
1034  if (nmdTitle0) grTitle=nmdTitle0->GetTitle();
1035  if (nmdTitle1) {
1036  grTitle+=":";
1037  grTitle+=nmdTitle1->GetTitle();
1038  }else{
1039  grTitle+=":";
1040  grTitle+=charray->At(1)->GetName();
1041  }
1042  if (nmdYAxis) {graph->GetYaxis()->SetTitle(nmdYAxis->GetTitle());}
1043  if (nmdXAxis) {graph->GetXaxis()->SetTitle(nmdXAxis->GetTitle());}
1044  graph->SetTitle(grTitle.Data());
1045  }
1046  delete charray;
1047  if (msize>0) graph->SetMarkerSize(msize);
1048  for(Int_t i=0;i<graph->GetN();i++) graph->GetX()[i]+=offset;
1049  //
1050  if (tree->GetVar(1)->IsString()){ // use string for axis description
1051  TAxis * axis = tree->GetHistogram()->GetXaxis();
1052  axis->Copy(*(graph->GetXaxis()));
1053  }
1054  if (tree->GetVar(0)->IsString()){
1055  TAxis * axis = tree->GetHistogram()->GetYaxis();
1056  axis->Copy(*(graph->GetYaxis()));
1057  }
1058  graph->Sort();
1059  return graph;
1060 
1061 }
1062 
1078 THashList* TStatToolkit::AddMetadata(TTree* tree, const char *varTagName,const char *varTagValue){
1079  if (!tree) return NULL;
1080  THashList * metaData = (THashList*) tree->GetUserInfo()->FindObject("metaTable");
1081  if (metaData == NULL){
1082  metaData=new THashList;
1083  metaData->SetName("metaTable");
1084  tree->GetUserInfo()->AddLast(metaData);
1085  }
1086  if (varTagName!=NULL && varTagValue!=NULL){
1087  TNamed * named = TStatToolkit::GetMetadata(tree, varTagName,NULL,kTRUE);
1088  if (named==NULL){
1089  metaData->AddLast(new TNamed(varTagName,varTagValue));
1090  }else{
1091  named->SetTitle(varTagValue);
1092  }
1093  }
1094  return metaData;
1095 }
1096 
1106 TNamed* TStatToolkit::GetMetadata(TTree* tree, const char *varTagName, TString * prefix,Bool_t fullMatch){
1107  //
1108 
1109  if (!tree) return 0;
1110  TTree * treeMeta=tree;
1111 
1112  THashList * metaData = (THashList*) treeMeta->GetUserInfo()->FindObject("metaTable");
1113  if (metaData == NULL){
1114  metaData=new THashList;
1115  metaData->SetName("metaTable");
1116  tree->GetUserInfo()->AddLast(metaData);
1117  return 0;
1118  }
1119  TNamed * named = (TNamed*)metaData->FindObject(varTagName);
1120  if (named || fullMatch) return named;
1121 
1122  TString metaName(varTagName);
1123  Int_t nDots= metaName.CountChar('.');
1124  if (prefix!=NULL) *prefix="";
1125  if (nDots>1){ //check if friend name expansion needed
1126  while (nDots>1){
1127  TList *fList= treeMeta->GetListOfFriends();
1128  if (fList!=NULL){
1129  Int_t nFriends= fList->GetEntries();
1130  for (Int_t kf=0; kf<nFriends; kf++){
1131  TPRegexp regFriend(TString::Format("^%s.",fList->At(kf)->GetName()).Data());
1132  if (metaName.Contains(regFriend)){
1133  treeMeta=treeMeta->GetFriend(fList->At(kf)->GetName());
1134  regFriend.Substitute(metaName,"");
1135  if (prefix!=NULL){
1136  (*prefix)+=fList->At(kf)->GetName();
1137  // (*prefix)+=""; /// FIX use prefix as it is
1138  }
1139  }
1140  }
1141  }
1142  if (nDots == metaName.CountChar('.')) break;
1143  nDots=metaName.CountChar('.');
1144  }
1145  }
1146 
1147  named = (TNamed*)metaData->FindObject(metaName.Data());
1148  return named;
1149 
1150 }
1151 
1152 
1153 TGraph * TStatToolkit::MakeGraphSparse(TTree * tree, const char * expr, const char * cut, Int_t mstyle, Int_t mcolor, Float_t msize, Float_t offset){
1154  //
1155  // Make a sparse draw of the variables
1156  // Format of expr : Var:Run or Var:Run:ErrorY or Var:Run:ErrorY:ErrorX
1157  // offset : points can slightly be shifted in x for better visibility with more graphs
1158  //
1159  // Patrick Reichelt and Marian Ivanov
1160  // maintained and updated by Marian Ivanov
1161  const Int_t entries = tree->Draw(expr,cut,"goff");
1162  if (entries<=0) {
1163  ::Error("TStatToolkit::MakeGraphSparse","Empty or Not valid expression (%s) or cut (%s)", expr, cut);
1164  return 0;
1165  }
1166 
1167  Double_t *graphY, *graphX;
1168  graphY = tree->GetV1();
1169  graphX = tree->GetV2();
1170 
1171  // sort according to run number
1172  Int_t *index = new Int_t[entries*4];
1173  TMath::Sort(entries,graphX,index,kFALSE);
1174 
1175  // define arrays for the new graph
1176  Double_t *unsortedX = new Double_t[entries];
1177  Int_t *runNumber = new Int_t[entries];
1178  Double_t count = 0.5;
1179 
1180  // evaluate arrays for the new graph according to the run-number
1181  Int_t icount=0;
1182  //first entry
1183  unsortedX[index[0]] = count;
1184  runNumber[0] = graphX[index[0]];
1185  // loop the rest of entries
1186  for(Int_t i=1;i<entries;i++)
1187  {
1188  if(graphX[index[i]]==graphX[index[i-1]])
1189  unsortedX[index[i]] = count;
1190  else if(graphX[index[i]]!=graphX[index[i-1]]){
1191  count++;
1192  icount++;
1193  unsortedX[index[i]] = count;
1194  runNumber[icount]=graphX[index[i]];
1195  }
1196  }
1197 
1198  // count the number of xbins (run-wise) for the new graph
1199  const Int_t newNbins = int(count+0.5);
1200  Double_t *newBins = new Double_t[newNbins+1];
1201  for(Int_t i=0; i<=count+1;i++){
1202  newBins[i] = i;
1203  }
1204 
1205  // define and fill the new graph
1206  TGraph *graphNew = 0;
1207  if (tree->GetV3()) {
1208  if (tree->GetV4()) {
1209  graphNew = new TGraphErrors(entries,unsortedX,graphY,tree->GetV4(),tree->GetV3());
1210  }
1211  else { graphNew = new TGraphErrors(entries,unsortedX,graphY,0,tree->GetV3()); }
1212  }
1213  else { graphNew = new TGraphErrors(entries,unsortedX,graphY,0,0); }
1214  // with "Set(...)", the x-axis is being sorted
1215  graphNew->GetXaxis()->Set(newNbins,newBins);
1216 
1217  // set the bins for the x-axis, apply shifting of points
1218  Char_t xName[50];
1219  for(Int_t i=0;i<count;i++){
1220  snprintf(xName,50,"%d",runNumber[i]);
1221  graphNew->GetXaxis()->SetBinLabel(i+1,xName);
1222  graphNew->GetX()[i]+=offset;
1223  }
1224  if (tree->GetVar(1)->IsInteger() && strlen(tree->GetHistogram()->GetXaxis()->GetBinLabel(1))>0){
1225  for(Int_t i=0;i<count;i++){
1226  graphNew->GetXaxis()->SetBinLabel(i+1,tree->GetHistogram()->GetXaxis()->GetBinLabel(i+1));
1227  }
1228  }
1229  if (tree->GetVar(0)->IsInteger() && strlen(tree->GetHistogram()->GetXaxis()->GetBinLabel(1))>0 ){
1230  for(Int_t i=0;i<count;i++){
1231  graphNew->GetYaxis()->SetBinLabel(i+1,tree->GetHistogram()->GetYaxis()->GetBinLabel(i+1));
1232  }
1233  }
1234 
1235  graphNew->GetHistogram()->SetTitle("");
1236  graphNew->SetMarkerStyle(mstyle);
1237  graphNew->SetMarkerColor(mcolor); graphNew->SetLineColor(mcolor);
1238  if (msize>0) { graphNew->SetMarkerSize(msize); graphNew->SetLineWidth(msize); }
1239  delete [] unsortedX;
1240  delete [] runNumber;
1241  delete [] index;
1242  delete [] newBins;
1243  //
1244  TString chstring(expr);
1245  if (cut) chstring+=TString::Format(" ( %s )", cut);
1246  graphNew->SetTitle(chstring);
1247 
1248  THashList * metaData = (THashList*) tree->GetUserInfo()->FindObject("metaTable");
1249  if (metaData != NULL){
1250  chstring=expr;
1251  TObjArray *charray = chstring.Tokenize(":");
1252  graphNew->GetXaxis()->SetTitle(charray->At(1)->GetName());
1253  graphNew->GetYaxis()->SetTitle(charray->At(0)->GetName());
1254  TNamed *nmdTitle0 = TStatToolkit::GetMetadata(tree,Form("%s.Title",charray->At(0)->GetName()));
1255  TNamed *nmdTitle1 = TStatToolkit::GetMetadata(tree,Form("%s.Title",charray->At(1)->GetName()));
1256  TNamed *nmdYAxis = TStatToolkit::GetMetadata(tree,Form("%s.AxisTitle",charray->At(0)->GetName()));
1257  TNamed *nmdXAxis = TStatToolkit::GetMetadata(tree,Form("%s.AxisTitle",charray->At(1)->GetName()));
1258  //
1259  TString grTitle=charray->At(0)->GetName();
1260  if (nmdTitle0) grTitle=nmdTitle0->GetTitle();
1261  if (nmdTitle1) {
1262  grTitle+=":";
1263  grTitle+=nmdTitle1->GetTitle();
1264  }else{
1265  grTitle+=":";
1266  grTitle+=charray->At(1)->GetName();
1267  }
1268  if (cut) grTitle+=TString::Format(" ( %s )", cut);
1269  graphNew->SetTitle(grTitle);
1270  if (nmdYAxis) {graphNew->GetYaxis()->SetTitle(nmdYAxis->GetTitle());}
1271  if (nmdXAxis) {graphNew->GetXaxis()->SetTitle(nmdXAxis->GetTitle());}
1272  delete charray;
1273  }
1274  return graphNew;
1275 }
1276 
1299 TMultiGraph * TStatToolkit::MakeMultGraph(TTree * tree, const char *groupName, const char* expr, const char * cut, const char * styles, const char *colors, Bool_t drawSparse, Float_t msize, Float_t sigmaRange, TLegend * legend, Bool_t comp){
1300 
1301  TMultiGraph *multiGraph=new TMultiGraph(groupName,groupName);
1302  TObjArray * exprVars=TString(expr).Tokenize(":");
1303 
1304  if (exprVars->GetEntries()<2) {
1305  ::Error("MakeMultGraph","NotValid expression %s",expr);
1306  delete exprVars;
1307  return 0;
1308  }
1309  TObjArray*exprVarArray = TString(exprVars->At(0)->GetName()).Tokenize(";");
1310  TObjArray*exprVarErrArray=(exprVars->GetEntries()>2)? TString(exprVars->At(2)->GetName()).Tokenize(";"):0;
1311  TString cutString=(cut!=0)?cut:"1";
1312  if (cutString.Length()==0) cutString="1";
1313  TObjArray*exprCutArray= cutString.Tokenize(";");
1314 
1315  //determine marker style and line style
1316  const char* markerstyles;
1317  const char* linestyles=0;
1318 
1319  if(TString(styles).Contains(",")){
1320  TString styls=TString(styles).ReplaceAll(" ","");
1321  TObjArray * stylarr=styls.Tokenize(",");
1322  markerstyles=TString(stylarr->At(0)->GetName());
1323  linestyles=TString(stylarr->At(1)->GetName());
1324  }
1325  else markerstyles = styles;
1326 
1327  //determine marker style and line style
1328  const char* markercolors;
1329  const char* linecolors=0;
1330 
1331  if(TString(colors).Contains(",")){
1332  TString cols=TString(colors).ReplaceAll(" ","");
1333  TObjArray * colarr=cols.Tokenize(",");
1334  markercolors=TString(colarr->At(0)->GetName());
1335  linecolors=TString(colarr->At(1)->GetName());
1336  }
1337  else markercolors = colors;
1338 // ::Info("MakeMultGraph","Linestyle %d", atoi(linestyles));
1339  //Int markerStyle=AliDrawStyle::GetMarkerStyles(markers);
1340  Int_t notOK=(exprVarArray->GetEntries()<1 && exprCutArray->GetEntries()<2);
1341  if (exprVarErrArray) notOK+=2*(exprVarArray->GetEntriesFast()!=exprVarErrArray->GetEntriesFast());
1342  if (notOK>0){
1343  ::Error("MakeMultGraph","Not compatible arrays of variables:err variables or cuts - Problem %d", notOK);
1344  exprVarArray->Print();
1345  exprCutArray->Print();
1346  if (exprVarErrArray) exprVarErrArray->Print();
1347  delete exprVars;
1348  return 0;
1349  }
1350  Int_t nCuts= TMath::Max(exprCutArray->GetEntries()-1,1);
1351  Int_t nExpr= exprVarArray->GetEntries();
1352  Int_t ngraphs = nCuts*nExpr;
1353  Double_t minValue=1;
1354  Double_t maxValue=-1;
1355  TVectorF vecMean(ngraphs);
1356  Bool_t flag = kFALSE;
1357  for (Int_t iGraph=0; iGraph<ngraphs; iGraph++){
1358  Int_t color=AliDrawStyle::GetMarkerColor(markercolors,iGraph);
1359  Int_t marker=AliDrawStyle::GetMarkerStyle(markerstyles, iGraph);
1360  Int_t iCut =iGraph%nCuts;
1361  Int_t iExpr=iGraph/nCuts;
1362  const char *expName=exprVarArray->At(iExpr)->GetName();
1363  const char *xName=exprVars->At(1)->GetName();
1364  const char *expErrName=(exprVarErrArray==NULL)? NULL:exprVarErrArray->At(iExpr)->GetName();
1365  TString cCut=exprCutArray->At(0)->GetName();
1366  if (nCuts>1){
1367  cCut+="&&";
1368  cCut+=exprCutArray->At(iCut+1)->GetName();
1369  }
1370  TGraph * gr = 0;
1371  if (drawSparse){
1372  if (exprVarErrArray==NULL){
1373  gr=TStatToolkit::MakeGraphSparse(tree, TString::Format("%s:%s",expName,xName).Data(),cCut.Data(), marker,color,msize);
1374  }else{
1375  gr=TStatToolkit::MakeGraphSparse(tree, TString::Format("%s:%s:%s",expName,xName,expErrName).Data(),cCut.Data(), marker,color,msize);
1376  }
1377  }else{
1378  if (exprVarErrArray==NULL){
1379  gr=TStatToolkit::MakeGraphErrors(tree, TString::Format("%s:%s",expName,xName).Data(),cCut.Data(), marker,color,msize);
1380  }else{
1381  gr=TStatToolkit::MakeGraphErrors(tree, TString::Format("%s:%s:%s",expName,xName,expErrName).Data(),cCut.Data(), marker,color,msize);
1382  }
1383  }
1384 
1385  if(linestyles!=0) gr->SetLineStyle(AliDrawStyle::GetLineStyle(linestyles,iGraph));
1386  if(linecolors!=0) gr->SetLineColor(AliDrawStyle::GetLineColor(linecolors,iGraph));
1387 // ::Info("MakeMultGraph","Linestyle2 %d", AliDrawStyle::GetLineStyle(linestyles,iGraph));
1388 
1389  if (gr) {
1390  if (marker<=0) { // explicitly specify draw options - try to avoid bug in TMultiGraph draw in xaxis definition
1391  multiGraph->Add(gr);
1392  }else{
1393  if (iGraph==0){
1394  multiGraph->Add(gr,"ap");
1395  }else{
1396  multiGraph->Add(gr,"p");
1397  }
1398  }
1399  }
1400  Double_t meanT,rmsT=0;
1401  if (gr==NULL){
1402  if(comp){
1403  ::Error("MakeMultGraph","Not valid expression %s or cut %s - return",expName,cCut.Data());
1404  delete exprVarArray;
1405  delete exprVarErrArray;
1406  delete exprCutArray;
1407  return 0;
1408  }
1409  else{
1410  ::Error("MakeMultGraph","Not valid sub-expression %s or cut %s - continue",expName,cCut.Data());
1411  continue;
1412  }
1413  }
1414  flag=kTRUE;
1415  if (gr->GetN()>2){
1416  TStatToolkit::EvaluateUni(gr->GetN(),gr->GetY(), meanT,rmsT, TMath::Max(0.75*gr->GetN(),1.));
1417  }else{
1418  meanT=TMath::Median(gr->GetN(), gr->GetY());
1419  rmsT=TMath::RMS(gr->GetN(), gr->GetY());
1420  }
1421  if (maxValue<minValue){
1422  maxValue=meanT+sigmaRange*rmsT;
1423  minValue=meanT-sigmaRange*rmsT;
1424  }
1425  vecMean[iGraph]=meanT;
1426  if (minValue>meanT-sigmaRange*rmsT) minValue=meanT-sigmaRange*rmsT;
1427  if (maxValue<meanT+sigmaRange*rmsT) maxValue=meanT+sigmaRange*rmsT;
1428  TString prefix="";
1429  gr->SetName(expName);
1430  TNamed*named = TStatToolkit::GetMetadata(tree,TString::Format("%s.Title",expName).Data(),&prefix);
1431  if (named) {
1432  gr->SetTitle(named->GetTitle());
1433  } else{
1434  gr->SetTitle(expName);
1435  }
1436  if (legend){
1437  TString legendName="";
1438  named = TStatToolkit::GetMetadata(tree,TString::Format("%s.Legend",expName).Data(),&prefix);
1439  if (named) {
1440  if (prefix.Length()>0){
1441  TString dummy=prefix+".Legend";
1442  if (TStatToolkit::GetMetadata(tree,dummy.Data())) {
1443  prefix=TStatToolkit::GetMetadata(tree,dummy.Data())->GetTitle();
1444  }
1445  legendName+=prefix.Data();
1446  legendName+=" ";
1447  }
1448  legendName+=named->GetTitle();
1449  }else{
1450  legendName=expName;
1451  }
1452  if (nCuts>1){
1453  legendName+=" ";
1454  named = TStatToolkit::GetMetadata(tree,TString::Format("%s.Legend",exprCutArray->At(iCut+1)->GetName()).Data(),&prefix);
1455  if (named) {
1456  legendName+=named->GetTitle();
1457  }else{
1458  legendName=exprCutArray->At(iCut+1)->GetName();
1459  }
1460  }
1461  legend->AddEntry(gr,legendName,"p");
1462  }
1463 
1464 
1465  }
1466 
1467  if(!flag){
1468  ::Error("Test::","Number of graphs 0 -return 0");
1469  delete exprVarArray;
1470  delete exprVarErrArray;
1471  //delete exprColors;
1472  //delete exprMarkers;
1473  return 0;
1474  }
1475  Double_t rmsGraphs = TMath::RMS(ngraphs, vecMean.GetMatrixArray());
1476  minValue-=sigmaRange*rmsGraphs;
1477  maxValue+=sigmaRange*rmsGraphs;
1478  Double_t xmin=0,xmax=0;
1479  for (Int_t igr=0; igr<multiGraph->GetListOfGraphs()->GetSize(); igr++){
1480  TGraph* gr = (TGraph*)(multiGraph->GetListOfGraphs()->At(igr));
1481  gr->SetMinimum(minValue);
1482  gr->SetMaximum(maxValue);
1483  if (igr==0){
1484  xmin=gr->GetXaxis()->GetXmin();
1485  xmax=gr->GetXaxis()->GetXmax();
1486  }
1487  if (xmin>gr->GetXaxis()->GetXmin()) xmin=gr->GetXaxis()->GetXmin();
1488  if (xmax<gr->GetXaxis()->GetXmax()) xmax=gr->GetXaxis()->GetXmax();
1489  }
1490  // multiGraph->GetXaxis()->SetLimits(xmin,xmax); // BUG/FEATURE mutli graph axis not defined in this moment (pointer==0)
1491  for (Int_t igr=0; igr<multiGraph->GetListOfGraphs()->GetSize(); igr++) {
1492  TGraph *gr = (TGraph *)(multiGraph->GetListOfGraphs()->At(igr));
1493  if (gr!=NULL) gr->GetXaxis()->SetLimits(xmin,xmax);
1494  }
1495  multiGraph->SetMinimum(minValue);
1496  multiGraph->SetMaximum(maxValue);
1497  delete exprVarArray;
1498  delete exprVarErrArray;
1499  delete exprCutArray;
1500  return multiGraph;
1501 }
1502 
1503 
1504 void TStatToolkit::DrawMultiGraph(TMultiGraph *graph, Option_t * option){
1505  //
1506  // TMultiGraph::Draw does not handle properly custom binning of the subcomponents
1507  // Using this function custom labels of the first graphs are used
1508  //
1509  static TPRegexp regAxis("^a");
1510  if (!graph) return;
1511  TList * grArray = graph->GetListOfGraphs();
1512  if (grArray==NULL) return;
1513  TGraph * gr0=(TGraph*)(grArray->At(0));
1514  if (gr0==NULL) return;
1515  gr0->Draw(option); // defaults frame/axis defined by first graph
1516  Int_t ngr=grArray->GetEntries();
1517  TString option2(option);
1518  regAxis.Substitute(option2,"");
1519  for (Int_t igr=1; igr<ngr; igr++){
1520  grArray->At(igr)->Draw(option2.Data());
1521  }
1522 
1523 
1524 }
1525 
1526 
1527 
1528 
1529 //
1530 // functions used for the trending
1531 //
1532 
1533 
1534 Int_t TStatToolkit::MakeStatAlias(TTree * tree, const char * expr, const char * cut, const char * alias)
1535 {
1536  //
1537  // Add alias using statistical values of a given variable.
1538  // (by MI, Patrick Reichelt)
1539  //
1540  // tree - input tree
1541  // expr - variable expression
1542  // cut - selection criteria
1543  // Output - return number of entries used to define variable
1544  // In addition mean, rms, median, and robust mean and rms (choosing fraction of data with smallest RMS)
1545  //
1546  /* Example usage:
1547  1.) create the robust estimators for variable expr="QA.TPC.CPass1.meanTPCncl" and create a corresponding
1548  aliases with the prefix alias[0]="ncl", calculated using fraction alias[1]="0.90"
1549 
1550  TStatToolkit::MakeStatAlias(tree,"QA.TPC.CPass1.meanTPCncl","QA.TPC.CPass1.status>0","ncl:0.9");
1551  root [4] tree->GetListOfAliases().Print()
1552  OBJ: TNamed ncl_Median (130.964333+0)
1553  OBJ: TNamed ncl_Mean (122.120387+0)
1554  OBJ: TNamed ncl_RMS (33.509623+0)
1555  OBJ: TNamed ncl_Mean90 (131.503862+0)
1556  OBJ: TNamed ncl_RMS90 (3.738260+0)
1557  */
1558  //
1559  Int_t entries = tree->Draw(expr,cut,"goff");
1560  if (entries<=1){
1561  printf("Expression or cut not valid:\t%s\t%s\n", expr, cut);
1562  return 0;
1563  }
1564  //
1565  TObjArray* oaAlias = TString(alias).Tokenize(":");
1566  if (oaAlias->GetEntries()<2) {
1567  printf("Alias must have 2 arguments:\t%s\n", alias);
1568  return 0;
1569  }
1570  Float_t entryFraction = atof( oaAlias->At(1)->GetName() );
1571  //
1572  Double_t median = TMath::Median(entries,tree->GetV1());
1573  Double_t mean = TMath::Mean(entries,tree->GetV1());
1574  Double_t rms = TMath::RMS(entries,tree->GetV1());
1575  Double_t meanEF=0, rmsEF=0;
1576  TStatToolkit::EvaluateUni(entries, tree->GetV1(), meanEF, rmsEF, entries*entryFraction);
1577  //
1578  tree->SetAlias(Form("%s_Median",oaAlias->At(0)->GetName()), Form("(%f+0)",median));
1579  tree->SetAlias(Form("%s_Mean",oaAlias->At(0)->GetName()), Form("(%f+0)",mean));
1580  tree->SetAlias(Form("%s_RMS",oaAlias->At(0)->GetName()), Form("(%f+0)",rms));
1581  tree->SetAlias(Form("%s_Mean%d",oaAlias->At(0)->GetName(),Int_t(entryFraction*100)), Form("(%f+0)",meanEF));
1582  tree->SetAlias(Form("%s_RMS%d",oaAlias->At(0)->GetName(),Int_t(entryFraction*100)), Form("(%f+0)",rmsEF));
1583  delete oaAlias;
1584  return entries;
1585 }
1586 
1587 Int_t TStatToolkit::SetStatusAlias(TTree * tree, const char * expr, const char * cut, const char * alias)
1588 {
1589  //
1590  // Add alias to trending tree using statistical values of a given variable.
1591  // (by MI, Patrick Reichelt)
1592  //
1593  // format of expr : varname (e.g. meanTPCncl)
1594  // format of cut : char like in TCut
1595  // format of alias: alias:query:entryFraction(EF) (fraction of entries used for uniformity evaluation)
1596  // e.g.: varname_Out:(abs(varname-meanEF)>6.*rmsEF):0.8
1597  // available internal variables are: 'varname, Median, Mean, MeanEF, RMS, RMSEF'
1598  // in the alias, 'varname' will be replaced by its content, and 'EF' by the percentage (e.g. MeanEF -> Mean80)
1599  //
1600  /* Example usage:
1601  1.) Define robust mean (possible, but easier done with TStatToolkit::MakeStatAlias(...))
1602  TStatToolkit::SetStatusAlias(tree, "meanTPCnclF", "meanTPCnclF>0", "meanTPCnclF_MeanEF:MeanEF:0.80") ;
1603  root [10] tree->GetListOfAliases()->Print()
1604  Collection name='TList', class='TList', size=1
1605  OBJ: TNamed meanTPCnclF_Mean80 0.899308
1606  2.) create alias outlyers - 6 sigma cut
1607  TStatToolkit::SetStatusAlias(tree, "meanTPCnclF", "meanTPCnclF>0", "meanTPCnclF_Out:(abs(meanTPCnclF-MeanEF)>6.*RMSEF):0.8")
1608  meanTPCnclF_Out ==> (abs(meanTPCnclF-0.899308)>6.*0.016590)
1609  3.) the same functionality as in 2.)
1610  TStatToolkit::SetStatusAlias(tree, "meanTPCnclF", "meanTPCnclF>0", "varname_Out2:(abs(varname-MeanEF)>6.*RMSEF):0.8")
1611  meanTPCnclF_Out2 ==> (abs(meanTPCnclF-0.899308)>6.*0.016590)
1612  */
1613  //
1614  Int_t entries = tree->Draw(expr,cut,"goff");
1615  if (entries<1){
1616  printf("Expression or cut not valid:\t%s\t%s\n", expr, cut);
1617  return 0;
1618  }
1619  //
1620  TObjArray* oaVar = TString(expr).Tokenize(":");
1621  char varname[50];
1622  snprintf(varname,50,"%s", oaVar->At(0)->GetName());
1623  Float_t entryFraction = 0.8;
1624  //
1625  TObjArray* oaAlias = TString(alias).Tokenize(":");
1626  if (oaAlias->GetEntries()<2) {
1627  printf("Alias must have at least 2 arguments:\t%s\n", alias);
1628  return 0;
1629  }
1630  else if (oaAlias->GetEntries()<3) {
1631  //printf("Using default entryFraction if needed:\t%f\n", entryFraction);
1632  }
1633  else entryFraction = atof( oaAlias->At(2)->GetName() );
1634  //
1635  Double_t median = TMath::Median(entries,tree->GetV1());
1636  Double_t mean = TMath::Mean(entries,tree->GetV1());
1637  Double_t rms = TMath::RMS(entries,tree->GetV1());
1638  Double_t meanEF=0, rmsEF=0;
1639  TStatToolkit::EvaluateUni(entries, tree->GetV1(), meanEF, rmsEF, entries*entryFraction);
1640  //
1641  TString sAlias( oaAlias->At(0)->GetName() );
1642  sAlias.ReplaceAll("varname",varname);
1643  sAlias.ReplaceAll("MeanEF", Form("Mean%1.0f",entryFraction*100) );
1644  sAlias.ReplaceAll("RMSEF", Form("RMS%1.0f",entryFraction*100) );
1645  TString sQuery( oaAlias->At(1)->GetName() );
1646  sQuery.ReplaceAll("varname",varname);
1647  sQuery.ReplaceAll("MeanEF", Form("%f",meanEF) );
1648  sQuery.ReplaceAll("RMSEF", Form("%f",rmsEF) ); //make sure to replace 'RMSEF' before 'RMS'...
1649  sQuery.ReplaceAll("Median", Form("%f",median) );
1650  sQuery.ReplaceAll("Mean", Form("%f",mean) );
1651  sQuery.ReplaceAll("RMS", Form("%f",rms) );
1652  printf("define alias:\t%s = %s\n", sAlias.Data(), sQuery.Data());
1653  //
1654  char query[200];
1655  char aname[200];
1656  snprintf(query,200,"%s", sQuery.Data());
1657  snprintf(aname,200,"%s", sAlias.Data());
1658  tree->SetAlias(aname, query);
1659  delete oaVar;
1660  delete oaAlias;
1661  return entries;
1662 }
1663 
1664 
1676 /*
1677  TString sTrendVars=";";
1678  // Ncl
1679  sTrendVars+="QA.TPC.meanTPCncl,TPC.Anchor.meanTPCncl,5,10,5;"; // delta Ncl warning 5 , error 10 (nominal ~ 100-140)
1680  sTrendVars+="QA.TPC.meanTPCnclF,TPC.Anchor.meanTPCnclF,0.02,0.05,0.05;"; // delta NclF warning 2%, error 5% (nominal ~ 90%)
1681  // dcaR resolution
1682  sTrendVars+="QA.TPC.dcarAP0,TPC.Anchor.dcarAP0,0.02,0.05,0.02;"; // dcarAP0; warning 0.02cm; error 0.05 cm (nominal ~ 0.2 cm)
1683 
1684 */
1685 void TStatToolkit::MakeAnchorAlias(TTree * tree, TString& sTrendVars, Int_t doCheck, Int_t verbose){
1686  const char* aType[4]={"Warning","Outlier","PhysAcc"};
1687  TObjArray *aTrendVars = sTrendVars.Tokenize(";");
1688  Int_t entries= aTrendVars->GetEntries();
1689  const Int_t kMaxEntries=100;
1690  for (Int_t ientry=0; ientry<entries; ientry++){
1691  TString variables[kMaxEntries];
1692  Bool_t isOK=kTRUE;
1693  TObjArray *descriptor=TString(aTrendVars->At(ientry)->GetName()).Tokenize(",");
1694  // check if valid sysntax
1695  if (descriptor->GetEntries()!=5){
1696  ::Error("makeAnchorAlias"," Invalid status descriptor. %s has %d members instead of 5 (variable, deltaWarning,deltaError, PhysAcc) ",aTrendVars->At(ientry)->GetName(),descriptor->GetEntries());
1697  continue;
1698  }
1699  // check individual variables
1700  for (Int_t ivar=0; ivar<5; ivar++){
1701  variables[ivar]=descriptor->At(ivar)->GetName();
1702  if ((doCheck&1)>0){ // check input formulas in case specified
1703  TTreeFormula *form=new TTreeFormula("dummy",descriptor->At(ivar)->GetName(),tree);
1704  if (form->GetTree()==NULL){
1705  isOK=kFALSE;
1706  ::Error("makeAnchorAlias"," Invalid element %d, %s in %s",ivar, descriptor->At(ivar)->GetName(),aTrendVars->At(ientry)->GetName());
1707  continue;
1708  }
1709  }
1710  }
1711  if (!isOK) continue;
1712  for (Int_t itype=0; itype<3; itype++) {
1713  TString aName = TString::Format("absDiff.%s_%s", variables[0].Data(), aType[itype]);
1714  TString aValue = "";
1715  if (itype < 2) {
1716  aValue = TString::Format("abs(%s-%s)>%s", variables[0].Data(), variables[1].Data(),
1717  variables[2 + itype].Data()); // Warning or Error
1718  } else {
1719  aValue = TString::Format("abs(%s-%s)<%s", variables[0].Data(), variables[1].Data(),
1720  variables[2 + itype].Data()); // Physics acceptable inside
1721  }
1722  tree->SetAlias(aName.Data(), aValue.Data());
1723  if ((doCheck & 2) > 0) {
1724  TTreeFormula *form = new TTreeFormula("dummy", aName.Data(), tree);
1725  if (form->GetTree() == NULL) {
1726  isOK = kFALSE;
1727  ::Error("makeAnchorAlias", "Alias not valid \t%s\t%s", aName.Data(), aValue.Data());
1728  }
1729  }
1730  if (verbose > 0) {
1731  ::Info("makeAnchorAlias", "SetAlias\t%s\t%s", aName.Data(), aValue.Data());
1732  }
1733  }
1734  TString aName = TString::Format("absDiff.%s_", variables[0].Data());
1735  tree->SetAlias((aName+"WarningBand").Data(), (TString("1.*")+variables[2+0]).Data());
1736  tree->SetAlias((aName+"OutlierBand").Data(), (TString("1.*")+variables[2+1]).Data());
1737  tree->SetAlias((aName+"PhysAccBand").Data(), (TString("1.*")+variables[2+2]).Data());
1738  if (verbose > 0) {
1739  ::Info("makeAnchorAlias", "SetAlias \t%s\t%s", (aName+"WarningBand").Data(), variables[2+0].Data());
1740  }
1741  delete descriptor;
1742  }
1743  delete aTrendVars;
1744 }
1745 
1759 void TStatToolkit::MakeCombinedAlias(TTree * tree, TString& sCombinedStatus, Bool_t doCheck, Int_t verbose ){
1760  const char* aType[3]={"_Warning","_Outlier","_PhysAcc"};
1761  TObjArray *aTrendStatus = sCombinedStatus.Tokenize(";");
1762  Int_t entries= aTrendStatus->GetEntries();
1763  const Int_t kMaxEntries=100;
1764  for (Int_t ientry=0; ientry<entries; ientry++){
1765  TString variables[kMaxEntries];
1766  Bool_t isOK=kTRUE;
1767  TObjArray *descriptor=TString(aTrendStatus->At(ientry)->GetName()).Tokenize(",");
1768  // check if valid sysntax
1769  Int_t nElems=descriptor->GetEntries();
1770  if (nElems<2){
1771  ::Error("makeCombinedAlias"," Invalid combined status descriptor. Too few entries %d in status %s", nElems, aTrendStatus->At(ientry)->GetName());
1772  continue;
1773  }
1774  // check individual variables
1775  TString aliasName[3], aliasContent[3];
1776  for (Int_t ivar=0; ivar<nElems; ivar++){
1777  if ((doCheck&1)>0 &&ivar>0) {
1778  TTreeFormula *form=new TTreeFormula("dummy",descriptor->At(ivar)->GetName(),tree);
1779  isOK=form->GetTree()!=NULL;
1780  }
1781  if (!isOK){
1782  isOK=kFALSE;
1783  ::Error("makeCombinedAlias"," Invalid element %d, %s in %s ",ivar, descriptor->At(ivar)->GetName(),aTrendStatus->At(ientry)->GetName());
1784  break;
1785  }else{
1786  variables[ivar]=descriptor->At(ivar)->GetName();
1787  for (Int_t iType=0; iType<3; iType++){ //type loop
1788  if (ivar==0) {
1789  aliasName[iType]=descriptor->At(ivar)->GetName();
1790  aliasName[iType]+=aType[iType];
1791  }
1792  else{
1793  if (ivar==1) aliasContent[iType]="(";
1794  aliasContent[iType]+=descriptor->At(ivar)->GetName();
1795  aliasContent[iType]+=aType[iType];
1796  if (ivar<nElems-1) aliasContent[iType]+="||";
1797  if (ivar==nElems-1) aliasContent[iType]+=")";
1798  }
1799  }
1800  }
1801  }
1802  if (isOK){
1803  for (Int_t iType=0; iType<3; iType++){
1804  tree->SetAlias(aliasName[iType].Data(), aliasContent[iType].Data());
1805  if (verbose>0){
1806  ::Info("makeCombinedAlias","SetAlias\t%s\t%s", aliasName[iType].Data(),aliasContent[iType].Data());
1807  }
1808  }
1809  }
1810  if (!isOK) continue;
1811  delete descriptor;
1812  }
1813  delete aTrendStatus;
1814 }
1815 
1816 
1852 
1853 TMultiGraph* TStatToolkit::MakeStatusMultGr(TTree * tree, const char * expr, const char * cut, const char * alias, Int_t igr)
1854 {
1855 
1856  TObjArray* oaVar = TString(expr).Tokenize(":");
1857  if (oaVar->GetEntries()<2) {
1858  printf("Expression has to be of type 'varname:xaxis':\t%s\n", expr);
1859  return 0;
1860  }
1861  char varname[50];
1862  char var_x[50];
1863  snprintf(varname,50,"%s", oaVar->At(0)->GetName());
1864  snprintf(var_x ,50,"%s", oaVar->At(1)->GetName());
1865  //
1866  TString sAlias(alias);
1867  sAlias.ReplaceAll("varname",varname);
1868  TObjArray* oaAlias = TString(sAlias.Data()).Tokenize(":");
1869  if (oaAlias->GetEntries()<2) {
1870  printf("Alias must have 2-6 arguments:\t%s\n", alias);
1871  return 0;
1872  }
1873  char query[200];
1874  TMultiGraph* multGr = new TMultiGraph();
1875  Int_t marArr[6] = {24+igr%2, 20+igr%2, 20+igr%2, 20+igr%2, 20+igr%2, 20+igr%2};
1876  Int_t colArr[6] = {kBlack, kBlack, kOrange, kRed, kGreen+1, kBlue};
1877  Double_t sizeArr[6] = {1.4, 1.1, 1.5, 1.1, 1.4, 0.8};
1878  Double_t shiftArr[6] = {0., 0., 0.25, 0.25, -0.25, -0.25};
1879  const Int_t ngr = oaAlias->GetEntriesFast();
1880  for (Int_t i=0; i<ngr; i++){
1881  snprintf(query,200, "%f*(%s-0.5):%s", 1.+igr, oaAlias->At(i)->GetName(), var_x);
1882  TGraphErrors * gr = (TGraphErrors*) TStatToolkit::MakeGraphSparse(tree,query,cut,marArr[i],colArr[i],sizeArr[i],shiftArr[i]);
1883  if (gr) {
1884  multGr->Add(gr);
1885  gr->SetName(oaAlias->At(i)->GetName());
1886  gr->SetTitle(oaAlias->At(i)->GetName());
1887  }
1888  else{
1889  ::Error("TStatToolkit::MakeStatusMultGr","MakeGraphSparse() returned with error when called with the following query: %s",query);
1890  return 0;
1891  }
1892  }
1893  //
1894  multGr->SetName(varname);
1895  multGr->SetTitle(varname); // used for y-axis labels of status bar, can be modified by calling function.
1896  delete oaVar;
1897  delete oaAlias;
1898  return multGr;
1899 }
1900 
1901 
1902 void TStatToolkit::AddStatusPad(TCanvas* c1, Float_t padratio, Float_t bottommargin)
1903 {
1904  //
1905  // add pad to bottom of canvas for Status graphs (by Patrick Reichelt)
1906  // call function "DrawStatusGraphs(...)" afterwards
1907  //
1908  TCanvas* c1_clone = (TCanvas*) c1->Clone("c1_clone");
1909  c1->Clear();
1910  // produce new pads
1911  c1->cd();
1912  TPad* pad1 = new TPad("pad1", "pad1", 0., padratio, 1., 1.);
1913  pad1->Draw();
1914  pad1->SetNumber(1); // so it can be called via "c1->cd(1);"
1915  c1->cd();
1916  TPad* pad2 = new TPad("pad2", "pad2", 0., 0., 1., padratio);
1917  pad2->Draw();
1918  pad2->SetNumber(2);
1919  // draw original canvas into first pad
1920  c1->cd(1);
1921  c1_clone->DrawClonePad();
1922  pad1->SetBottomMargin(0.001);
1923  pad1->SetRightMargin(0.01);
1924  // set up second pad
1925  c1->cd(2);
1926  pad2->SetGrid(3);
1927  pad2->SetTopMargin(0);
1928  pad2->SetBottomMargin(bottommargin); // for the long x-axis labels (runnumbers)
1929  pad2->SetRightMargin(0.01);
1930 }
1931 
1932 
1934 {
1935  //
1936  // draw Status graphs into active pad of canvas (by MI, Patrick Reichelt)
1937  // ...into bottom pad, if called after "AddStatusPad(...)"
1938  //
1939  const Int_t nvars = oaMultGr->GetEntriesFast();
1940  TGraph* grAxis = (TGraph*) ((TMultiGraph*) oaMultGr->At(0))->GetListOfGraphs()->At(0);
1941  grAxis->SetMaximum(0.5*nvars+0.5);
1942  grAxis->SetMinimum(0);
1943  grAxis->GetYaxis()->SetLabelSize(0);
1944  grAxis->GetYaxis()->SetTitle("");
1945  grAxis->SetTitle("");
1946  Int_t entries = grAxis->GetN();
1947  grAxis->GetXaxis()->SetLabelSize(5.7*TMath::Min(TMath::Max(5./entries,0.01),0.03));
1948  grAxis->GetXaxis()->LabelsOption("v");
1949  grAxis->Draw("ap");
1950  //
1951  // draw multigraphs & names of status variables on the y axis
1952  for (Int_t i=0; i<nvars; i++){
1953  ((TMultiGraph*) oaMultGr->At(i))->Draw("p");
1954  TLatex* ylabel = new TLatex(-0.1, 0.5*i+0.5, ((TMultiGraph*) oaMultGr->At(i))->GetTitle());
1955  ylabel->SetTextAlign(32); //hor:right & vert:centered
1956  ylabel->SetTextSize(0.025/gPad->GetHNDC());
1957  ylabel->Draw();
1958  }
1959 }
1960 
1961 
1962 TTree* TStatToolkit::WriteStatusToTree(TObject* oStatusGr)
1963 {
1964  //
1965  // Create Tree with Integers for each status variable flag (warning, outlier, physacc).
1966  // (by Patrick Reichelt)
1967  //
1968  // input: either a TMultiGraph with status of single variable, which
1969  // was computed by TStatToolkit::MakeStatusMultGr(),
1970  // or a TObjArray which contains up to 10 of such variables.
1971  // example: TTree* statusTree = WriteStatusToTree( TStatToolkit::MakeStatusMultGr(tree, "tpcItsMatch:run", "", sCriteria.Data(), 0) );
1972  // or : TTree* statusTree = TStatToolkit::WriteStatusToTree(oaMultGr);
1973  //
1974  // output tree: 1=flag is true, 0=flag is false, -1=flag was not computed.
1975  // To be rewritten to the pcstream
1976 
1977  TObjArray* oaMultGr = NULL;
1978  Bool_t needDeletion=kFALSE;
1979  if (oStatusGr->IsA() == TObjArray::Class()) {
1980  oaMultGr = (TObjArray*) oStatusGr;
1981  }
1982  else if (oStatusGr->IsA() == TMultiGraph::Class()) {
1983  oaMultGr = new TObjArray(); needDeletion=kTRUE;
1984  oaMultGr->Add((TMultiGraph*) oStatusGr);
1985  }
1986  else {
1987  Printf("WriteStatusToTree(): Error! 'oStatusGr' must be a TMultiGraph or a TObjArray of them!");
1988  return 0;
1989  }
1990  // variables for output tree
1991  const int nvarsMax=10;
1992  const int ncritMax=5;
1993  Int_t currentRun;
1994  Int_t treevars[nvarsMax*ncritMax];
1995  TString varnames[nvarsMax*ncritMax];
1996  for (int i=0; i<nvarsMax*ncritMax; i++) treevars[i]=-1;
1997 
1998  Printf("WriteStatusToTree(): writing following variables to TTree (maybe only subset of listed criteria filled)");
1999  for (Int_t vari=0; vari<nvarsMax; vari++)
2000  {
2001  if (vari < oaMultGr->GetEntriesFast()) {
2002  varnames[vari*ncritMax+0] = Form("%s_statisticOK", ((TMultiGraph*) oaMultGr->At(vari))->GetName());
2003  varnames[vari*ncritMax+1] = Form("%s_Warning", ((TMultiGraph*) oaMultGr->At(vari))->GetName());
2004  varnames[vari*ncritMax+2] = Form("%s_Outlier", ((TMultiGraph*) oaMultGr->At(vari))->GetName());
2005  varnames[vari*ncritMax+3] = Form("%s_PhysAcc", ((TMultiGraph*) oaMultGr->At(vari))->GetName());
2006  varnames[vari*ncritMax+4] = Form("%s_Extra", ((TMultiGraph*) oaMultGr->At(vari))->GetName());
2007  }
2008  else {
2009  varnames[vari*ncritMax+0] = Form("dummy");
2010  varnames[vari*ncritMax+1] = Form("dummy");
2011  varnames[vari*ncritMax+2] = Form("dummy");
2012  varnames[vari*ncritMax+3] = Form("dummy");
2013  varnames[vari*ncritMax+4] = Form("dummy");
2014  }
2015  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;
2016  }
2017 
2018  TTree* statusTree = new TTree("statusTree","statusTree");
2019  statusTree->Branch("run", &currentRun );
2020  statusTree->Branch(varnames[ 0].Data(), &treevars[ 0]);
2021  statusTree->Branch(varnames[ 1].Data(), &treevars[ 1]);
2022  statusTree->Branch(varnames[ 2].Data(), &treevars[ 2]);
2023  statusTree->Branch(varnames[ 3].Data(), &treevars[ 3]);
2024  statusTree->Branch(varnames[ 4].Data(), &treevars[ 4]);
2025  statusTree->Branch(varnames[ 5].Data(), &treevars[ 5]);
2026  statusTree->Branch(varnames[ 6].Data(), &treevars[ 6]);
2027  statusTree->Branch(varnames[ 7].Data(), &treevars[ 7]);
2028  statusTree->Branch(varnames[ 8].Data(), &treevars[ 8]);
2029  statusTree->Branch(varnames[ 9].Data(), &treevars[ 9]);
2030  statusTree->Branch(varnames[10].Data(), &treevars[10]);
2031  statusTree->Branch(varnames[11].Data(), &treevars[11]);
2032  statusTree->Branch(varnames[12].Data(), &treevars[12]);
2033  statusTree->Branch(varnames[13].Data(), &treevars[13]);
2034  statusTree->Branch(varnames[14].Data(), &treevars[14]);
2035  statusTree->Branch(varnames[15].Data(), &treevars[15]);
2036  statusTree->Branch(varnames[16].Data(), &treevars[16]);
2037  statusTree->Branch(varnames[17].Data(), &treevars[17]);
2038  statusTree->Branch(varnames[18].Data(), &treevars[18]);
2039  statusTree->Branch(varnames[19].Data(), &treevars[19]);
2040  statusTree->Branch(varnames[20].Data(), &treevars[20]);
2041  statusTree->Branch(varnames[21].Data(), &treevars[21]);
2042  statusTree->Branch(varnames[22].Data(), &treevars[22]);
2043  statusTree->Branch(varnames[23].Data(), &treevars[23]);
2044  statusTree->Branch(varnames[24].Data(), &treevars[24]);
2045  statusTree->Branch(varnames[25].Data(), &treevars[25]);
2046  statusTree->Branch(varnames[26].Data(), &treevars[26]);
2047  statusTree->Branch(varnames[27].Data(), &treevars[27]);
2048  statusTree->Branch(varnames[28].Data(), &treevars[28]);
2049  statusTree->Branch(varnames[29].Data(), &treevars[29]);
2050  statusTree->Branch(varnames[30].Data(), &treevars[30]);
2051  statusTree->Branch(varnames[31].Data(), &treevars[31]);
2052  statusTree->Branch(varnames[32].Data(), &treevars[32]);
2053  statusTree->Branch(varnames[33].Data(), &treevars[33]);
2054  statusTree->Branch(varnames[34].Data(), &treevars[34]);
2055  statusTree->Branch(varnames[35].Data(), &treevars[35]);
2056  statusTree->Branch(varnames[36].Data(), &treevars[36]);
2057  statusTree->Branch(varnames[37].Data(), &treevars[37]);
2058  statusTree->Branch(varnames[38].Data(), &treevars[38]);
2059  statusTree->Branch(varnames[39].Data(), &treevars[39]);
2060  statusTree->Branch(varnames[40].Data(), &treevars[40]);
2061  statusTree->Branch(varnames[41].Data(), &treevars[41]);
2062  statusTree->Branch(varnames[42].Data(), &treevars[42]);
2063  statusTree->Branch(varnames[43].Data(), &treevars[43]);
2064  statusTree->Branch(varnames[44].Data(), &treevars[44]);
2065  statusTree->Branch(varnames[45].Data(), &treevars[45]);
2066  statusTree->Branch(varnames[46].Data(), &treevars[46]);
2067  statusTree->Branch(varnames[47].Data(), &treevars[47]);
2068  statusTree->Branch(varnames[48].Data(), &treevars[48]);
2069  statusTree->Branch(varnames[49].Data(), &treevars[49]);
2070 
2071  // run loop
2072  Double_t graphX; // x-position of marker (0.5, 1.5, ...)
2073  Double_t graphY; // if >0 -> warning/outlier/physacc! if =-0.5 -> no warning/outlier/physacc
2074  TList* arrRuns = (TList*) ((TGraph*) ((TMultiGraph*) oaMultGr->At(0))->GetListOfGraphs()->At(0))->GetXaxis()->GetLabels();
2075  //'TAxis->GetLabels()' returns THashList of TObjString, but using THashList gives compilation error "... incomplete type 'struct THashList' "
2076  for (Int_t runi=0; runi<arrRuns->GetSize(); runi++)
2077  {
2078  currentRun = atoi( arrRuns->At(runi)->GetName() );
2079  //Printf(" runi=%2i, name: %s \t run number: %i", runi, arrRuns->At(runi)->GetName(), currentRun);
2080 
2081  // status variable loop
2082  for (Int_t vari=0; vari<oaMultGr->GetEntriesFast(); vari++)
2083  {
2084  TMultiGraph* multGr = (TMultiGraph*) oaMultGr->At(vari);
2085 
2086  // criteria loop
2087  // the order is given by TStatToolkit::MakeStatusMultGr().
2088  // criterion #1 is 'statisticOK' and mandatory, the rest is optional. (#0 is always True, thus skipped)
2089  for (Int_t criti=1; criti<multGr->GetListOfGraphs()->GetEntries(); criti++)
2090  {
2091  TGraph* grCriterion = (TGraph*) multGr->GetListOfGraphs()->At(criti);
2092  graphX = -1, graphY = -1;
2093  grCriterion->GetPoint(runi, graphX, graphY);
2094  treevars[(vari)*ncritMax+(criti-1)] = (graphY>0)?1:0;
2095  }
2096  }
2097  statusTree->Fill();
2098  }
2099 
2100  if (needDeletion) delete oaMultGr;
2101 
2102  return statusTree;
2103 }
2104 
2105 
2106 void TStatToolkit::MakeSummaryTree(TTree* treeIn, TTreeSRedirector *pcstream, TObjString & sumID, TCut &selection){
2107  //
2108  // Make a summary tree for the input tree
2109  // For the moment statistic works only for the primitive branches (Float/Double/Int)
2110  // Extension recursive version planned for graphs a and histograms
2111  //
2112  // Following statistics are exctracted:
2113  // - Standard: mean, meadian, rms
2114  // - LTM robust statistic: mean60, rms60, mean90, rms90
2115  // Parameters:
2116  // treeIn - input tree
2117  // pctream - Output redirector
2118  // sumID - ID as will be used in output tree
2119  // selection - selection criteria define the set of entries used to evaluat statistic
2120  //
2121  // Curently only predefined statistic used to fill summary information
2122  // Future plans an option user defined statistic descriptor instead of the defualt (if exist)
2123  //
2124  // e.g
2125  // default.Branches=median:mean90:rms90:mean60:rms60
2126  // interactionRate.Branches mean90:median:rms90:mean95:rms95:mean60:rms60
2127  //
2128  TObjArray * brArray = treeIn->GetListOfBranches();
2129  Int_t tEntries= treeIn->GetEntries();
2130  Int_t nBranches=brArray->GetEntries();
2131  TString treeName = treeIn->GetName();
2132 
2133  (*pcstream)<<treeName.Data()<<"entries="<<tEntries;
2134  (*pcstream)<<treeName.Data()<<"ID.="<<&sumID;
2135 
2136  TMatrixD valBranch(nBranches,7);
2137  for (Int_t iBr=0; iBr<nBranches; iBr++){
2138  TString brName= brArray->At(iBr)->GetName();
2139  Int_t entries=treeIn->Draw(TString::Format("%s>>dummy(10,0,1)",brArray->At(iBr)->GetName()).Data(),selection,"goff");
2140  if (entries==0) continue;
2141  Double_t median, mean, rms, mean60,rms60, mean90, rms90;
2142  mean = TMath::Mean(entries,treeIn->GetV1());
2143  median= TMath::Median(entries,treeIn->GetV1());
2144  rms = TMath::RMS(entries,treeIn->GetV1());
2145  TStatToolkit::EvaluateUni(entries, treeIn->GetV1(), mean60,rms60,TMath::Min(TMath::Max(2., 0.60*entries),Double_t(entries)));
2146  TStatToolkit::EvaluateUni(entries, treeIn->GetV1(), mean90,rms90,TMath::Min(TMath::Max(2., 0.90*entries),Double_t(entries)));
2147  valBranch(iBr,0)=mean;
2148  valBranch(iBr,1)=median;
2149  valBranch(iBr,2)=rms;
2150  valBranch(iBr,3)=mean60;
2151  valBranch(iBr,4)=rms60;
2152  valBranch(iBr,5)=mean90;
2153  valBranch(iBr,6)=rms90;
2154  (*pcstream)<<treeName.Data()<<
2155  brName+"="<<valBranch(iBr,1)<< // use as an default median estimator
2156  brName+"_Mean="<<valBranch(iBr,0)<<
2157  brName+"_Median="<<valBranch(iBr,1)<<
2158  brName+"_RMS="<<valBranch(iBr,2)<<
2159  brName+"_Mean60="<<valBranch(iBr,3)<<
2160  brName+"_RMS60="<<valBranch(iBr,4)<<
2161  brName+"_Mean90="<<valBranch(iBr,5)<<
2162  brName+"_RMS90="<<valBranch(iBr,6);
2163  }
2164  (*pcstream)<<treeName.Data()<<"\n";
2165 }
2166 
2167 
2168 
2169 TMultiGraph* TStatToolkit::MakeStatusLines(TTree * tree, const char * expr, const char * cut, const char * alias)
2170 {
2171  //
2172  // Create status lines for trending using MakeGraphSparse(), very similar to MakeStatusMultGr().
2173  // (by Patrick Reichelt)
2174  //
2175  // format of expr : varname:xaxis (e.g. meanTPCncl:run, but 'varname' can be any string that you need for seach-and-replace)
2176  // format of cut : char like in TCut
2177  // format of alias: varname_OutlierMin:varname_OutlierMax:varname_WarningMin:varname_WarningMax:varname_PhysAccMin:varname_PhysAccMax:varname_RobustMean
2178  //
2179  TObjArray* oaVar = TString(expr).Tokenize(":");
2180  if (oaVar->GetEntries()<2) {
2181  printf("Expression has to be of type 'varname:xaxis':\t%s\n", expr);
2182  return 0;
2183  }
2184  char varname[50];
2185  char var_x[50];
2186  snprintf(varname,50,"%s", oaVar->At(0)->GetName());
2187  snprintf(var_x ,50,"%s", oaVar->At(1)->GetName());
2188  //
2189  TString sAlias(alias);
2190  if (sAlias.IsNull()) { // alias for default usage set here:
2191  sAlias = "varname_OutlierMin:varname_OutlierMax:varname_WarningMin:varname_WarningMax:varname_PhysAccMin:varname_PhysAccMax:varname_RobustMean";
2192  }
2193  sAlias.ReplaceAll("varname",varname);
2194  TObjArray* oaAlias = TString(sAlias.Data()).Tokenize(":");
2195  if (oaAlias->GetEntries()<2) {
2196  printf("Alias must have 2-7 arguments:\t%s\n", alias);
2197  return 0;
2198  }
2199  char query[200];
2200  TMultiGraph* multGr = new TMultiGraph();
2201  Int_t colArr[7] = {kRed, kRed, kOrange, kOrange, kGreen+1, kGreen+1, kGray+2};
2202  const Int_t ngr = oaAlias->GetEntriesFast();
2203  for (Int_t i=0; i<ngr; i++){
2204  snprintf(query,200, "%s:%s", oaAlias->At(i)->GetName(), var_x);
2205  multGr->Add( (TGraphErrors*) TStatToolkit::MakeGraphSparse(tree,query,cut,29,colArr[i],1.5) );
2206  }
2207  //
2208  multGr->SetName(varname);
2209  multGr->SetTitle(varname);
2210  delete oaVar;
2211  delete oaAlias;
2212  return multGr;
2213 }
2214 
2215 
2216 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 )
2217 {
2218  //
2219  // Draw histogram from TTree with robust range
2220  // Only for 1D so far!
2221  //
2222  // Parameters:
2223  // - histoname: name of histogram
2224  // - histotitle: title of histgram
2225  // - fraction: fraction of data to define the robust mean
2226  // - nsigma: nsigma value for range
2227  //
2228  // To add:
2229  // automatic ranges - separatelly for X, Y and Z nbins - as string
2230  // names for the variables
2231  // option, entries, first entry like in tree draw
2232  //
2233  TString drawStr(drawCommand);
2234  TString cutStr(cuts);
2235  Int_t dim = 1;
2236 
2237  if(!tree) {
2238  ::Error("TStatToolkit::DrawHistogram","Tree pointer is NULL!");
2239  return 0;
2240  }
2241 
2242  // get entries
2243  Int_t entries = tree->Draw(drawStr.Data(), cutStr.Data(), "goff");
2244  if (entries == -1) {
2245  ::Error("TStatToolkit::DrawHistogram","Tree draw returns -!");
2246  return 0;
2247  }
2248  TObjArray *charray = drawStr.Tokenize(":");
2249 
2250  // get dimension
2251  if(tree->GetV1()) dim = 1;
2252  if(tree->GetV2()) dim = 2;
2253  if(tree->GetV3()) dim = 3;
2254  if(dim > 2){
2255  cerr<<"TTree has more than 2 dimensions (not yet supported)"<<endl;
2256  return 0;
2257  }
2258 
2259  // draw robust
2260  // Get estimators
2261  Double_t mean1=0, rms1=0, min1=0, max1=0;
2262  Double_t mean2=0, rms2=0, min2=0, max2=0;
2263  Double_t mean3=0, rms3=0, min3=0, max3=0;
2264 
2265  TStatToolkit::GetMinMaxMean( tree->GetV1(),entries, min1,max1, mean1);
2266  TStatToolkit::EvaluateUni(entries, tree->GetV1(),mean1,rms1, fraction*entries);
2267  if(dim>1){
2268  TStatToolkit::GetMinMaxMean( tree->GetV2(),entries, min2,max2, mean2);
2269  TStatToolkit::EvaluateUni(entries, tree->GetV1(),mean2,rms2, fraction*entries);
2270  }
2271  if(dim>2){
2272  TStatToolkit::GetMinMaxMean( tree->GetV3(),entries, min3,max3, mean3);
2273  TStatToolkit::EvaluateUni(entries, tree->GetV3(),mean3,rms3, fraction*entries);
2274  }
2275 
2276  TH1* hOut=NULL;
2277  if(dim==1){
2278  hOut = new TH1F(histoname, histotitle, 200, mean1-nsigma*rms1, mean1+nsigma*rms1);
2279  for (Int_t i=0; i<entries; i++) hOut->Fill(tree->GetV1()[i]);
2280  hOut->GetXaxis()->SetTitle(tree->GetHistogram()->GetXaxis()->GetTitle());
2281  }
2282  else if(dim==2){
2283  hOut = new TH2F(histoname, histotitle, 200, min2, max2,200, mean1-nsigma*rms1, mean1+nsigma*rms1);
2284  for (Int_t i=0; i<entries; i++) hOut->Fill(tree->GetV2()[i],tree->GetV1()[i]);
2285  hOut->GetXaxis()->SetTitle(tree->GetHistogram()->GetXaxis()->GetTitle());
2286  hOut->GetYaxis()->SetTitle(tree->GetHistogram()->GetYaxis()->GetTitle());
2287  }
2288  THashList * metaData = (THashList*) tree->GetUserInfo()->FindObject("metaTable");
2289 
2290  if (metaData != NULL){
2291  TNamed *nmdTitle0 = TStatToolkit::GetMetadata(tree,Form("%s.Title",charray->At(0)->GetName()));
2292  TNamed *nmdXAxis = TStatToolkit::GetMetadata(tree,Form("%s.AxisTitle",charray->At(1)->GetName()));
2293  TNamed *nmdTitle1 = TStatToolkit::GetMetadata(tree,Form("%s.Title",charray->At(1)->GetName()));
2294  TNamed *nmdYAxis = TStatToolkit::GetMetadata(tree,Form("%s.AxisTitle",charray->At(0)->GetName()));
2295  //
2296  TString hisTitle=charray->At(0)->GetName();
2297  if (nmdTitle0) hisTitle=nmdTitle0->GetTitle();
2298  if (nmdTitle1) {
2299  hisTitle+=":";
2300  hisTitle+=nmdTitle1->GetTitle();
2301  }else{
2302  hisTitle+=":";
2303  hisTitle+=charray->At(1)->GetName();
2304  }
2305  if (nmdYAxis) {hOut->GetYaxis()->SetTitle(nmdYAxis->GetTitle());}
2306  if (nmdXAxis) {hOut->GetXaxis()->SetTitle(nmdXAxis->GetTitle());}
2307  hOut->SetTitle(hisTitle);
2308  }
2309  delete charray;
2310  // if (option) hOut->Draw(option);
2311  return hOut;
2312 }
2313 
2314 void TStatToolkit::CheckTreeAliases(TTree * tree, Int_t ncheck){
2315  //
2316  // Check consistency of tree aliases
2317  //
2318  Int_t nCheck=100;
2319  TList * aliases = (TList*)tree->GetListOfAliases();
2320  Int_t entries = aliases->GetEntries();
2321  for (Int_t i=0; i<entries; i++){
2322  TObject * object= aliases->At(i);
2323  if (!object) continue;
2324  Int_t ndraw=tree->Draw(aliases->At(i)->GetName(),"1","goff",nCheck);
2325  if (ndraw==0){
2326  ::Error("Alias:\tProblem","%s",aliases->At(i)->GetName());
2327  }else{
2328  ::Info("Alias:\tOK","%s",aliases->At(i)->GetName());
2329  }
2330  }
2331 }
2332 
2333 
2334 
2335 
2336 Double_t TStatToolkit::GetDefaultStat(TTree * tree, const char * var, const char * selection, TStatType statType){
2337  //
2338  //
2339  //
2340  Int_t entries = tree->Draw(var,selection,"goff");
2341  if (entries==0) return 0;
2342  switch(statType){
2343  case kEntries:
2344  return entries;
2345  case kSum:
2346  return entries*TMath::Mean(entries, tree->GetV1());
2347  case kMean:
2348  return TMath::Mean(entries, tree->GetV1());
2349  case kRMS:
2350  return TMath::RMS(entries, tree->GetV1());
2351  case kMedian:
2352  return TMath::Median(entries, tree->GetV1());
2353  default:;
2354  }
2355  return 0;
2356 }
2357 
2358 //_____________________________________________________________________________
2359 void TStatToolkit::CombineArray(TTree *tree, TVectorD &values)
2360 {
2370  const Int_t numberOfDimensions = tree->GetPlayer()->GetDimension();
2371  if (numberOfDimensions==1) {
2372  values.Use(tree->GetSelectedRows(), tree->GetVal(0));
2373  return;
2374  }
2375 
2376  const Int_t numberOfSelectedRows = tree->GetSelectedRows();
2377  values.ResizeTo(numberOfDimensions * numberOfSelectedRows);
2378 
2379  Int_t nfill=0;
2380  for (Int_t idim=0; idim<numberOfDimensions; ++idim) {
2381  const Double_t *arr = tree->GetVal(idim);
2382  if (!arr) continue;
2383 
2384  for (Int_t ival=0; ival<numberOfSelectedRows; ++ival) {
2385  values.GetMatrixArray()[nfill++] = arr[ival];
2386  }
2387  }
2388 
2389 }
2390 
2391 //_____________________________________________________________________________
2392 Double_t TStatToolkit::GetDistance(const TVectorD &values, const ENormType normType,
2393  const Bool_t normaliseToEntries/*=kFALSE*/, const Double_t pvalue/*=1*/)
2394 {
2401 
2402  Double_t norm=0.;
2403 
2404  switch (normType) {
2405  case kL1:
2406  norm=values.Norm1();
2407  break;
2408  case kL2:
2409  norm=TMath::Sqrt(values.Norm2Sqr());
2410  break;
2411  case kLp:
2412  {
2413  if (pvalue<1.) {
2414  ::Error("TStatToolkit::GetDistance","Lp norm: p-value=%5.3g not valid. Only p-value>=1 is allowed", pvalue);
2415  break;
2416  }
2417  Double_t sum=0.;
2418  for (Int_t ival=0; ival<values.GetNrows(); ++ival) {
2419  sum+=TMath::Power(TMath::Abs(values.GetMatrixArray()[ival]), pvalue);
2420  }
2421  norm=TMath::Power(sum, 1./pvalue);
2422  }
2423  break;
2424  case kMax:
2425  norm=values.NormInf();
2426  break;
2427  case kHamming:
2428  {
2429  Double_t sum=0.;
2430  for (Int_t ival=0; ival<values.GetNrows(); ++ival) {
2431  if (TMath::Abs(values.GetMatrixArray()[ival])>1e-30) ++sum;
2432  }
2433  norm=sum;
2434  }
2435  break;
2436  default:
2437  return 0;
2438  }
2439  if (normaliseToEntries && values.GetNrows()>0) {
2440  norm/=values.GetNrows();
2441  }
2442  return norm;
2443 }
2444 
2445 //_____________________________________________________________________________
2446 Double_t TStatToolkit::GetDistance(const Int_t size, const Double_t *values, const ENormType normType,
2447  const Bool_t normaliseToEntries/*=kFALSE*/, const Double_t pvalue/*=1*/)
2448 {
2451  TVectorD vecvalues;
2452  vecvalues.Use(size, values);
2453  return GetDistance(vecvalues, normType, normaliseToEntries, pvalue);
2454 }
2455 
2456 //_____________________________________________________________________________
2457 Double_t TStatToolkit::GetDistance(TTree * tree, const char* var, const char * selection,
2458  const ENormType normType, const Bool_t normaliseToEntries/*=kFALSE*/, const Double_t pvalue/*=1*/)
2459 {
2475  Int_t entries = tree->Draw(var,selection,"goff");
2476  if (entries==0) return 0.;
2477 
2478  TVectorD values;
2479  CombineArray(tree, values);
2480  return GetDistance(values, normType, normaliseToEntries, pvalue);
2481 }
2482 
2483 
2484 
2485 void TStatToolkit::MakeDistortionMap(Int_t iter, THnBase * histo, TTreeSRedirector *pcstream, TMatrixD &projectionInfo,Int_t dumpHisto, Int_t verbose){
2486  //
2487  // Recursive function to calculate Distortion maps from the residual histograms
2488  // Input:
2489  // iter - ndim..0
2490  // histo - THn histogram
2491  // pcstream -
2492  // projectionInfo - TMatrix speicifiing distortion map cration setup
2493  // user specify columns:
2494  // 0.) sequence of dimensions
2495  // 1.) grouping in dimensions (how many bins will be groupd in specific dimension - 0 means onl specified bin 1, curren +-1 bin ...)
2496  // 2.) step in dimension ( in case >1 some n(projectionInfo(<dim>,2) bins will be not exported
2497  // internally used collumns (needed to pass current bin index and bin center to the recursive function)
2498  // 3.) current bin value
2499  // 4.) current bin center
2500  //
2501  // Output:
2502  // pcstream - file with output distortion tree
2503  // 1.) distortion characteristic: mean, rms, gaussian fit parameters, meang, rmsG chi2 ... at speciefied bin
2504  // 2.) specidfied bins (tree branches) are defined by the name of the histogram axis in input histograms
2505  //
2506  //
2507  // Example projection info
2508  /*
2509  TFile *f = TFile::Open("/hera/alice/hellbaer/alice-tpc-notes/SpaceChargeDistortion/data/ATO-108/fullMerge/SCcalibMergeLHC12d.root");
2510  THnF* histof= (THnF*) f->Get("deltaY_ClTPC_ITSTOF");
2511  histof->SetName("deltaRPhiTPCTISTOF");
2512  histof->GetAxis(4)->SetName("qpt");
2513  TH1::SetDirectory(0);
2514  TTreeSRedirector * pcstream = new TTreeSRedirector("distortion.root","recreate");
2515  TMatrixD projectionInfo(5,5);
2516  projectionInfo(0,0)=0; projectionInfo(0,1)=0; projectionInfo(0,2)=0;
2517  projectionInfo(1,0)=4; projectionInfo(1,1)=0; projectionInfo(1,2)=1;
2518  projectionInfo(2,0)=3; projectionInfo(2,1)=3; projectionInfo(2,2)=2;
2519  projectionInfo(3,0)=2; projectionInfo(3,1)=0; projectionInfo(3,2)=5;
2520  projectionInfo(4,0)=1; projectionInfo(4,1)=5; projectionInfo(4,2)=20;
2521  MakeDistortionMap(4, histof, pcstream, projectionInfo);
2522  delete pcstream;
2523  */
2524  //
2525  static TF1 fgaus("fgaus","gaus",-10,10);
2526  const Double_t kMinEntries=50;
2527  Int_t ndim=histo->GetNdimensions();
2528 // Int_t axis[ndim];
2529  Double_t meanVector[ndim];
2530  Int_t binVector[ndim];
2531  Double_t centerVector[ndim];
2532 // for (Int_t idim=0; idim<ndim; idim++) axis[idim]=idim;
2533  //
2534  if (iter==0){
2535  char tname[100];
2536  char aname[100];
2537  char bname[100];
2538  char cname[100];
2539 
2540  snprintf(tname, 100, "%sDist",histo->GetName());
2541  //
2542  //
2543  // 1.) Calculate properties - mean, rms, gaus fit, chi2, entries
2544  // 2.) Dump properties to tree 1D properties - plus dimension descriptor f
2545  Int_t dimProject = TMath::Nint(projectionInfo(iter,0));
2546  TH1 *his1DFull = histo->Projection(dimProject);
2547  Double_t mean= his1DFull->GetMean();
2548  Double_t rms= his1DFull->GetRMS();
2549  Int_t entries= his1DFull->GetEntries();
2550  TString hname="his_";
2551  for (Int_t idim=0; idim<ndim; idim++) {hname+="_"; hname+=TMath::Nint(projectionInfo(idim,3));}
2552  Double_t meanG=0, rmsG=0, chi2G=0;
2553  if (entries>kMinEntries){
2554  fgaus.SetParameters(entries,mean,rms);
2555  his1DFull->Fit(&fgaus,"qnr","qnr");
2556  meanG = fgaus.GetParameter(1);
2557  rmsG = fgaus.GetParameter(2);
2558  chi2G = fgaus.GetChisquare()/fgaus.GetNumberFreeParameters();
2559  }
2560  if (dumpHisto>=0) {
2561  static Int_t histoCounter=0;
2562  if ((histoCounter%dumpHisto)==0) his1DFull->Write(hname.Data());
2563  histoCounter++;
2564  }
2565  delete his1DFull;
2566  (*pcstream)<<tname<<
2567  "entries="<<entries<< // number of entries
2568  "mean="<<mean<< // mean value of the last dimension
2569  "rms="<<rms<< // rms value of the last dimension
2570  "meanG="<<meanG<< // mean of the gaus fit
2571  "rmsG="<<rmsG<< // rms of the gaus fit
2572  "chi2G="<<chi2G; // chi2 of the gaus fit
2573 
2574  for (Int_t idim=0; idim<ndim; idim++){
2575  TH1 *his1DAxis = histo->Projection(idim);
2576  meanVector[idim] = his1DAxis->GetMean();
2577  snprintf(aname, 100, "%sMean=",histo->GetAxis(idim)->GetName());
2578  (*pcstream)<<tname<<
2579  aname<<meanVector[idim]; // current bin means
2580  delete his1DAxis;
2581  }
2582  for (Int_t iIter=0; iIter<ndim; iIter++){
2583  Int_t idim = TMath::Nint(projectionInfo(iIter,0));
2584  binVector[idim] = TMath::Nint(projectionInfo(iIter,3));
2585  centerVector[idim] = projectionInfo(iIter,4);
2586  snprintf(bname, 100, "%sBin=",histo->GetAxis(idim)->GetName());
2587  snprintf(cname, 100, "%sCenter=",histo->GetAxis(idim)->GetName());
2588  (*pcstream)<<tname<<
2589  bname<<binVector[idim]<< // current bin values
2590  cname<<centerVector[idim]; // current bin centers
2591  }
2592  (*pcstream)<<tname<<"\n";
2593  }else{
2594  // loop over the diminsion of interest
2595  // project selecting bin+-deltabin histoProj
2596  // MakeDistortionMap(histoProj ...)
2597  //
2598  Int_t dimProject = TMath::Nint(projectionInfo(iter,0));
2599  Int_t groupProject = TMath::Nint(projectionInfo(iter,1));
2600  Int_t stepProject = TMath::Nint(projectionInfo(iter,2));
2601  if (stepProject<1) stepProject=1;
2602  Int_t nbins = histo->GetAxis(dimProject)->GetNbins();
2603 
2604  for (Int_t ibin=1; ibin<=nbins; ibin+=stepProject){
2605  if (iter>1 && verbose){
2606  for (Int_t idim=0; idim<ndim; idim++){
2607  printf("\t%d(%d,%d)",TMath::Nint(projectionInfo(idim,3)),TMath::Nint(projectionInfo(idim,0)),TMath::Nint(projectionInfo(idim,1) ));
2608  }
2609  printf("\n");
2610  }
2611  Int_t bin0=TMath::Max(ibin-groupProject,1);
2612  Int_t bin1=TMath::Min(ibin+groupProject,nbins);
2613  histo->GetAxis(dimProject)->SetRange(bin0,bin1);
2614  projectionInfo(iter,3)=ibin;
2615  projectionInfo(iter,4)=histo->GetAxis(dimProject)->GetBinCenter(ibin);
2616  Int_t iterProject=iter-1;
2617  MakeDistortionMap(iterProject, histo, pcstream, projectionInfo);
2618  }
2619  }
2620  //
2621 }
2622 
2623 void TStatToolkit::MakeDistortionMapFast(THnBase * histo, TTreeSRedirector *pcstream, TMatrixD &projectionInfo,Int_t verbose, Double_t fractionCut, const char * estimators)
2624 {
2625  //
2626  // Function to calculate Distortion maps from the residual histograms
2627  // Input:
2628  // histo - THn histogram
2629  // pcstream -
2630  // projectionInfo - TMatrix specifying distortion map creation setup
2631  // user specify columns:
2632  // 0.) sequence of dimensions
2633  // 1.) grouping in dimensions (how many bins will be grouped in specific dimension - 0 means onl specified bin 1, current +-1 bin ...)
2634  // 2.) step in dimension ( in case >1 some n(projectionInfo(<dim>,2) bins will be not exported
2635  //
2636  // Output:
2637  // pcstream - file with output distortion tree
2638  // 1.) distortion characteristic: mean, rms, gaussian fit parameters, meanG, rmsG chi2 ... at specified bin
2639  // 2.) specified bins (tree branches) are defined by the name of the histogram axis in input histograms
2640  // 3.) in debug mode - controlled by env variable "gDumpHistoFraction" fraction of histogram + fits dumped to the file
2641  //
2642  // Example projection info
2643  /*
2644  TFile *f = TFile::Open("/hera/alice/hellbaer/alice-tpc-notes/SpaceChargeDistortion/data/ATO-108/fullMerge/SCcalibMergeLHC12d.root");
2645  THnF* histoF= (THnF*) f->Get("deltaY_ClTPC_ITSTOF");
2646  histoF->SetName("deltaRPhiTPCTISTOF");
2647  histoF->GetAxis(4)->SetName("qpt");
2648  TH1::SetDirectory(0);
2649  TTreeSRedirector * pcstream = new TTreeSRedirector("distortion.root","recreate");
2650  TMatrixD projectionInfo(5,3);
2651  projectionInfo(0,0)=0; projectionInfo(0,1)=0; projectionInfo(0,2)=0;
2652  projectionInfo(1,0)=4; projectionInfo(1,1)=0; projectionInfo(1,2)=1;
2653  projectionInfo(2,0)=3; projectionInfo(2,1)=3; projectionInfo(2,2)=2;
2654  projectionInfo(3,0)=2; projectionInfo(3,1)=0; projectionInfo(3,2)=5;
2655  projectionInfo(4,0)=1; projectionInfo(4,1)=5; projectionInfo(4,2)=20;
2656  MakeDistortionMap(histoF, pcstream, projectionInfo);
2657  delete pcstream;
2658  */
2659  //
2660  if (histo->GetNdimensions()<=1) {
2661  ::Error("TStatToolkit::MakeDistortionMapFast","Invalid dimension of input histogram");
2662  return;
2663  }
2664  if (fractionCut<0 || fractionCut>0.4){
2665  ::Error("TStatToolkit::MakeDistortionMapFast","Invalid input fraction cut %f\r. Should be in range <0,0.4>", fractionCut);
2666  return;
2667  }
2668  const Double_t kMinEntries=30, kUseLLFrom=20;
2669  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
2670  char tname[100];
2671  char aname[100];
2672  char bname[100];
2673  char cname[100];
2674  Float_t fractionLTM[100]={0.8};
2675  TVectorF *vecLTM[100]={0};
2676  Int_t nestimators=1;
2677  if (estimators!=NULL){
2678  TObjArray * array=TString(estimators).Tokenize(":");
2679  nestimators=array->GetEntries();
2680  for (Int_t iest=0; iest<nestimators; iest++){
2681  fractionLTM[iest]=TString(array->At(iest)->GetName()).Atof();
2682  }
2683  }
2684  for (Int_t iest=0; iest<nestimators; iest++) {
2685  vecLTM[iest]=new TVectorF(10);
2686  (*(vecLTM[iest]))[9]= fractionLTM[iest];
2687  }
2688 
2689  //
2690  int ndim = histo->GetNdimensions();
2691  int nbins[ndim],idx[ndim],idxmin[ndim],idxmax[ndim],idxSav[ndim];
2692  for (int id=0;id<ndim;id++) nbins[id] = histo->GetAxis(id)->GetNbins();
2693  //
2694  int axOrd[ndim],binSt[ndim],binGr[ndim];
2695  for (int i=0;i<ndim;i++) {
2696  axOrd[i] = TMath::Nint(projectionInfo(i,0));
2697  binGr[i] = TMath::Nint(projectionInfo(i,1));
2698  binSt[i] = TMath::Max(1,TMath::Nint(projectionInfo(i,2)));
2699  }
2700  int tgtDim = axOrd[0],tgtStep=binSt[0],tgtNb=nbins[tgtDim],tgtNb1=tgtNb+1;
2701  double binY[tgtNb],binX[tgtNb],meanVector[ndim],centerVector[ndim], meanVectorMI[ndim+2];
2702  Int_t binVector[ndim];
2703  // prepare X axis
2704  TAxis* xax = histo->GetAxis(tgtDim);
2705  for (int i=tgtNb;i--;) binX[i] = xax->GetBinCenter(i+1);
2706  for (int i=ndim;i--;) idx[i]=1;
2707  Bool_t grpOn = kFALSE;
2708  for (int i=1;i<ndim;i++) if (binGr[i]) grpOn = kTRUE;
2709  //
2710  // estimate number of output fits
2711  histo->GetListOfAxes()->Print();
2712  ULong64_t nfits = 1, fitCount=0;
2713  printf("index\tdim\t|\tnbins\tgrouping\tstep\tnfits\n");
2714  for (int i=1;i<ndim;i++) {
2715  int idim = axOrd[i];
2716  nfits *= TMath::Max(1,nbins[idim]/binSt[idim]);
2717  printf("%d %d | %d %d %d %lld\n",i,idim,nbins[idim],binGr[idim], binSt[idim],nfits);
2718  }
2719  printf("Expect %lld nfits\n",nfits);
2720  ULong64_t fitProgress = nfits/100;
2721  //
2722  // setup fit function, at the moment full root fit
2723  static TF1 fgaus("fgaus","gaus",-10,10);
2724  fgaus.SetRange(xax->GetXmin(),xax->GetXmax());
2725  // TGraph grafFit(tgtNb);
2726  TH1F* hfit = new TH1F("hfit","hfit",tgtNb,xax->GetXmin(),xax->GetXmax());
2727  //
2728  snprintf(tname, 100, "%sDist",histo->GetName());
2729  TStopwatch sw;
2730  sw.Start();
2731  int dimVar=1, dimVarID = axOrd[dimVar];
2732  //
2733  // TVectorF vecLTM(9);
2734  while(1) {
2735  //
2736  //double dimVarCen = histo->GetAxis(dimVarID)->GetBinCenter(idx[dimVarID]); // center of currently varied bin
2737  //
2738  if (grpOn) { //>> grouping requested?
2739  memset(binY,0,tgtNb*sizeof(double)); // need to accumulate
2740  //
2741  for (int idim=1;idim<ndim;idim++) {
2742  int grp = binGr[idim];
2743  int idimR = axOrd[idim]; // real axis id
2744  idxSav[idimR]=idx[idimR]; // save central bins
2745  idxmax[idimR] = TMath::Min(idx[idimR]+grp,nbins[idimR]);
2746  idx[idimR] = idxmin[idimR] = TMath::Max(1,idx[idimR]-grp);
2747  //
2748  // effective bin center
2749  meanVector[idimR] = 0;
2750  TAxis* ax = histo->GetAxis(idimR);
2751  if (grp>0) {
2752  for (int k=idxmin[idimR];k<=idxmax[idimR];k++) meanVector[idimR] += ax->GetBinCenter(k);
2753  meanVector[idimR] /= (1+(grp<<1));
2754  }
2755  else meanVector[idimR] = ax->GetBinCenter(idxSav[idimR]);
2756  } // set limits for grouping
2757  if (verbose>0) {
2758  printf("output bin: ");
2759  for (int i=0;i<ndim;i++) if (i!=tgtDim) printf("[D:%d]:%d ",i,idxSav[i]);
2760  printf("\n");
2761  printf("integrates: ");
2762  for (int i=0;i<ndim;i++) if (i!=tgtDim) printf("[D:%d]:%d-%d ",i,idxmin[i],idxmax[i]);
2763  printf("\n");
2764  }
2765  //
2766  for (Int_t jDim=0; jDim<ndim+2; jDim++) meanVectorMI[jDim]=0;
2767  while(1) {
2768  // loop over target dimension: accumulation
2769  int &it = idx[tgtDim];
2770  for (it=1;it<tgtNb1;it+=tgtStep) {
2771  Double_t content=histo->GetBinContent(idx);
2772  binY[it-1] += content;
2773  for (Int_t jDim=0; jDim<ndim; jDim++) meanVectorMI[jDim]+=content*histo->GetAxis(jDim)->GetBinCenter(idx[jDim]);
2774  meanVectorMI[ndim]+=content;
2775  meanVectorMI[ndim+1]+=1.;
2776  if (verbose>1) {for (int i=0;i<ndim;i++) printf("%d ",idx[i]); printf(" | accumulation\n");}
2777  }
2778  //
2779  int idim;
2780  for (idim=1;idim<ndim;idim++) { // dimension being groupped
2781  int idimR = axOrd[idim]; // real axis id in the histo
2782  if ( (++idx[idimR]) > idxmax[idimR] ) idx[idimR]=idxmin[idimR];
2783  else break;
2784  }
2785  if (idim==ndim) break;
2786  }
2787  } // <<grouping requested
2788  else {
2789  int &it = idx[tgtDim];
2790  for (it=1;it<tgtNb1;it+=tgtStep) {
2791  binY[it-1] = histo->GetBinContent(idx);
2792  //
2793  //for (int i=0;i<ndim;i++) printf("%d ",idx[i]); printf(" | \n");
2794  }
2795  for (int idim=1;idim<ndim;idim++) {
2796  int idimR = axOrd[idim]; // real axis id
2797  meanVector[idimR] = histo->GetAxis(idimR)->GetBinCenter(idx[idimR]);
2798  }
2799  }
2800  if (grpOn) for (int i=ndim;i--;) idx[i]=idxSav[i]; // restore central bins
2801  idx[tgtDim] = 0;
2802  if (verbose>0) {for (int i=0;i<ndim;i++) printf("%d ",idx[i]); printf(" | central bin fit\n");}
2803  if (meanVectorMI[ndim]>1) { //TODO check valgrind
2804  for (Int_t jDim=0; jDim<ndim; jDim++) meanVectorMI[jDim]/= meanVectorMI[ndim];
2805  }else{
2806  for (Int_t jDim=0; jDim<ndim; jDim++) meanVectorMI[jDim]= meanVector[ndim];
2807  }
2808  //
2809  // >> ------------- do fit
2810  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;
2811  hfit->Reset();
2812  for (int ip=tgtNb;ip--;) {
2813  //grafFit.SetPoint(ip,binX[ip],binY[ip]);
2814  hfit->SetBinContent(ip+1,binY[ip]);
2815  nrm += binY[ip];
2816  mean += binX[ip]*binY[ip];
2817  mom2 += binX[ip]*binX[ip]*binY[ip];
2818  if (maxVal<binY[ip]) maxVal = binY[ip];
2819  }
2820  if (nrm>0) {
2821  mean /= nrm;
2822  mom2 /= nrm;
2823  rms = mom2 - mean*mean;
2824  rms = rms>0 ? TMath::Sqrt(rms):0;
2825  }
2826  mean0=mean;
2827  rms0=rms;
2828 
2829 
2830  Int_t nbins1D=hfit->GetNbinsX();
2831  Float_t binMedian=0;
2832  Double_t limits[2]={hfit->GetBinCenter(1), hfit->GetBinCenter(nbins1D)};
2833  if (nrm>5) {
2834  for (Int_t iest=0; iest<nestimators; iest++){
2835  TStatToolkit::LTMHisto(hfit, *(vecLTM[iest]), fractionLTM[iest]);
2836  }
2837  Double_t* integral=hfit->GetIntegral();
2838  for (Int_t i=1; i<nbins1D-1; i++){
2839  if (integral[i-1]<0.5 && integral[i]>=0.5){
2840  if (hfit->GetBinContent(i-1)+hfit->GetBinContent(i)>0){
2841  binMedian=hfit->GetBinCenter(i);
2842  Double_t dIdx=-(integral[i-1]-integral[i]);
2843  Double_t dx=(0.5+(0.5-integral[i])/dIdx)*hfit->GetBinWidth(i);
2844  binMedian+=dx;
2845  }
2846  }
2847  if (integral[i-1]<fractionCut && integral[i]>=fractionCut){
2848  limits[0]=hfit->GetBinCenter(i-1)-hfit->GetBinWidth(i);
2849  }
2850  if (integral[i]<1-fractionCut && integral[i+1]>=1-fractionCut){
2851  limits[1]=hfit->GetBinCenter(i+1)+hfit->GetBinWidth(i);
2852  }
2853  }
2854  }
2855  if (nrm>5&&fractionCut>0 &&rms>0) {
2856  hfit->GetXaxis()->SetRangeUser(limits[0], limits[1]);
2857  mean=hfit->GetMean();
2858  rms=hfit->GetRMS();
2859  if (nrm>0 && rms>0) {
2860  m3=hfit->GetSkewness();
2861  m4=hfit->GetKurtosis();
2862  }
2863  // fgaus.SetRange(limits[0]-rms, limits[1]+rms); //TODO - to fix bug in range limits
2864  }else{
2865  fgaus.SetRange(xax->GetXmin(),xax->GetXmax());
2866  }
2867 
2868 
2869  Bool_t isFitValid=kFALSE;
2870  if (nrm>=kMinEntries && rms>hfit->GetBinWidth(nbins1D)/TMath::Sqrt(12.)) {
2871  fgaus.SetParameters(nrm/(rms/hfit->GetBinWidth(nbins1D)),mean,rms);
2872  fgaus.SetParError(0,nrm/(rms/hfit->GetBinWidth(nbins1D)));
2873  fgaus.SetParError(1,rms);
2874  fgaus.SetParError(2,rms);
2875  //grafFit.Fit(&fgaus,/*maxVal<kUseLLFrom ? "qnrl":*/"qnr");
2876  TFitResultPtr fitPtr= hfit->Fit(&fgaus,maxVal<kUseLLFrom ? "qnrlS+":"qnrS+");
2877  //TFitResultPtr fitPtr= hfit->Fit(&fgaus,"qnrlS");
2878  entriesG = fgaus.GetParameter(0);
2879  meanG = fgaus.GetParameter(1);
2880  rmsG = fgaus.GetParameter(2);
2881  chi2G = fgaus.GetChisquare()/fgaus.GetNumberFreeParameters();
2882  TFitResult * result = fitPtr.Get();
2883  if (result!=NULL){
2884  isFitValid = result->IsValid();
2885  }
2886  //
2887  }
2888  TH1 * hDump=0;
2889  if (nrm>=kMinEntries&& kDumpHistoFraction>0 && (gRandom->Rndm()<kDumpHistoFraction || isFitValid!=kTRUE)){
2890  hDump=hfit;
2891  }
2892  if (kDumpHistoFraction>=1.){
2893  hDump=hfit;
2894  }
2895  if (hDump){
2896  hfit->GetListOfFunctions()->AddLast(&fgaus);
2897  (*pcstream)<<TString::Format("%sDump", tname).Data()<<
2898  "entries="<<nrm<< // number of entries
2899  "isFitValid="<<isFitValid<< // true if the gaus fit converged
2900  "hDump.="<<hDump<< // histogram - by default not filled
2901  "mean0="<<mean0<< // mean value of the last dimension - without fraction cut
2902  "rms0="<<rms0<< // rms value of the last dimension - without fraction cut
2903  "mean="<<mean<< // mean value of the last dimension
2904  "rms="<<rms<< // rms value of the last dimension
2905  "m3="<<m3<< // m3 (skewnes) of the last dimension
2906  "m4="<<m4<< // m4 (kurtosis) of the last dimension
2907  "binMedian="<<binMedian<< //binned median value of 1D histogram
2908  "entriesG="<<entriesG<<
2909  "meanG="<<meanG<< // mean of the gaus fit
2910  "rmsG="<<rmsG<< // rms of the gaus fit
2911  "vecLTM.="<<vecLTM[0]<< // LTM frac% statistic
2912  "chi2G="<<chi2G; // chi2 of the gaus fit
2913  for (Int_t iest=1; iest<nestimators; iest++)
2914  (*pcstream)<<TString::Format("%sDump", tname).Data()<<TString::Format("vecLTM%d.=",iest)<<vecLTM[iest]; // LTM frac% statistic
2915 
2916  }
2917 
2918  (*pcstream)<<tname<<
2919  "entries="<<nrm<< // number of entries
2920  "isFitValid="<<isFitValid<< // true if the gaus fit converged
2921  "mean0="<<mean0<< // mean value of the last dimension - without fraction cut
2922  "rms0="<<rms0<< // rms value of the last dimension - without fraction cut
2923  "mean="<<mean<< // mean value of the last dimension
2924  "rms="<<rms<< // rms value of the last dimension
2925  "m3="<<m3<< // m3 (skewnes) of the last dimension
2926  "m4="<<m4<< // m4 (kurtosis) of the last dimension
2927  "binMedian="<<binMedian<< //binned median value of 1D histogram
2928  "entriesG="<<entriesG<< //
2929  "meanG="<<meanG<< // mean of the gaus fit
2930  "rmsG="<<rmsG<< // rms of the gaus fit
2931  "vecLTM.="<<vecLTM[0]<< // LTM frac% statistic
2932  "chi2G="<<chi2G; // chi2 of the gaus fit
2933  for (Int_t iest=1; iest<nestimators; iest++)
2934  (*pcstream)<<tname<<TString::Format("vecLTM%d.=",iest)<<vecLTM[iest]; // LTM frac% statistic
2935 
2936 
2937  //
2938  meanVector[tgtDim] = mean; // what's a point of this?
2939  for (Int_t idim=0; idim<ndim; idim++){
2940  snprintf(aname, 100, "%sMean=",histo->GetAxis(idim)->GetName());
2941  (*pcstream)<<tname<<
2942  aname<<meanVectorMI[idim]; // current bin means
2943  }
2944  //
2945  for (Int_t iIter=0; iIter<ndim; iIter++){
2946  Int_t idim = axOrd[iIter];
2947  binVector[idim] = idx[idim];
2948  centerVector[idim] = histo->GetAxis(idim)->GetBinCenter(idx[idim]);
2949  snprintf(bname, 100, "%sBin=",histo->GetAxis(idim)->GetName());
2950  snprintf(cname, 100, "%sCenter=",histo->GetAxis(idim)->GetName());
2951  (*pcstream)<<tname<<
2952  bname<<binVector[idim]<< // current bin values
2953  cname<<centerVector[idim]; // current bin centers
2954  if (hDump){
2955  (*pcstream)<<TString::Format("%sDump", tname).Data()<<
2956  bname<<binVector[idim]<< // current bin values
2957  cname<<centerVector[idim]; // current bin centers
2958  }
2959  }
2960  (*pcstream)<<tname<<"\n";
2961  if (hDump) {
2962  (*pcstream)<<TString::Format("%sDump", tname).Data()<<"\n";
2963  hfit->GetListOfFunctions()->RemoveLast();
2964  }
2965  // << ------------- do fit
2966  //
2967  if ( fitProgress>0 && nfits>0) {
2968  if (((++fitCount)%fitProgress)==0) {
2969  printf("fit %lld %4.1f%% done\n",fitCount,100*double(fitCount)/nfits);
2970  }
2971  }
2972  //
2973  //next global bin in which target dimention will be looped
2974  for (dimVar=1;dimVar<ndim;dimVar++) { // varying dimension
2975  dimVarID = axOrd[dimVar]; // real axis id in the histo
2976  if ( (idx[dimVarID]+=binSt[dimVar]) > nbins[dimVarID] ) idx[dimVarID]=1;
2977  else break;
2978  }
2979  if (dimVar==ndim) break;
2980  }
2981  delete hfit;
2982  sw.Stop();
2983  sw.Print();
2984  /*
2985  int nb = histo->GetNbins();
2986  int prc = nb/100;
2987  for (int i=0;i<nb;i++) {
2988  histo->GetBinContent(i);
2989  if (i && (i%prc)==0) printf("Done %d%%\n",int(float(100*i)/prc));
2990  }
2991  */
2992 }
2993 
2994 
2999 
3010 void TStatToolkit::RebinSparseGraph(TGraph * graph0, TGraph *graph1, Option_t * option){
3011  if (graph0==NULL) throw std::invalid_argument( "RebinSparseGraph.graph0");
3012  if (graph1==NULL) throw std::invalid_argument( "RebinSparseGraph.graph1");
3013  TString opt = option;
3014  opt.ToLower();
3015  map<string,int> mapStrInt0, mapStrInt1;
3016  map<int,string> mapIntStr0, mapIntStr1;
3017  for (Int_t i=1; i<=graph0->GetXaxis()->GetNbins(); i++){
3018  mapStrInt0[graph0->GetXaxis()->GetBinLabel(i)]=i;
3019  mapIntStr0[i]=graph0->GetXaxis()->GetBinLabel(i);
3020  }
3021  for (Int_t i=1; i<=graph1->GetXaxis()->GetNbins(); i++){
3022  mapStrInt1[graph1->GetXaxis()->GetBinLabel(i)]=i;
3023  mapIntStr1[i]=graph1->GetXaxis()->GetBinLabel(i);
3024  }
3025  if (opt.Contains("merge")) {
3026  ::Error("Not supported option","%s",opt.Data());
3027  }
3028  for (Int_t i=1; i<=graph0->GetXaxis()->GetNbins(); i++){
3029  Int_t indexNew=mapStrInt1[mapIntStr0[i]];
3030  Double_t offset=graph0->GetX()[i-1]-int(graph0->GetX()[i-1]);
3031  graph0->GetX()[i-1]=indexNew+offset-1;
3032  }
3033  graph0->GetXaxis()->Set(graph1->GetXaxis()->GetNbins(),graph1->GetXaxis()->GetXbins()->GetArray());
3034  for (Int_t i=1; i<=graph0->GetXaxis()->GetNbins(); ++i){
3035  graph0->GetXaxis()->SetBinLabel(i,graph1->GetXaxis()->GetBinLabel(i));
3036  }
3037 }
3038 
3045 
3053 void TStatToolkit::RebinSparseMultiGraph(TMultiGraph *multiGraph, TGraph *graphRef){
3054  if (multiGraph==NULL) throw std::invalid_argument( "RebinSparseMultyGraph.multiGraph");
3055  if (multiGraph->GetListOfGraphs()==NULL) throw std::invalid_argument( "RebinSparseMultyGraph.multiGraph empty");
3056  if (graphRef==NULL) {
3057  graphRef=(TGraph*)multiGraph->GetListOfGraphs()->At(0);
3058  }
3059  for (Int_t i=0; i<multiGraph->GetListOfGraphs()->GetEntries(); i++){
3060  try {
3061  RebinSparseGraph((TGraph*)multiGraph->GetListOfGraphs()->At(i),graphRef);
3062  }catch(const std::invalid_argument& error){
3063  ::Error("RebinSparseMultiGraph","%s",error.what());
3064  }
3065  }
3066 }
3067 
3080 
3088 void TStatToolkit::MakeMultiGraphSparse(TMultiGraph *multiGraph) {
3089  //map<string,int> mapStrInt0, mapStrInt1;
3090  if (multiGraph == NULL) throw std::invalid_argument("MakeSparseMultiGraphInt.multiGraph");
3091  map<int, int> intCounter;
3092  map<int, int> intMap;
3093  vector<int> valueArray;
3094  const TList *grArray = multiGraph->GetListOfGraphs();
3095  for (Int_t iGr = 0; iGr < grArray->GetEntries(); ++iGr) {
3096  TGraph *iGraph = (TGraph *) grArray->At(iGr);
3097  for (Int_t iPoint = 0; iPoint < iGraph->GetN(); ++iPoint) {
3098  Int_t value = TMath::Nint(iGraph->GetX()[iPoint]);
3099  intCounter[value]++;
3100  }
3101  }
3102  for (std::map<int, int>::iterator iterator = intCounter.begin(); iterator != intCounter.end(); iterator++)
3103  valueArray.push_back(iterator->first);
3104  std::sort(valueArray.begin(), valueArray.begin());
3105  //stdsort(valueArray);
3106  for (UInt_t iValue = 0; iValue < valueArray.size(); iValue++) intMap[valueArray[iValue]] = iValue;
3107  //
3108  for (Int_t iGr = 0; iGr < grArray->GetEntries(); ++iGr) {
3109  TGraph *iGraph = (TGraph *) grArray->At(iGr);
3110  iGraph->GetXaxis()->Set(valueArray.size(), 0, valueArray.size());
3111  for (UInt_t iValue = 0; iValue < valueArray.size(); iValue++)
3112  iGraph->GetXaxis()->SetBinLabel(iValue + 1, TString::Format("%d", valueArray[iValue]).Data());
3113  for (Int_t iPoint = 0; iPoint < iGraph->GetN(); ++iPoint) {
3114  iGraph->GetX()[iPoint] = intMap[TMath::Nint(iGraph->GetX()[iPoint]) + 0.5];
3115  }
3116  }
3117  if (multiGraph->GetXaxis()) {
3118  multiGraph->GetXaxis()->Set(valueArray.size(), 0, valueArray.size());
3119  for (UInt_t iValue = 0; iValue < valueArray.size(); iValue++)
3120  multiGraph->GetXaxis()->SetBinLabel(iValue + 1, TString::Format("%d", valueArray[iValue]).Data());
3121  }
3122 }
3123 
3124 
3131 Int_t TStatToolkit::AdaptHistoMetadata(TTree* tree, TH1 *histogram, TString option){
3132  if (histogram==NULL) histogram= tree->GetHistogram();
3133  if (histogram==NULL) return -1;
3134  TNamed *named=0;
3135  named = TStatToolkit::GetMetadata(tree,TString(histogram->GetXaxis()->GetTitle())+".AxisTitle");
3136  if (named) histogram->GetXaxis()->SetTitle(named->GetTitle());
3137  named = TStatToolkit::GetMetadata(tree,TString(histogram->GetYaxis()->GetTitle())+".AxisTitle");
3138  if (named) histogram->GetYaxis()->SetTitle(named->GetTitle());
3139  if (histogram->GetZaxis()){
3140  named = TStatToolkit::GetMetadata(tree,TString(histogram->GetZaxis()->GetTitle())+".AxisTitle");
3141  if (named) histogram->GetZaxis()->SetTitle(named->GetTitle());
3142  }
3144  TPaletteAxis *palette = (TPaletteAxis*)histogram->GetListOfFunctions()->FindObject("palette");
3145  if (palette) palette->SetX2NDC(0.92);
3146  histogram->Draw(option.Data());
3147  return 1;
3148 }
3149 
3150 THashList * TStatToolkit::GetMetadata(TTree *tree) {
3151  TList *list=(TList*)(tree->GetUserInfo());
3152  return (THashList*)((list!=nullptr)? list->FindObject("metaTable"):0);
3153 }
Int_t colors[3]
Definition: CalibCosmic.C:62
TNamed * GetMetadata(TTree *tree, const char *vartagName, TString *prefix=0, Bool_t fullMatch=kFALSE)
printf("Chi2/npoints = %f\n", TMath::Sqrt(chi2/npoints))
TString FilterFit(const TString &input, const TString filter, TVectorD &vec, TMatrixD &covar)
void TestGausFit(Int_t nhistos=5000)
void DrawStatusGraphs(TObjArray *oaMultGr)
TMultiGraph * MakeStatusMultGr(TTree *tree, const char *expr, const char *cut, const char *alias, Int_t igr=0)
void MakeDistortionMapFast(THnBase *histo, TTreeSRedirector *pcstream, TMatrixD &projectionInfo, Int_t verbose=0, Double_t fractionCut=0.1, const char *estimators=0)
Double_t FitGaus(TH1 *his, TVectorT< T > *param=0, TMatrixT< T > *matrix=0, Float_t xmin=0, Float_t xmax=0, Bool_t verbose=kFALSE)
Definition: TStatToolkit.h:347
#define TObjArray
TVectorD vec
Definition: AnalyzeLaser.C:8
void Constrain1D(const TString &input, const TString filter, TVectorD &param, TMatrixD &covar, Double_t mean, Double_t sigma)
Double_t RobustBinMedian(TH1 *hist, Double_t fractionCut=1.)
void MakeCombinedAlias(TTree *tree, TString &sCombinedStatus, Bool_t doCheck, Int_t verbose)
void MakeDistortionMap(Int_t iter, THnBase *histo, TTreeSRedirector *pcstream, TMatrixD &projectionInfo, Int_t dumpHisto=100, Int_t verbose=kFALSE)
Int_t SetStatusAlias(TTree *tree, const char *expr, const char *cut, const char *alias)
TGraph2D * MakeStat2D(TH3 *his, Int_t delta0, Int_t delta1, Int_t type)
TH1 * DrawHistogram(TTree *tree, const char *drawCommand, const char *cuts="1", const char *hname="histo", const char *htitle="histo", Int_t nsigma=4, Float_t fraction=0.75, TObjArray *description=0)
void MakeMultiGraphSparse(TMultiGraph *multiGraph)
THashList * AddMetadata(TTree *, const char *vartagName, const char *varTagValue)
TMultiGraph * MakeStatusLines(TTree *tree, const char *expr, const char *cut, const char *alias)
TTreeSRedirector * pcstream
void MakeAnchorAlias(TTree *tree, TString &sTrendVars, Int_t doCheck, Int_t verbose)
strP3 Tokenize("+") -> Print()
TGraphErrors * MakeGraphErrors(TTree *tree, const char *expr="Entry", const char *cut="1", Int_t mstyle=25, Int_t mcolor=1, Float_t msize=-1, Float_t offset=0.0, Int_t entries=10000000, Int_t firstEntry=0)
npoints
Definition: driftITSTPC.C:85
AliSplineFit fit
Definition: CalibTime.C:30
void Update1D(Double_t delta, Double_t sigma, Int_t s1, TMatrixD &param, TMatrixD &covar)
TStatToolkit stat
Definition: AnalyzeLaser.C:5
Bool_t LTMHisto(TH1 *his, TVectorT< T > &param, Float_t fraction=1)
Definition: TStatToolkit.h:228
void LTM(TH1 *his, TVectorT< T > *param=0, Float_t fraction=1, Bool_t verbose=kFALSE)
Definition: TStatToolkit.h:444
TObjArray * array
Definition: AnalyzeLaser.C:12
void sum()
Int_t AdaptHistoMetadata(TTree *tree, TH1 *histogram, TString option)
Int_t Freq(Int_t n, const Int_t *inlist, Int_t *outlist, Bool_t down)
TTree * tree
Double_t chi2
Definition: AnalyzeLaser.C:7
void CheckTreeAliases(TTree *tree, Int_t ncheck)
static Int_t GetLineColor(const char *style, Int_t index)
char * prefix
TString * 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=-1, Int_t start=0, Int_t stop=10000000, Bool_t fix0=kFALSE)
TGraph * gr
Definition: CalibTime.C:25
TGraphErrors * MakeStat1D(TH2 *his, Int_t deltaBin, Double_t fraction, Int_t returnType, Int_t markerStyle, Int_t markerColor)
AliComparisonDraw comp
Definition: TestAnalisys.C:74
TString MakeFitString(const TString &input, const TVectorD &param, const TMatrixD &covar, Bool_t verbose=kFALSE)
Float_t range[5]
void EvaluateUni(Int_t nvectors, Double_t *data, Double_t &mean, Double_t &sigma, Int_t hh)
TTree * treeIn
TString * 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=-1, Int_t start=0, Int_t stop=10000000, Double_t constrain=-1)
void RebinSparseMultiGraph(TMultiGraph *multiGraph, TGraph *graphRef)
TTree * WriteStatusToTree(TObject *oStatusGr)
Int_t nBins
static Int_t GetLineStyle(const char *style, Int_t index)
TString * fitString
Definition: MakeGlobalFit.C:71
TVectorD errors
Definition: driftITSTPC.C:97
void CombineArray(TTree *tree, TVectorD &values)
void RebinSparseGraph(TGraph *graph0, TGraph *graph1, Option_t *option="")
gr Draw("same*")
void MakeSummaryTree(TTree *treeIn, TTreeSRedirector *pcstream, TObjString &sumID, TCut &selection)
TGraph * MakeGraphSparse(TTree *tree, const char *expr="Entry", const char *cut="1", Int_t mstyle=25, Int_t mcolor=1, Float_t msize=-1, Float_t offset=0.0)
Int_t grp(UInt_t run, TString &gdc)
Definition: onlineReco.C:70
Double_t GetDefaultStat(TTree *tree, const char *var, const char *selection, TStatType statType)
void GetMinMaxMean(const T *arr, Long64_t n, double &minVal, double &maxVal, double &meanVal)
Definition: TStatToolkit.h:180
static Int_t runNumber
Definition: pdc06_config.C:126
void EvaluateUniExternal(Int_t nvectors, Double_t *data, Double_t &mean, Double_t &sigma, Int_t hh, Float_t externalfactor=1)
void MakeDistortionMap()
void MedianFilter(TH1 *his1D, Int_t nmedian)
TCut cut
Definition: MakeGlobalFit.C:75
Float_t GetCOG(const Short_t *arr, Int_t nBins, Float_t xMin, Float_t xMax, Float_t *rms=0, Float_t *sum=0)
Double_t GetDistance(const TVectorD &values, const ENormType normType, const Bool_t normaliseToEntries=kFALSE, const Double_t pvalue=1.)
Int_t MakeStatAlias(TTree *tree, const char *expr, const char *cut, const char *alias)
static Int_t GetMarkerColor(const char *style, Int_t index)
void DrawMultiGraph(TMultiGraph *graph, Option_t *option)
static Int_t GetMarkerStyle(const char *style, Int_t index)
TMultiGraph * MakeMultGraph(TTree *tree, const char *groupName, const char *expr, const char *cut, const char *markers, const char *colors, Bool_t drawSparse, Float_t msize, Float_t sigmaRange, TLegend *legend, Bool_t comp=kTRUE)
TString * 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=-1, Int_t start=0, Int_t stop=10000000)
Int_t GetFitIndex(const TString fString, const TString subString)
void AddStatusPad(TCanvas *c1, Float_t padratio, Float_t bottommargin)