]>
Commit | Line | Data |
---|---|---|
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 | ||
18 | TObjArray arrayFit(3); | |
19 | Int_t kMarkers[10]={25,24,20,21}; | |
20 | Int_t kColors[10]={1,2,4,3}; | |
21 | const char * names[10] = {"LTM","LThisto","LTMhisto1"}; | |
22 | const char * distNames[10] = {"Gauss","Gauss+flat","Gauss(m)+0.2*Gauss(m+5#sigma)","Gauss(m)+0.3*Gauss(m+3#sigma)"}; | |
23 | ||
24 | ||
25 | void TestLTM(Int_t nevents=10000){ | |
26 | // | |
27 | // Goal test and benchamerk numerical stability and precission of the LTM method | |
28 | // Binned data and not binned data: | |
29 | // Here we assume that binwidth<sigma | |
30 | // Distribution types examples used in test: | |
31 | // | |
32 | // 0 - gaussian | |
33 | // 1 - gauss+uniform background | |
34 | // 2 - gauss+second gaus 5 sigma away | |
35 | // | |
36 | Int_t npointsMax=1000; | |
37 | TTreeSRedirector * pcstream = new TTreeSRedirector("TStatToolkit_LTMtest.root","recreate"); | |
38 | // | |
39 | TVectorD values(npointsMax); | |
40 | TVectorD vecLTM(6); | |
41 | TVectorD meanLTM(6); | |
42 | TVectorD sigmaLTM(6); | |
43 | TVectorD meanLTMHisto(6); | |
44 | TVectorD sigmaLTMHisto(6); | |
45 | TVectorD meanLTMHisto1(6); | |
46 | TVectorD sigmaLTMHisto1(6); | |
47 | TVectorD paramLTM(10); | |
48 | // | |
49 | for (Int_t iltm=0; iltm<6; iltm++) vecLTM[iltm]=0.99999*(0.5+iltm/10.); | |
50 | TF1 fg("fg","gaus"); | |
51 | for (Int_t ievent = 0; ievent<nevents; ievent++){ | |
52 | if (ievent%1000==0) printf("%d\n",ievent); | |
53 | Int_t distType=Int_t(gRandom->Rndm()*4); | |
54 | Int_t npoints= 50+(npointsMax-50)*gRandom->Rndm(); | |
55 | Double_t mean = 0.5+(gRandom->Rndm()-0.5)*0.2; | |
56 | Double_t sigma = mean*0.2*(1+2*gRandom->Rndm()); | |
57 | TH1F histo("histo","histo",100,-0.5,1.5); | |
58 | // | |
59 | for (Int_t ipoint=0; ipoint<npoints; ipoint++){ | |
60 | Double_t value=0; | |
61 | if (distType==0) value=gRandom->Gaus(mean,sigma); // gauss | |
62 | if (distType==1) { | |
63 | if (gRandom->Rndm()>0.2) { // gauss + background | |
64 | value=gRandom->Gaus(mean,sigma); | |
65 | }else{ | |
66 | value=mean+(gRandom->Rndm()-0.5)*2; | |
67 | } | |
68 | } | |
69 | if (distType==2) { | |
70 | if (gRandom->Rndm()>0.2) { // gauss + second gaus 5 sigma away | |
71 | value=gRandom->Gaus(mean,sigma); | |
72 | }else{ | |
73 | value=gRandom->Gaus(mean+5*sigma,sigma); | |
74 | } | |
75 | } | |
76 | if (distType==3) { | |
77 | if (gRandom->Rndm()>0.3) { // gauss + second gaus 4 sigma away | |
78 | value=gRandom->Gaus(mean,sigma); | |
79 | }else{ | |
80 | value=gRandom->Gaus(mean+5*sigma,sigma); | |
81 | } | |
82 | } | |
83 | values[ipoint]=value; | |
84 | histo.Fill(value); | |
85 | } | |
86 | // | |
87 | histo.Fit(&fg,"QN","QN"); | |
88 | Double_t meanG = fg.GetParameter(1); | |
89 | Double_t rmsG = fg.GetParameter(2); | |
90 | Double_t meanA = TMath::Mean(npoints,values.GetMatrixArray()); | |
91 | Double_t rmsA = TMath::Mean(npoints,values.GetMatrixArray()); | |
92 | Double_t meanH = histo.GetMean(); | |
93 | Double_t rmsH = histo.GetRMS(); | |
94 | // | |
95 | for (Int_t iltm=0; iltm<6; iltm++){ | |
96 | // | |
97 | Double_t meanV,sigmaV=0; | |
98 | TStatToolkit::EvaluateUni(npoints,values.GetMatrixArray(), meanV,sigmaV, vecLTM[iltm]*npoints); | |
99 | meanLTM[iltm]=meanV; | |
100 | sigmaLTM[iltm]=sigmaV; | |
101 | // | |
102 | TStatToolkit::LTM(&histo, ¶mLTM, vecLTM[iltm]); | |
103 | meanLTMHisto1[iltm]=paramLTM[1]; | |
104 | sigmaLTMHisto1[iltm]=paramLTM[2]; | |
105 | // | |
106 | TStatToolkit::LTMHisto(&histo, paramLTM, vecLTM[iltm]); | |
107 | meanLTMHisto[iltm]=paramLTM[1]; | |
108 | sigmaLTMHisto[iltm]=paramLTM[2]; | |
109 | } | |
110 | (*pcstream)<<"ltm"<< | |
111 | //Input | |
112 | "npoints="<<npoints<< | |
113 | "distType="<<distType<< | |
114 | "mean="<<mean<< | |
115 | "sigma="<<sigma<< | |
116 | // "Standart" statistic output | |
117 | "meanA="<<meanA<< | |
118 | "rmsA="<<rmsA<< | |
119 | "meanH="<<meanH<< | |
120 | "rmsH="<<rmsH<< | |
121 | "meanG="<<meanG<< | |
122 | "rmsG="<<rmsG<< | |
123 | // "LTM" output | |
124 | "vecLTM.="<<&vecLTM<< | |
125 | "meanLTM.="<<&meanLTM<< | |
126 | "meanLTMHisto.="<<&meanLTMHisto<< | |
127 | "meanLTMHisto1.="<<&meanLTMHisto1<< | |
128 | // | |
129 | "sigmaLTM.="<<&sigmaLTM<< | |
130 | "sigmaLTMHisto.="<<&sigmaLTMHisto<< | |
131 | "sigmaLTMHisto1.="<<&sigmaLTMHisto1<< | |
132 | "\n"; | |
133 | } | |
134 | delete pcstream; | |
135 | // | |
136 | // | |
137 | // | |
138 | TFile *fltm = TFile::Open("TStatToolkit_LTMtest.root","update"); | |
139 | TTree * tree = (TTree*)fltm->Get("ltm"); | |
140 | tree->SetMarkerSize(0.5); | |
141 | tree->SetMarkerStyle(25); | |
142 | // | |
143 | // 1. Get numerical error estimate of the LTM method for gaussian distribution | |
144 | // | |
145 | TH2 * hisLTMTrunc[10]={0}; | |
146 | TH1 * hisResol[10]={0}; | |
147 | TGraphErrors * grSigma[10]={0}; | |
148 | TCanvas *canvasLTM= new TCanvas("canvasLTM","canvasLTM",800,700); | |
149 | canvasLTM->Divide(2,2); | |
150 | for (Int_t itype=0; itype<4; itype++){ | |
151 | canvasLTM->cd(itype+1); | |
152 | TCut cutType=TString::Format("distType==0%d",itype).Data(); | |
153 | 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"); | |
154 | hisLTMTrunc[0]= (TH2*)(tree->GetHistogram()->Clone()); | |
155 | 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"); | |
156 | hisLTMTrunc[1]= (TH2*)tree->GetHistogram()->Clone(); | |
157 | tree->Draw("sqrt(npoints-2)*(meanLTMHisto1.fElements-mean)/sigma:vecLTM.fElements>>hisLTMHist1(6,0.45,1.05,100,-10,10)",cutType+"npoints>50&&sigma>0.1","colzgoff"); | |
158 | hisLTMTrunc[2]= (TH2*)tree->GetHistogram()->Clone(); | |
159 | TLegend * legend = new TLegend(0.5,0.7,0.9,0.9,distNames[itype]); | |
160 | legend->SetBorderSize(0); | |
161 | for (Int_t ihis=0; ihis<3; ihis++){ | |
162 | // MakeStat1D(TH2 * his, Int_t deltaBin, Double_t fraction, Int_t returnType, Int_t markerStyle, Int_t markerColor); | |
163 | grSigma[ihis]=TStatToolkit::MakeStat1D( hisLTMTrunc[ihis],0,0.99,5,kMarkers[ihis],kColors[ihis]); | |
164 | grSigma[ihis]->GetXaxis()->SetTitle("LTM fraction"); | |
165 | grSigma[ihis]->GetYaxis()->SetTitle("#sigma_{meas}/sigma_{gauss}"); | |
166 | if (ihis==0)grSigma[ihis]->Draw("alp"); | |
167 | if (ihis>0)grSigma[ihis]->Draw("lp"); | |
168 | legend->AddEntry(grSigma[ihis],names[ihis],"p"); | |
169 | } | |
170 | legend->Draw(); | |
171 | } | |
172 | canvasLTM->SaveAs("robustStatLTM_Performance.pdf"); | |
173 | ||
174 | ||
175 | } |