Update master to aliroot
[u/mrichter/AliRoot.git] / STAT / Macros / TStatToolkitTest.C
CommitLineData
8ecd39c2 1/*
2 .L $ALICE_ROOT/STAT/Macros/TStatToolkitTest.C+
3
4*/
5
6#include "TH1.h"
7#include "TF1.h"
8#include "TMath.h"
9#include "TRandom.h"
10#include "TVectorD.h"
11#include "TStatToolkit.h"
12#include "TTreeStream.h"
13#include "TLegend.h"
14#include "TGraphErrors.h"
15#include "TCut.h"
16#include "TCanvas.h"
71be42ab 17#include "TLinearFitter.h"
8ecd39c2 18
19TObjArray arrayFit(3);
71be42ab 20Int_t kMarkers[10]={25,24,20,21,22};
21Int_t kColors[10]={1,2,4,3,6};
22const char * names[10] = {"LTM","LThisto","LTMhistoPar","Gaus fit", "Gaus in range"};
23const char * distNames[10] = {"Gauss","Gauss+flat","Gauss(m)+0.2*Gauss(m+5#sigma)","Gauss(m)+0.4*Gauss(m+5#sigma)"};
8ecd39c2 24
25
26void TestLTM(Int_t nevents=10000){
27 //
28 // Goal test and benchamerk numerical stability and precission of the LTM method
29 // Binned data and not binned data:
30 // Here we assume that binwidth<sigma
31 // Distribution types examples used in test:
32 //
33 // 0 - gaussian
34 // 1 - gauss+uniform background
35 // 2 - gauss+second gaus 5 sigma away
36 //
37 Int_t npointsMax=1000;
38 TTreeSRedirector * pcstream = new TTreeSRedirector("TStatToolkit_LTMtest.root","recreate");
39 //
40 TVectorD values(npointsMax);
41 TVectorD vecLTM(6);
42 TVectorD meanLTM(6);
43 TVectorD sigmaLTM(6);
44 TVectorD meanLTMHisto(6);
45 TVectorD sigmaLTMHisto(6);
46 TVectorD meanLTMHisto1(6);
47 TVectorD sigmaLTMHisto1(6);
71be42ab 48 TVectorD meanLTMHistoPar(6);
49 TVectorD sigmaLTMHistoPar(6);
50 TVectorD meanLTMHistoGausLim(6);
51 TVectorD sigmaLTMHistoGausLim(6);
8ecd39c2 52 TVectorD paramLTM(10);
53 //
54 for (Int_t iltm=0; iltm<6; iltm++) vecLTM[iltm]=0.99999*(0.5+iltm/10.);
55 TF1 fg("fg","gaus");
56 for (Int_t ievent = 0; ievent<nevents; ievent++){
57 if (ievent%1000==0) printf("%d\n",ievent);
58 Int_t distType=Int_t(gRandom->Rndm()*4);
59 Int_t npoints= 50+(npointsMax-50)*gRandom->Rndm();
60 Double_t mean = 0.5+(gRandom->Rndm()-0.5)*0.2;
61 Double_t sigma = mean*0.2*(1+2*gRandom->Rndm());
62 TH1F histo("histo","histo",100,-0.5,1.5);
63 //
64 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
65 Double_t value=0;
66 if (distType==0) value=gRandom->Gaus(mean,sigma); // gauss
67 if (distType==1) {
68 if (gRandom->Rndm()>0.2) { // gauss + background
69 value=gRandom->Gaus(mean,sigma);
70 }else{
71 value=mean+(gRandom->Rndm()-0.5)*2;
72 }
73 }
74 if (distType==2) {
75 if (gRandom->Rndm()>0.2) { // gauss + second gaus 5 sigma away
76 value=gRandom->Gaus(mean,sigma);
77 }else{
78 value=gRandom->Gaus(mean+5*sigma,sigma);
79 }
80 }
81 if (distType==3) {
71be42ab 82 if (gRandom->Rndm()>0.4) { // gauss + second gaus 5 sigma away
8ecd39c2 83 value=gRandom->Gaus(mean,sigma);
84 }else{
85 value=gRandom->Gaus(mean+5*sigma,sigma);
86 }
87 }
88 values[ipoint]=value;
89 histo.Fill(value);
90 }
91 //
92 histo.Fit(&fg,"QN","QN");
93 Double_t meanG = fg.GetParameter(1);
94 Double_t rmsG = fg.GetParameter(2);
95 Double_t meanA = TMath::Mean(npoints,values.GetMatrixArray());
96 Double_t rmsA = TMath::Mean(npoints,values.GetMatrixArray());
97 Double_t meanH = histo.GetMean();
98 Double_t rmsH = histo.GetRMS();
99 //
100 for (Int_t iltm=0; iltm<6; iltm++){
101 //
102 Double_t meanV,sigmaV=0;
103 TStatToolkit::EvaluateUni(npoints,values.GetMatrixArray(), meanV,sigmaV, vecLTM[iltm]*npoints);
104 meanLTM[iltm]=meanV;
105 sigmaLTM[iltm]=sigmaV;
106 //
107 TStatToolkit::LTM(&histo, &paramLTM, vecLTM[iltm]);
108 meanLTMHisto1[iltm]=paramLTM[1];
109 sigmaLTMHisto1[iltm]=paramLTM[2];
110 //
111 TStatToolkit::LTMHisto(&histo, paramLTM, vecLTM[iltm]);
112 meanLTMHisto[iltm]=paramLTM[1];
113 sigmaLTMHisto[iltm]=paramLTM[2];
71be42ab 114 //
115 // graph fit
116 //
117 Int_t binC=histo.GetXaxis()->FindBin(paramLTM[1]);
118 Int_t bin0=TMath::Max(histo.GetXaxis()->FindBin(paramLTM[1]-5*paramLTM[2]),1);
119 Int_t bin1=TMath::Min(histo.GetXaxis()->FindBin(paramLTM[1]+5*paramLTM[2]),histo.GetXaxis()->GetNbins()) ;
120 if ( (bin0-bin1) <3){
121 bin0=TMath::Max(binC-2,1);
122 bin1=TMath::Min(binC+2, histo.GetXaxis()->GetNbins());
123 }
124 //
125 /* test input:
126 TH1F histo("aaa","aaa",100,-10,10);for (Int_t i=0; i<10000; i++) histo.Fill(gRandom->Gaus(2,1));
127 */
128
129 TLinearFitter fitter(3,"hyp2");
130 for (Int_t ibin=bin0; ibin<bin1; ibin++){
131 Double_t x=histo.GetXaxis()->GetBinCenter(ibin);
132 Double_t y=histo.GetBinContent(ibin);
133 Double_t sy=histo.GetBinError(ibin);
134 Double_t xxx[3]={x,x*x};
135 if (y>3.*sy){
136 fitter.AddPoint(xxx,TMath::Log(y),sy/y);
137 }
138 }
139 if (fitter.GetNpoints()>3){
140 fitter.Eval();
141 // f(x)=[0]+[1]*x+[2]*x*x;
142 // f(x)=s*(x0-x0)^2+s
143 Double_t parSigma=TMath::Sqrt(2*TMath::Abs(fitter.GetParameter(2)));
144 Double_t parMean=-fitter.GetParameter(1)/(2.*fitter.GetParameter(2));
145 histo.Fit(&fg,"QN","QN", paramLTM[1]-5*paramLTM[2], paramLTM[1]+5*paramLTM[2]);
146 meanLTMHistoPar[iltm]=parMean;
147 sigmaLTMHistoPar[iltm]=parSigma;
148 meanLTMHistoGausLim[iltm]=fg.GetParameter(1);
149 sigmaLTMHistoGausLim[iltm]=fg.GetParameter(2);
150 }else{
151 meanLTMHistoPar[iltm]=paramLTM[1];
152 sigmaLTMHistoPar[iltm]=paramLTM[2];
153 meanLTMHistoGausLim[iltm]=meanG;
154 sigmaLTMHistoGausLim[iltm]=rmsG;
155 }
156
8ecd39c2 157 }
158 (*pcstream)<<"ltm"<<
159 //Input
160 "npoints="<<npoints<<
161 "distType="<<distType<<
162 "mean="<<mean<<
163 "sigma="<<sigma<<
164 // "Standart" statistic output
165 "meanA="<<meanA<<
166 "rmsA="<<rmsA<<
167 "meanH="<<meanH<<
168 "rmsH="<<rmsH<<
169 "meanG="<<meanG<<
170 "rmsG="<<rmsG<<
171 // "LTM" output
172 "vecLTM.="<<&vecLTM<<
173 "meanLTM.="<<&meanLTM<<
174 "meanLTMHisto.="<<&meanLTMHisto<<
175 "meanLTMHisto1.="<<&meanLTMHisto1<<
71be42ab 176 "meanLTMHistoPar.="<<&meanLTMHistoPar<<
177 "meanLTMHistoGausLim.="<<&meanLTMHistoGausLim<<
8ecd39c2 178 //
179 "sigmaLTM.="<<&sigmaLTM<<
180 "sigmaLTMHisto.="<<&sigmaLTMHisto<<
71be42ab 181 "sigmaLTMHistoPar.="<<&sigmaLTMHistoPar<<
182 "sigmaLTMHistoGausLim.="<<&sigmaLTMHistoGausLim<<
8ecd39c2 183 "\n";
184 }
185 delete pcstream;
186 //
187 //
188 //
189 TFile *fltm = TFile::Open("TStatToolkit_LTMtest.root","update");
190 TTree * tree = (TTree*)fltm->Get("ltm");
191 tree->SetMarkerSize(0.5);
192 tree->SetMarkerStyle(25);
193 //
194 // 1. Get numerical error estimate of the LTM method for gaussian distribution
195 //
196 TH2 * hisLTMTrunc[10]={0};
197 TH1 * hisResol[10]={0};
198 TGraphErrors * grSigma[10]={0};
199 TCanvas *canvasLTM= new TCanvas("canvasLTM","canvasLTM",800,700);
71be42ab 200 canvasLTM->Divide(2,2,0,0);
8ecd39c2 201 for (Int_t itype=0; itype<4; itype++){
202 canvasLTM->cd(itype+1);
203 TCut cutType=TString::Format("distType==0%d",itype).Data();
204 tree->Draw("sqrt(npoints-2)*(meanLTM.fElements-mean)/sigma:vecLTM.fElements>>hisLTM(6,0.45,1.05,100,-10,10)",cutType+"npoints>50&&sigma>0.1","colzgof");
205 hisLTMTrunc[0]= (TH2*)(tree->GetHistogram()->Clone());
206 tree->Draw("sqrt(npoints-2)*(meanLTMHisto.fElements-mean)/sigma:vecLTM.fElements>>hisLTMHisto(6,0.45,1.05,100,-10,10)",cutType+"npoints>50&&sigma>0.1","colzgoff");
207 hisLTMTrunc[1]= (TH2*)tree->GetHistogram()->Clone();
71be42ab 208 tree->Draw("sqrt(npoints-2)*(meanLTMHistoPar.fElements-mean)/sigma:vecLTM.fElements>>hisLTMHist1(6,0.45,1.05,100,-10,10)",cutType+"npoints>50&&sigma>0.1","colzgoff");
8ecd39c2 209 hisLTMTrunc[2]= (TH2*)tree->GetHistogram()->Clone();
71be42ab 210 tree->Draw("sqrt(npoints-2)*(meanG-mean)/sigma:vecLTM.fElements>>hisLTMHist1(6,0.45,1.05,100,-10,10)",cutType+"npoints>50&&sigma>0.1","colzgoff");
211 hisLTMTrunc[3]= (TH2*)tree->GetHistogram()->Clone();
212 tree->Draw("sqrt(npoints-2)*(meanLTMHistoGausLim.fElements-mean)/sigma:vecLTM.fElements>>hisLTMHist1(6,0.45,1.05,100,-10,10)",cutType+"npoints>50&&sigma>0.1","colzgoff");
213 hisLTMTrunc[4]= (TH2*)tree->GetHistogram()->Clone();
214 TLegend * legend = new TLegend(0.5,0.7,0.88,0.88,distNames[itype]);
8ecd39c2 215 legend->SetBorderSize(0);
71be42ab 216 for (Int_t ihis=0; ihis<5; ihis++){
8ecd39c2 217 // MakeStat1D(TH2 * his, Int_t deltaBin, Double_t fraction, Int_t returnType, Int_t markerStyle, Int_t markerColor);
71be42ab 218 grSigma[ihis]=TStatToolkit::MakeStat1D( hisLTMTrunc[ihis],0,0.99,1,kMarkers[ihis],kColors[ihis]);
219 grSigma[ihis]->SetMaximum(5.0);
220 grSigma[ihis]->SetMinimum(0.5);
8ecd39c2 221 grSigma[ihis]->GetXaxis()->SetTitle("LTM fraction");
71be42ab 222 grSigma[ihis]->GetYaxis()->SetTitle("#sigma_{meas}/#sigma_{Theor.}");
8ecd39c2 223 if (ihis==0)grSigma[ihis]->Draw("alp");
224 if (ihis>0)grSigma[ihis]->Draw("lp");
225 legend->AddEntry(grSigma[ihis],names[ihis],"p");
226 }
227 legend->Draw();
228 }
229 canvasLTM->SaveAs("robustStatLTM_Performance.pdf");
230
231
232}