AliRoot Core  edcc906 (edcc906)
TStatToolkit.h
Go to the documentation of this file.
1 #ifndef TSTATTOOLKIT_H
2 #define TSTATTOOLKIT_H
3 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4  * See cxx source for full Copyright notice */
5 /* $Id$ */
6 
11 
12 #include "TMath.h"
13 #include "Riostream.h"
14 #include "TH1F.h"
15 #include "TH2F.h"
16 #include "TH3.h"
17 #include "THnBase.h"
18 #include "TF1.h"
19 #include "TTree.h"
20 #include "TChain.h"
21 #include "TObjString.h"
22 #include "TLinearFitter.h"
23 #include "TGraph2D.h"
24 #include "TGraph.h"
25 #include "TGraphErrors.h"
26 #include "TMultiGraph.h"
27 #include "TCanvas.h"
28 #include "TLatex.h"
29 #include "TCut.h"
30 #include "THashList.h"
31 #include "TFitResultPtr.h"
32 #include "TFitResult.h"
33 #include "TPaletteAxis.h"
34 //
35 // includes necessary for test functions
36 //
37 #include "TSystem.h"
38 #include "TRandom.h"
39 #include "TStopwatch.h"
40 #include "TTreeStream.h"
41 #include "AliSysInfo.h"
42 
43 #include "TObject.h"
44 #include "TVectorD.h"
45 #include "TVectorF.h"
46 #include "TMatrixD.h"
47 #include "TMatrixF.h"
48 #include <float.h>
49 //#include "TGraph2D.h"
50 //#include "TGraph.h"
51 class THashList;
52 
53 namespace TStatToolkit {
55  enum ENormType {kL1, kL2, kLp, kMax, kHamming, kNNormType }; // http://en.wikipedia.org/w/index.php?title=Norm_(mathematics)&oldid=655824636
56  //
57  //
58  void EvaluateUni(Int_t nvectors, Double_t *data, Double_t &mean, Double_t &sigma, Int_t hh);
59  void EvaluateUniExternal(Int_t nvectors, Double_t *data, Double_t &mean, Double_t &sigma, Int_t hh, Float_t externalfactor=1);
60  Int_t Freq(Int_t n, const Int_t *inlist, Int_t *outlist, Bool_t down);
61  //
62  // HISTOGRAMS TOOLS
63  //
64  template <typename T>
65  void TruncatedMean(const TH1 * his, TVectorT<T> *param, Float_t down=0, Float_t up=1.0, Bool_t verbose=kFALSE);
66  void MedianFilter(TH1 * his1D, Int_t nmedian);
67 
68  template <typename T> Bool_t LTMHisto(TH1 * his, TVectorT<T> &param , Float_t fraction=1);
69  template <typename T> Int_t* LTMUnbinned(int np, const T *arr, TVectorT<T> &params , Float_t keep=1.0);
70 
71  template <typename T> void Reorder(int np, T *arr, const int *idx);
72  //
73  template <typename T>
74  void LTM(TH1 * his, TVectorT<T> *param=0 , Float_t fraction=1, Bool_t verbose=kFALSE);
75 
76  template <typename T>
77  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);
78 
79  template <typename T>
80  Double_t FitGaus(Float_t *arr, Int_t nBins, Float_t xMin, Float_t xMax, TVectorT<T> *param=0, TMatrixT<T> *matrix=0, Bool_t verbose=kFALSE);
81  Float_t GetCOG(const Short_t *arr, Int_t nBins, Float_t xMin, Float_t xMax, Float_t *rms=0, Float_t *sum=0);
82 
83  TGraph2D * MakeStat2D(TH3 * his, Int_t delta0, Int_t delta1, Int_t type);
84  TGraphErrors * MakeStat1D(TH2 * his, Int_t deltaBin, Double_t fraction, Int_t returnType, Int_t markerStyle, Int_t markerColor);
85 
86  Double_t RobustBinMedian(TH1* hist, Double_t fractionCut=1.);
87  //
88  // Graph tools
89  //
90  THashList *AddMetadata(TTree*, const char *vartagName,const char *varTagValue);
91  TNamed *GetMetadata(TTree* tree, const char *vartagName, TString *prefix=0, Bool_t fullMatch=kFALSE);
92  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);
93  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);
94  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 );
95  void RebinSparseGraph(TGraph * graph0, TGraph *graph1, Option_t * option="");
96  void RebinSparseMultiGraph(TMultiGraph *multiGraph, TGraph *graphRef);
97  void MakeMultiGraphSparse(TMultiGraph *multiGraph);
98  //
99  // Fitting function
100  //
101  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);
102  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);
103  //
104  //Linear fitter helper function
105  //
106  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);
107  Int_t GetFitIndex(const TString fString, const TString subString);
108  TString FilterFit(const TString &input, const TString filter, TVectorD &vec, TMatrixD &covar);
109  void Update1D(Double_t delta, Double_t sigma, Int_t s1, TMatrixD &param, TMatrixD &covar);
110  void Constrain1D(const TString &input, const TString filter, TVectorD &param, TMatrixD & covar, Double_t mean, Double_t sigma);
111  TString MakeFitString(const TString &input, const TVectorD &param, const TMatrixD & covar, Bool_t verbose=kFALSE);
112  //
113  // TTree function for the trending
114  //
115  Int_t MakeStatAlias(TTree * tree, const char * expr, const char * cut, const char * alias);
116  void MakeAnchorAlias(TTree * tree, TString& sTrendVars, Int_t doCheck, Int_t verbose);
117  void MakeCombinedAlias(TTree * tree, TString& sCombinedStatus, Bool_t doCheck, Int_t verbose );
118  Int_t SetStatusAlias(TTree * tree, const char * expr, const char * cut, const char * alias);
119 
120  TMultiGraph* MakeStatusMultGr(TTree * tree, const char * expr, const char * cut, const char * alias, Int_t igr=0);
121  void AddStatusPad(TCanvas* c1, Float_t padratio, Float_t bottommargin);
122  void DrawMultiGraph(TMultiGraph *graph, Option_t * option); //
123  void DrawStatusGraphs(TObjArray* oaMultGr);
124  TTree* WriteStatusToTree(TObject* oStatusGr);
125  TMultiGraph* MakeStatusLines(TTree * tree, const char * expr, const char * cut, const char * alias);
126  void MakeSummaryTree(TTree* treeIn, TTreeSRedirector *pcstream, TObjString& sumID, TCut &selection);
127  Double_t GetDefaultStat(TTree * tree, const char * var, const char * selection, TStatType statType);
128  //
129  //
130  void MakeDistortionMap(Int_t iter, THnBase * histo, TTreeSRedirector *pcstream, TMatrixD &projectionInfo, Int_t dumpHisto=100,Int_t verbose=kFALSE);
131  void MakeDistortionMapFast(THnBase * histo, TTreeSRedirector *pcstream, TMatrixD &projectionInfo, Int_t verbose=0, Double_t fractionCut=0.1, const char * estimators=0);
132 
133  //
134  // norm (distance) functions
135  //
136  void CombineArray(TTree *tree, TVectorD &values);
137  Double_t GetDistance(const TVectorD &values, const ENormType normType,
138  const Bool_t normaliseToEntries=kFALSE, const Double_t pvalue=1.);
139  Double_t GetDistance(const Int_t size, const Double_t *values, const ENormType normType,
140  const Bool_t normaliseToEntries=kFALSE, const Double_t pvalue=1.);
141  Double_t GetDistance(TTree * tree, const char * var, const char * selection,
142  const ENormType normType, const Bool_t normaliseToEntries=kFALSE, const Double_t pvalue=1.);
143  //
144  // TTree function for robust draw
145  //
146  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 );
147  Int_t AdaptHistoMetadata(TTree* tree, TH1 *histogram, TString option);
148  //
149  // TestFunctions:
150  //
151  void TestGausFit(Int_t nhistos=5000);
152  void CheckTreeAliases(TTree * tree, Int_t ncheck);
153  //
154  // min, max, mean ...
155  template <typename T>
156  void GetMinMax(const T* arr, Long64_t n, double &minVal, double &maxVal);
157  template <typename T>
158  void GetMinMaxMean(const T* arr, Long64_t n, double &minVal, double &maxVal, double &meanVal);
159 
160 };
161 
162 //___________________________________________________________
163 template <typename T>
164 void TStatToolkit::GetMinMax(const T* arr, Long64_t n, double &minVal, double &maxVal)
165 {
166  // find min, max entries in the array in a single loop
167  minVal = DBL_MAX;
168  maxVal = -DBL_MAX;
169  for (int i=n;i--;) {
170  double val = arr[i];
171  if (val<minVal) minVal = val;
172  if (val>maxVal) maxVal = val;
173  }
174 }
175 
176 //___________________________________________________________
177 template <typename T>
178 void TStatToolkit::GetMinMaxMean(const T* arr, Long64_t n, double &minVal, double &maxVal, double &meanVal)
179 {
180  // find min, max entries in the array in a single loop
181  minVal = DBL_MAX;
182  maxVal = -DBL_MAX;
183  meanVal = 0;
184  for (int i=n;i--;) {
185  double val = arr[i];
186  if (val<minVal) minVal = val;
187  if (val>maxVal) maxVal = val;
188  meanVal += val;
189  }
190  if (n) meanVal /= n;
191 }
192 
193 //___TStatToolkit__________________________________________________________________________
194 template <typename T>
195 void TStatToolkit::TruncatedMean(const TH1 * his, TVectorT<T> *param, Float_t down, Float_t up, Bool_t verbose){
196  //
197  //
198  //
199  Int_t nbins = his->GetNbinsX();
200  Float_t nentries = his->GetEntries();
201  Float_t sum =0;
202  Float_t mean = 0;
203  Float_t sigma2 = 0;
204  Float_t ncumul=0;
205  for (Int_t ibin=1;ibin<nbins; ibin++){
206  ncumul+= his->GetBinContent(ibin);
207  Float_t fraction = Float_t(ncumul)/Float_t(nentries);
208  if (fraction>down && fraction<up){
209  sum+=his->GetBinContent(ibin);
210  mean+=his->GetBinCenter(ibin)*his->GetBinContent(ibin);
211  sigma2+=his->GetBinCenter(ibin)*his->GetBinCenter(ibin)*his->GetBinContent(ibin);
212  }
213  }
214  mean/=sum;
215  sigma2= TMath::Sqrt(TMath::Abs(sigma2/sum-mean*mean));
216  if (param){
217  (*param)[0] = his->GetMaximum();
218  (*param)[1] = mean;
219  (*param)[2] = sigma2;
220 
221  }
222  if (verbose) printf("Mean\t%f\t Sigma2\t%f\n", mean,sigma2);
223 }
224 
225 template <typename T>
226 Bool_t TStatToolkit::LTMHisto(TH1 *his1D, TVectorT<T> &params , Float_t fraction){
227  //
228  // LTM : Trimmed mean on histogram - Modified version for binned data
229  //
230  // Robust statistic to estimate properties of the distribution
231  // To handle binning error special treatment
232  // for definition of unbinned data see:
233  // http://en.wikipedia.org/w/index.php?title=Trimmed_estimator&oldid=582847999
234  //
235  // Function parameters:
236  // his1D - input histogram
237  // params - vector with parameters
238  // - 0 - area
239  // - 1 - mean
240  // - 2 - rms
241  // - 3 - error estimate of mean
242  // - 4 - error estimate of RMS
243  // - 5 - first accepted bin position
244  // - 6 - last accepted bin position
245  //
246  Int_t nbins = his1D->GetNbinsX();
247  Int_t nentries = (Int_t)his1D->GetEntries();
248  const Double_t kEpsilon=0.0000000001;
249 
250  if (nentries<=0) return 0;
251  if (fraction>1) fraction=0;
252  if (fraction<0) return 0;
253  TVectorD vectorX(nbins);
254  TVectorD vectorMean(nbins);
255  TVectorD vectorRMS(nbins);
256  Double_t sumCont=0;
257  for (Int_t ibin0=1; ibin0<=nbins; ibin0++) sumCont+=his1D->GetBinContent(ibin0);
258  //
259  Double_t minRMS=his1D->GetRMS()*10000;
260  Int_t maxBin=0;
261  //
262  for (Int_t ibin0=1; ibin0<nbins; ibin0++){
263  Double_t sum0=0, sum1=0, sum2=0;
264  Int_t ibin1=ibin0;
265  for ( ibin1=ibin0; ibin1<=nbins; ibin1++){
266  Double_t cont=his1D->GetBinContent(ibin1);
267  Double_t x= his1D->GetBinCenter(ibin1);
268  sum0+=cont;
269  sum1+=cont*x;
270  sum2+=cont*x*x;
271  if ( (ibin0!=ibin1) && sum0>=fraction*sumCont) break;
272  }
273  vectorX[ibin0]=his1D->GetBinCenter(ibin0);
274  if (sum0<fraction*sumCont) continue;
275  //
276  // substract fractions of bin0 and bin1 to keep sum0=fration*sumCont
277  //
278  Double_t diff = sum0-fraction*sumCont;
279  Double_t mean = (sum0>0) ? sum1/sum0:0;
280  //
281  Double_t x0=his1D->GetBinCenter(ibin0);
282  Double_t x1=his1D->GetBinCenter(ibin1);
283  Double_t y0=his1D->GetBinContent(ibin0);
284  Double_t y1=his1D->GetBinContent(ibin1);
285  //
286  Double_t d = y0+y1-diff; //enties to keep
287  Double_t w0=0,w1=0;
288  if (y0<=kEpsilon&&y1>kEpsilon){
289  w1=d/y1;
290  }
291  if (y1<=kEpsilon&&y0>kEpsilon){
292  w0=d/y0;
293  }
294  if (y0>kEpsilon && y1>kEpsilon && x1>x0 ){
295  w0 = (d*(x1-mean))/((x1-x0)*y0);
296  w1 = (d-y0*w0)/y1;
297  //
298  if (w0>1) {w1+=(w0-1)*y0/y1; w0=1;}
299  if (w1>1) {w0+=(w1-1)*y1/y0; w1=1;}
300  }
301  if ( (x1>x0) &&TMath::Abs(y0*w0+y1*w1-d)>kEpsilon*sum0){
302  printf(" TStatToolkit::LTMHisto error\n");
303  }
304  sum0-=y0+y1;
305  sum1-=y0*x0;
306  sum1-=y1*x1;
307  sum2-=y0*x0*x0;
308  sum2-=y1*x1*x1;
309  //
310  Double_t xx0=his1D->GetXaxis()->GetBinUpEdge(ibin0)-0.5*w0*his1D->GetBinWidth(ibin0);
311  Double_t xx1=his1D->GetXaxis()->GetBinLowEdge(ibin1)+0.5*w1*his1D->GetBinWidth(ibin1);
312  sum0+=y0*w0+y1*w1;
313  sum1+=y0*w0*xx0;
314  sum1+=y1*w1*xx1;
315  sum2+=y0*w0*xx0*xx0;
316  sum2+=y1*w1*xx1*xx1;
317 
318  //
319  // choose the bin with smallest rms
320  //
321  if (sum0>0){
322  vectorMean[ibin0]=sum1/sum0;
323  vectorRMS[ibin0]=TMath::Sqrt(TMath::Abs(sum2/sum0-vectorMean[ibin0]*vectorMean[ibin0]));
324  if (vectorRMS[ibin0]<minRMS){
325  minRMS=vectorRMS[ibin0];
326  params[0]=sum0;
327  params[1]=vectorMean[ibin0];
328  params[2]=vectorRMS[ibin0];
329  params[3]=vectorRMS[ibin0]/TMath::Sqrt(sumCont*fraction);
330  params[4]=0; // what is the formula for error of RMS???
331  params[5]=ibin0;
332  params[6]=ibin1;
333  params[7]=his1D->GetBinCenter(ibin0);
334  params[8]=his1D->GetBinCenter(ibin1);
335  maxBin=ibin0;
336  }
337  }else{
338  break;
339  }
340  }
341  return kTRUE;
342 }
343 
344 template <typename T>
345 Double_t TStatToolkit::FitGaus(TH1* his, TVectorT<T> *param, TMatrixT<T> */*matrix*/, Float_t xmin, Float_t xmax, Bool_t verbose){
346  //
347  // Fit histogram with gaussian function
348  //
349  // Prameters:
350  // return value- chi2 - if negative ( not enough points)
351  // his - input histogram
352  // param - vector with parameters
353  // xmin, xmax - range to fit - if xmin=xmax=0 - the full histogram range used
354  // Fitting:
355  // 1. Step - make logarithm
356  // 2. Linear fit (parabola) - more robust - always converge
357  // 3. In case of small statistic bins are averaged
358  //
359  static TLinearFitter fitter(3,"pol2");
360  TVectorD par(3);
361  TVectorD sigma(3);
362  TMatrixD mat(3,3);
363  if (his->GetMaximum()<4) return -1;
364  if (his->GetEntries()<12) return -1;
365  if (his->GetRMS()<mat.GetTol()) return -1;
366  Float_t maxEstimate = his->GetEntries()*his->GetBinWidth(1)/TMath::Sqrt((TMath::TwoPi()*his->GetRMS()));
367  Int_t dsmooth = TMath::Nint(6./TMath::Sqrt(maxEstimate));
368 
369  if (maxEstimate<1) return -1;
370  Int_t nbins = his->GetNbinsX();
371  Int_t npoints=0;
372  //
373 
374 
375  if (xmin>=xmax){
376  xmin = his->GetXaxis()->GetXmin();
377  xmax = his->GetXaxis()->GetXmax();
378  }
379  for (Int_t iter=0; iter<2; iter++){
380  fitter.ClearPoints();
381  npoints=0;
382  for (Int_t ibin=1;ibin<nbins+1; ibin++){
383  Int_t countB=1;
384  Float_t entriesI = his->GetBinContent(ibin);
385  for (Int_t delta = -dsmooth; delta<=dsmooth; delta++){
386  if (ibin+delta>1 &&ibin+delta<nbins-1){
387  entriesI += his->GetBinContent(ibin+delta);
388  countB++;
389  }
390  }
391  entriesI/=countB;
392  Double_t xcenter= his->GetBinCenter(ibin);
393  if (xcenter<xmin || xcenter>xmax) continue;
394  Double_t error=1./TMath::Sqrt(countB);
395  Float_t cont=2;
396  if (iter>0){
397  if (par[0]+par[1]*xcenter+par[2]*xcenter*xcenter>20) return 0;
398  cont = TMath::Exp(par[0]+par[1]*xcenter+par[2]*xcenter*xcenter);
399  if (cont>1.) error = 1./TMath::Sqrt(cont*Float_t(countB));
400  }
401  if (entriesI>1&&cont>1){
402  fitter.AddPoint(&xcenter,TMath::Log(Float_t(entriesI)),error);
403  npoints++;
404  }
405  }
406  if (npoints>3){
407  fitter.Eval();
408  fitter.GetParameters(par);
409  }else{
410  break;
411  }
412  }
413  if (npoints<=3){
414  return -1;
415  }
416  fitter.GetParameters(par);
417  fitter.GetCovarianceMatrix(mat);
418  if (TMath::Abs(par[1])<mat.GetTol()) return -1;
419  if (TMath::Abs(par[2])<mat.GetTol()) return -1;
420  Double_t chi2 = fitter.GetChisquare()/Float_t(npoints);
421  //fitter.GetParameters();
422  if (!param) param = new TVectorT<T>(3);
423  // if (!matrix) matrix = new TMatrixD(3,3); // Covariance matrix to be implemented
424  (*param)[1] = par[1]/(-2.*par[2]);
425  (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
426  (*param)[0] = TMath::Exp(par[0]+ par[1]* (*param)[1] + par[2]*(*param)[1]*(*param)[1]);
427  if (verbose){
428  par.Print();
429  mat.Print();
430  param->Print();
431  printf("Chi2=%f\n",chi2);
432  TF1 * f1= new TF1("f1","[0]*exp(-(x-[1])^2/(2*[2]*[2]))",his->GetXaxis()->GetXmin(),his->GetXaxis()->GetXmax());
433  f1->SetParameter(0, (*param)[0]);
434  f1->SetParameter(1, (*param)[1]);
435  f1->SetParameter(2, (*param)[2]);
436  f1->Draw("same");
437  }
438  return chi2;
439 }
440 
441 template <typename T>
442 void TStatToolkit::LTM(TH1 * his, TVectorT<T> *param , Float_t fraction, Bool_t verbose){
443  //
444  // LTM : Trimmed mean on histogram - Modified version for binned data
445  //
446  // Robust statistic to estimate properties of the distribution
447  // See http://en.wikipedia.org/w/index.php?title=Trimmed_estimator&oldid=582847999
448  //
449  // New faster version is under preparation
450  //
451  if (!param) return;
452  (*param)[0]=0;
453  (*param)[1]=0;
454  (*param)[2]=0;
455  Int_t nbins = his->GetNbinsX();
456  Int_t nentries = (Int_t)his->GetEntries();
457  if (nentries<=0) return;
458  Double_t *data = new Double_t[nentries];
459  Int_t npoints=0;
460  for (Int_t ibin=1;ibin<nbins; ibin++){
461  Double_t entriesI = his->GetBinContent(ibin);
462  //Double_t xcenter= his->GetBinCenter(ibin);
463  Double_t x0 = his->GetXaxis()->GetBinLowEdge(ibin);
464  Double_t w = his->GetXaxis()->GetBinWidth(ibin);
465  for (Int_t ic=0; ic<entriesI; ic++){
466  if (npoints<nentries){
467  data[npoints]= x0+w*Double_t((ic+0.5)/entriesI);
468  npoints++;
469  }
470  }
471  }
472  Double_t mean, sigma;
473  Int_t npoints2=TMath::Min(Int_t(fraction*Float_t(npoints)),npoints-1);
474  npoints2=TMath::Max(Int_t(0.5*Float_t(npoints)),npoints2);
475  TStatToolkit::EvaluateUni(npoints, data, mean,sigma,npoints2);
476  delete [] data;
477  if (verbose) printf("Mean\t%f\t Sigma2\t%f\n", mean,sigma);if (param){
478  (*param)[0] = his->GetMaximum();
479  (*param)[1] = mean;
480  (*param)[2] = sigma;
481  }
482 }
483 
484 
485 template <typename T>
486 Double_t TStatToolkit::FitGaus(Float_t *arr, Int_t nBins, Float_t xMin, Float_t xMax, TVectorT<T> *param, TMatrixT<T> */*matrix*/, Bool_t verbose){
487  //
488  // Fit histogram with gaussian function
489  //
490  // Prameters:
491  // nbins: size of the array and number of histogram bins
492  // xMin, xMax: histogram range
493  // param: paramters of the fit (0-Constant, 1-Mean, 2-Sigma)
494  // matrix: covariance matrix -- not implemented yet, pass dummy matrix!!!
495  //
496  // Return values:
497  // >0: the chi2 returned by TLinearFitter
498  // -3: only three points have been used for the calculation - no fitter was used
499  // -2: only two points have been used for the calculation - center of gravity was uesed for calculation
500  // -1: only one point has been used for the calculation - center of gravity was uesed for calculation
501  // -4: invalid result!!
502  //
503  // Fitting:
504  // 1. Step - make logarithm
505  // 2. Linear fit (parabola) - more robust - always converge
506  //
507  static TLinearFitter fitter(3,"pol2");
508  static TMatrixD mat(3,3);
509  static Double_t kTol = mat.GetTol();
510  fitter.StoreData(kFALSE);
511  fitter.ClearPoints();
512  TVectorD par(3);
513  TVectorD sigma(3);
514  TMatrixD matA(3,3);
515  TMatrixD b(3,1);
516  Float_t rms = TMath::RMS(nBins,arr);
517  Float_t max = TMath::MaxElement(nBins,arr);
518  Float_t binWidth = (xMax-xMin)/(Float_t)nBins;
519 
520  Float_t meanCOG = 0;
521  Float_t rms2COG = 0;
522  Float_t sumCOG = 0;
523 
524  Float_t entries = 0;
525  Int_t nfilled=0;
526 
527  for (Int_t i=0; i<nBins; i++){
528  entries+=arr[i];
529  if (arr[i]>0) nfilled++;
530  }
531 
532  if (max<4) return -4;
533  if (entries<12) return -4;
534  if (rms<kTol) return -4;
535 
536  Int_t npoints=0;
537  //
538 
539  //
540  for (Int_t ibin=0;ibin<nBins; ibin++){
541  Float_t entriesI = arr[ibin];
542  if (entriesI>1){
543  Double_t xcenter = xMin+(ibin+0.5)*binWidth;
544 
545  Float_t error = 1./TMath::Sqrt(entriesI);
546  Float_t val = TMath::Log(Float_t(entriesI));
547  fitter.AddPoint(&xcenter,val,error);
548  if (npoints<3){
549  matA(npoints,0)=1;
550  matA(npoints,1)=xcenter;
551  matA(npoints,2)=xcenter*xcenter;
552  b(npoints,0)=val;
553  meanCOG+=xcenter*entriesI;
554  rms2COG +=xcenter*entriesI*xcenter;
555  sumCOG +=entriesI;
556  }
557  npoints++;
558  }
559  }
560 
561 
562  Double_t chi2 = 0;
563  if (npoints>=3){
564  if ( npoints == 3 ){
565  //analytic calculation of the parameters for three points
566  matA.Invert();
567  TMatrixD res(1,3);
568  res.Mult(matA,b);
569  par[0]=res(0,0);
570  par[1]=res(0,1);
571  par[2]=res(0,2);
572  chi2 = -3.;
573  } else {
574  // use fitter for more than three points
575  fitter.Eval();
576  fitter.GetParameters(par);
577  fitter.GetCovarianceMatrix(mat);
578  chi2 = fitter.GetChisquare()/Float_t(npoints);
579  }
580  if (TMath::Abs(par[1])<kTol) return -4;
581  if (TMath::Abs(par[2])<kTol) return -4;
582 
583  if (!param) param = new TVectorT<T>(3);
584  //if (!matrix) matrix = new TMatrixT<T>(3,3); // !!!!might be a memory leek. use dummy matrix pointer to call this function! // Covariance matrix to be implemented
585 
586  (*param)[1] = par[1]/(-2.*par[2]);
587  (*param)[2] = 1./TMath::Sqrt(TMath::Abs(-2.*par[2]));
588  Double_t lnparam0 = par[0]+ par[1]* (*param)[1] + par[2]*(*param)[1]*(*param)[1];
589  if ( lnparam0>307 ) return -4;
590  (*param)[0] = TMath::Exp(lnparam0);
591  if (verbose){
592  par.Print();
593  mat.Print();
594  param->Print();
595  printf("Chi2=%f\n",chi2);
596  TF1 * f1= new TF1("f1","[0]*exp(-(x-[1])^2/(2*[2]*[2]))",xMin,xMax);
597  f1->SetParameter(0, (*param)[0]);
598  f1->SetParameter(1, (*param)[1]);
599  f1->SetParameter(2, (*param)[2]);
600  f1->Draw("same");
601  }
602  return chi2;
603  }
604 
605  if (npoints == 2){
606  //use center of gravity for 2 points
607  meanCOG/=sumCOG;
608  rms2COG /=sumCOG;
609  (*param)[0] = max;
610  (*param)[1] = meanCOG;
611  (*param)[2] = TMath::Sqrt(TMath::Abs(meanCOG*meanCOG-rms2COG));
612  chi2=-2.;
613  }
614  if ( npoints == 1 ){
615  meanCOG/=sumCOG;
616  (*param)[0] = max;
617  (*param)[1] = meanCOG;
618  (*param)[2] = binWidth/TMath::Sqrt(12);
619  chi2=-1.;
620  }
621  return chi2;
622 
623 }
624 
625 template <typename T>
626 Int_t* TStatToolkit::LTMUnbinned(int np, const T *arr, TVectorT<T> &params , Float_t keep)
627 {
628  //
629  // LTM : Trimmed mean of unbinned array
630  //
631  // Robust statistic to estimate properties of the distribution
632  // To handle binning error special treatment
633  // for definition of unbinned data see:
634  // http://en.wikipedia.org/w/index.php?title=Trimmed_estimator&oldid=582847999
635  //
636  // Function parameters:
637  // np - number of points in the array
638  // arr - data array (unsorted)
639  // params - vector with parameters
640  // - 0 - area
641  // - 1 - mean
642  // - 2 - rms
643  // - 3 - error estimate of mean
644  // - 4 - error estimate of RMS
645  // - 5 - first accepted element (of sorted array)
646  // - 6 - last accepted element (of sorted array)
647  //
648  // On success returns index of sorted events
649  //
650  static int book = 0;
651  static int *index = 0;
652  static float* w = 0;
653  int keepN = np*keep;
654  if (keepN>np) keepN = np;
655  if (keepN<2) return 0;
656  params[0] = 0.0f;
657  if (book<np) {
658  delete[] index;
659  book = np;
660  index = new int[book];
661  delete[] w;
662  w = new float[book+book];
663  }
664  //
665  float *wx1 = w, *wx2 = wx1+np;
666  TMath::Sort(np,arr,index,kFALSE); // sort in increasing order
667  // build cumulants
668  double sum1=0.0,sum2=0.0;
669  for (int i=0;i<np;i++) {
670  double x = arr[index[i]];
671  wx1[i] = (sum1+=x);
672  wx2[i] = (sum2+=x*x);
673  }
674  double minRMS = sum2+1e6;
675  params[0] = keepN;
676  int limI = np - keepN+1;
677  for (int i=0;i<limI;i++) {
678  int limJ = i+keepN-1;
679  Double_t sum1 = double(wx1[limJ]) - double(i ? wx1[i-1] : 0.0);
680  Double_t sum2 = double(wx2[limJ]) - double(i ? wx2[i-1] : 0.0);
681  double mean = sum1/keepN;
682  double rms2 = sum2/keepN - mean*mean;
683  // printf("%d : %d %e %e\n",i,limJ, mean, TMath::Sqrt(rms2));
684  if (rms2>minRMS) continue;
685  minRMS = rms2;
686  params[1] = mean;
687  params[2] = rms2;
688  params[5] = i;
689  params[6] = limJ;
690  }
691  //
692  if (!params[0]) return 0;
693  if (params[2]<0) {
694  ::Error("TStatToolkit::LTMUnbinned","rounding error in RMS<0");
695  // current fast algorithm N - Sometimes rounding errors can lead to negative RMS
696  // more precise algorithm N2
697  throw 1;
698  }
699  params[2] = TMath::Sqrt(params[2]);
700  params[3] = params[2]/TMath::Sqrt(params[0]); // error on mean
701  params[4] = params[3]/TMath::Sqrt(2.0); // error on RMS
702  return index;
703 }
704 
705 template <typename T>
706 void TStatToolkit::Reorder(int np, T *arr, const int *idx)
707 {
708  // rearange arr in order given by idx
709  const int kMaxOnStack = 10000;
710  // don't abuse stack
711  T *arrCHeap=0, arrCstack[np<kMaxOnStack ? np:1], *arrC=np<kMaxOnStack ? &arrCstack[0] : (arrCHeap=new T[np]);
712  memcpy(arrC,arr,np*sizeof(T));
713  for (int i=np;i--;) arr[i] = arrC[idx[i]];
714  delete[] arrCHeap;
715  //
716 }
717 
718 #endif
TBrowser b
Definition: RunAnaESD.C:12
Int_t colors[3]
Definition: CalibCosmic.C:62
void GetMinMax(const T *arr, Long64_t n, double &minVal, double &maxVal)
Definition: TStatToolkit.h:164
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)
const Double_t kEpsilon
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)
Summary of statistics functions.
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:345
#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)
TMatrixD mat
Definition: AnalyzeLaser.C:9
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)
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)
TTree * tree
npoints
Definition: driftITSTPC.C:85
void Update1D(Double_t delta, Double_t sigma, Int_t s1, TMatrixD &param, TMatrixD &covar)
Bool_t LTMHisto(TH1 *his, TVectorT< T > &param, Float_t fraction=1)
Definition: TStatToolkit.h:226
void LTM(TH1 *his, TVectorT< T > *param=0, Float_t fraction=1, Bool_t verbose=kFALSE)
Definition: TStatToolkit.h:442
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)
Double_t chi2
Definition: AnalyzeLaser.C:7
void CheckTreeAliases(TTree *tree, Int_t ncheck)
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)
TGraphErrors * MakeStat1D(TH2 *his, Int_t deltaBin, Double_t fraction, Int_t returnType, Int_t markerStyle, Int_t markerColor)
void TruncatedMean(const TH1 *his, TVectorT< T > *param, Float_t down=0, Float_t up=1.0, Bool_t verbose=kFALSE)
Definition: TStatToolkit.h:195
AliComparisonDraw comp
Definition: TestAnalisys.C:74
TString MakeFitString(const TString &input, const TVectorD &param, const TMatrixD &covar, Bool_t verbose=kFALSE)
void Reorder(int np, T *arr, const int *idx)
Definition: TStatToolkit.h:706
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
void CombineArray(TTree *tree, TVectorD &values)
void RebinSparseGraph(TGraph *graph0, TGraph *graph1, Option_t *option="")
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)
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:178
void EvaluateUniExternal(Int_t nvectors, Double_t *data, Double_t &mean, Double_t &sigma, Int_t hh, Float_t externalfactor=1)
class TVectorT< Double_t > TVectorD
void MedianFilter(TH1 *his1D, Int_t nmedian)
Int_t * LTMUnbinned(int np, const T *arr, TVectorT< T > &params, Float_t keep=1.0)
Definition: TStatToolkit.h:626
TCut cut
Definition: MakeGlobalFit.C:75
void res(Char_t i)
Definition: Resolution.C:2
Float_t GetCOG(const Short_t *arr, Int_t nBins, Float_t xMin, Float_t xMax, Float_t *rms=0, Float_t *sum=0)
const Int_t markers[5]
Definition: makeBBFit.C:58
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)
void DrawMultiGraph(TMultiGraph *graph, Option_t *option)
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)