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