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