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