ATO-65 - replace FitSlicesY with robust fit function (commit after one run test ...
[u/mrichter/AliRoot.git] / STAT / Macros / TStatToolkitTest.C
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"
17 #include "TLinearFitter.h"
18
19 TObjArray arrayFit(3);
20 Int_t kMarkers[10]={25,24,20,21,22};
21 Int_t kColors[10]={1,2,4,3,6};
22 const char * names[10] = {"LTM","LThisto","LTMhistoPar","Gaus fit", "Gaus in range"};
23 const char * distNames[10] = {"Gauss","Gauss+flat","Gauss(m)+0.2*Gauss(m+5#sigma)","Gauss(m)+0.4*Gauss(m+5#sigma)"};
24
25
26 void 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);
48   TVectorD meanLTMHistoPar(6);
49   TVectorD sigmaLTMHistoPar(6);
50   TVectorD meanLTMHistoGausLim(6);
51   TVectorD sigmaLTMHistoGausLim(6);
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) {
82         if (gRandom->Rndm()>0.4) {                     //  gauss + second gaus 5 sigma away
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];
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
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<<
176       "meanLTMHistoPar.="<<&meanLTMHistoPar<<
177       "meanLTMHistoGausLim.="<<&meanLTMHistoGausLim<<
178       //
179       "sigmaLTM.="<<&sigmaLTM<<
180       "sigmaLTMHisto.="<<&sigmaLTMHisto<<
181       "sigmaLTMHistoPar.="<<&sigmaLTMHistoPar<<      
182       "sigmaLTMHistoGausLim.="<<&sigmaLTMHistoGausLim<<      
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);
200   canvasLTM->Divide(2,2,0,0);
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();
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");
209     hisLTMTrunc[2]= (TH2*)tree->GetHistogram()->Clone();
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]);
215     legend->SetBorderSize(0);
216     for (Int_t ihis=0; ihis<5; ihis++){
217       // MakeStat1D(TH2 * his, Int_t deltaBin, Double_t fraction, Int_t returnType, Int_t markerStyle, Int_t markerColor);
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);
221       grSigma[ihis]->GetXaxis()->SetTitle("LTM fraction");
222       grSigma[ihis]->GetYaxis()->SetTitle("#sigma_{meas}/#sigma_{Theor.}");
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 }