Fork me on GitHub

source: svn/trunk/routines/resolutions.C@ 258

Last change on this file since 258 was 258, checked in by severine ovyn, 15 years ago

resolution from delphes output

File size: 9.9 KB
Line 
1#include "TROOT.h"
2#include "TFile.h"
3#include "TTree.h"
4#include "TCanvas.h"
5#include "TProfile.h"
6#include "TF1.h"
7#include "TGraph.h"
8#include "TLegend.h"
9
10#include "interface/FuncDef.h"
11
12void JetResol()
13{
14 setTDRStyle();
15 gROOT->Reset();
16
17 TFile *f1 = new TFile("JET2.root","read");
18 TTree *Analyze = (TTree*)f1->Get("Analysis");
19
20 const Int_t numBin=16;
21 double bins[numBin]={0,10,20,30,40,50,60,70,80,100,120,140,180,220,300,1200};
22 TProfile *ETSoverET = new TProfile("ETSoverET","Jet resolution: E_{T}^{rec}/E_{T}^{mc}",(numBin-1),bins,-10,10);
23
24 double rms[numBin];
25 double mean[numBin];
26
27 TCanvas *c1 = new TCanvas("c1","JET resol",0,0,1000,650);
28 c1->cd(); int frame=0;
29 c1->Divide(6,2);
30
31 float x[numBin-1];
32 float y[numBin-1];
33 float ex[numBin-1];
34 float ey[numBin-1];
35
36 float finval=0;//valeur a remplir puis a fitter
37
38 for ( int i=0; i<(numBin-1); i++) // premiÚre bin : i ==1 et pas i == 0
39 {
40 TAxis *xaxis = ETSoverET->GetXaxis();
41 float binCenter = xaxis->GetBinCenter(i+1);
42 int binMin = (int)xaxis->GetBinLowEdge(i+1);
43 int binMax = (int)xaxis->GetBinUpEdge(i+1);
44 char tempMin[500];
45 if(i==0)binMin=5;
46 sprintf(tempMin,"JetPTResol.PT > %d",binMin);
47 string mystringMin(tempMin);
48 char tempMax[500];
49 sprintf(tempMax,"JetPTResol.PT < %d",binMax);
50 string mystringMax(tempMax);
51 char tempName[500];
52 sprintf(tempName,"(JetPTResol.SmearedPT)>>hETSoverET%d",i);
53 string mystringName(tempName);
54
55 c1->cd(++frame);
56 GaussValues(Analyze,tempName,rms[i],mean[i],mystringMin,mystringMax);
57 x[i]=binCenter;
58 finval=rms[i]/mean[i];
59 y[i]=(finval*100);
60 ex[i]=0;
61 ey[i]=0;
62 }
63
64 TCanvas *c2 = new TCanvas("c2","JET resol",100,100,600,450);
65 c2->cd();
66
67 TF1 *fitfun = new TF1("user","sqrt(pow([0]/x,2)+pow([1]/sqrt(x),2)+pow([2],2))",10,800);
68 TF1 *fitfunCMS = new TF1("userCMS","sqrt(pow(1.05*100/sqrt(x),2)+pow(0.06*100,2)+pow(0*100/x,2))",10,800);
69 //TF1 *fitfunCMS = new TF1("userCMS","sqrt(pow(1.17*100/sqrt(x),2)+pow(0.039*100,2))",10,800);
70
71 TGraph *gr11 = new TGraph((numBin-1),x,y);
72 gr11->Draw("AP");
73 gr11->SetTitle("");
74 gr11->GetXaxis()->SetTitle("E_{T}^{MC} [GeV]");
75 gr11->GetYaxis()->SetRangeUser(0,50);
76 gr11->GetYaxis()->SetTitle("#sigma(E_{T}^{rec}/E_{T}^{MC})_{fit}/<E_{T}^{rec}/E_{T}^{MC}>_{fit}");
77
78 Double_t* params = fitfun->GetParameters();
79
80 fitfun->SetLineColor(6);
81 gr11->Fit("user","QR");
82 gr11->Fit("user","QR");
83 gr11->Fit("user","QR");
84 gr11->Fit("user","QR");
85 char tempResol[100];
86 sprintf(tempResol,"Delphes resolution: #frac{#sigma(E_{T}^{rec}/E_{T}^{MC})}{<E_{T}^{rec}/E_{T}^{MC}>} = #frac{%f}{#sqrt{E_{T}^{MC}}} #oplus %f",params[1],params[2]);
87 char tempResol2[100];
88 sprintf(tempResol2,"sqrt(pow(%f/sqrt(x),2)+pow(%f,2))",params[1],params[2]);
89
90
91 TF1 *fitfunDelphes = new TF1("userDelphes",tempResol2,7,1000);
92 fitfunDelphes->SetLineColor(6);
93 fitfunDelphes->Draw("same");
94
95 //CMS values
96 fitfunCMS->SetLineColor(9);
97 fitfunCMS->Draw("same");
98
99 TPaveText *events = MakeTPave(0.2,0.75,0.35,0.8,"Events: pp #rightarrow gg ");
100 events->Draw();
101
102 TPaveText *Delphes = MakeTPave(0.2,0.15,0.35,0.2,"MG/ME + Delphes");
103 Delphes->Draw();
104
105 TLegend *legend = new TLegend(0.2,0.6,0.9,0.85,NULL,"NDC");
106 legend->AddEntry(fitfunCMS,"CMS resolution","l");
107 legend->AddEntry(fitfunDelphes,tempResol,"l");
108 legend->SetFillColor(10);
109 legend->SetBorderSize(0);
110 legend->Draw();
111
112 delete fitfun;
113}
114
115void ElecResol()
116{
117
118 setTDRStyle();
119 gROOT->Reset();
120
121 TFile *f1 = new TFile("ETMIS2.root","read");
122 TTree *Analyze = (TTree*)f1->Get("Analysis");
123
124 const Int_t numBin=11;
125 double bins[numBin]={0,10,20,30,40,50,60,70,80,100,120};
126 TProfile *ETSminusET = new TProfile("ETSminusET","Electron resolution: E_{T}^{rec}/E_{T}^{mc}",(numBin-1),bins,-10,10);
127
128 double rms[numBin];
129 double mean[numBin];
130
131 TCanvas *c3 = new TCanvas("c3","ELEC resol",0,0,1000,650);
132 c3->cd(); int frame=0;
133 c3->Divide(6,2);
134
135 float x[numBin-1];
136 float y[numBin-1];
137 float ex[numBin-1];
138 float ey[numBin-1];
139
140 float finval=0;//valeur a remplir puis a fitter
141
142 for ( int i=0; i<(numBin-1); i++) // premiÚre bin : i ==1 et pas i == 0
143 {
144 TAxis *xaxis = ETSminusET->GetXaxis();
145 float binCenter = xaxis->GetBinCenter(i+1);
146 int binMin = (int)xaxis->GetBinLowEdge(i+1);
147 int binMax = (int)xaxis->GetBinUpEdge(i+1);
148 char tempMin[500];
149 if(i==0)binMin=5;
150 sprintf(tempMin,"ElecEResol.E > %d",binMin);
151 string mystringMin(tempMin);
152 char tempMax[500];
153 sprintf(tempMax,"ElecEResol.E < %d",binMax);
154 string mystringMax(tempMax);
155 char tempName[500];
156 sprintf(tempName,"(ElecEResol.E-ElecEResol.SmearedE)>>hETSoverET%d",i);
157 string mystringName(tempName);
158
159 c3->cd(++frame);
160 GaussValuesElec(Analyze,tempName,rms[i],mean[i],mystringMin,mystringMax);
161 //GaussValues(Analyze,tempName,rms[i],mean[i],mystringMin,mystringMax);
162 x[i]=binCenter;
163 finval=rms[i]/binCenter;
164 y[i]=(finval*100);
165 ex[i]=0;
166 ey[i]=0;
167 }
168
169 TCanvas *c4 = new TCanvas("c4","ELEC resol",100,100,600,450);
170 c4->cd();
171
172 TF1 *fitfun = new TF1("user","sqrt(pow([0],2)+pow([1]/sqrt(x),2)+pow([2]/x,2))",1,400);
173
174 TGraphErrors *gr11 = new TGraphErrors((numBin-1),x,y,ex,ey);
175 gr11->Draw("AP");
176 gr11->SetTitle("");
177 gr11->GetXaxis()->SetTitle("E [GeV]");
178 gr11->GetYaxis()->SetRangeUser(0,5);
179 gr11->GetYaxis()->SetTitle("#sigma/E");
180
181 Double_t* params = fitfun->GetParameters();
182
183 gr11->Fit("user","QR");
184 gr11->Fit("user","QRI");
185 gr11->Fit("user","QRI");
186 gr11->Fit("user","QRI");
187 char tempResol[500];
188 sprintf(tempResol,"#frac{#sigma}{E} = #frac{%f}{#sqrt{E}} #oplus #frac{%f}{E} #oplus %f",params[1]/100,params[2]/100,params[0]/100);
189
190 TPaveText *leg1 = MakeTPave(0.4,0.6,0.8,0.65,tempResol);
191 leg1->Draw();
192
193 TPaveText *Delphes = MakeTPave(0.2,0.15,0.35,0.2,"MG/ME + Delphes");
194 Delphes->Draw();
195
196 TPaveText *events = MakeTPave(0.2,0.75,0.35,0.8,"Events: WHq'#rightarrow W#tau#tauq'#rightarrowjjl#tauq', m_{H}=150 GeV ");
197 events->Draw();
198
199 delete fitfun;
200}
201
202void TauJetInfo()
203{
204
205 setTDRStyle();
206 gROOT->Reset();
207
208 TFile *f1 = new TFile("TAUJET2.root","read");
209 TTree *Analyze = (TTree*)f1->Get("Analysis");
210
211 TCanvas *ct1 = new TCanvas("ct1","Tau information",100,100,600,450);
212 ct1->cd();
213
214 TH1F *tauEnergy =MakeNormTH1F(20,0.8,1,Analyze,"TauJetPTResol.EnergieCen>>tauEnergy",1, 0, 1,2,false);
215 tauEnergy->Draw();
216 tauEnergy->SetTitle("");
217 tauEnergy->GetYaxis()->SetTitle("Fraction of events");
218 tauEnergy->GetXaxis()->SetTitle("C_{#tau}");
219
220
221 TPaveText *Delphes1 = MakeTPave(0.3,0.85,0.45,0.9,"MG/ME + Delphes");
222 Delphes1->Draw();
223 TPaveText *ctau = MakeTPave(0.3,0.75,0.45,0.8,"C_{#tau} = #frac{#sum E^{towers} (#DeltaR = 0.15)}{E^{jet}}");
224 ctau->Draw();
225 TPaveText *events1 = MakeTPave(0.3,0.65,0.45,0.7,"Events: WHq'#rightarrow WWWq'#rightarrow lllq', m_{H}=150 GeV ");
226 events1->Draw();
227
228
229 TCanvas *ct2 = new TCanvas("ct2","Tau information",100,100,600,450);
230 ct2->cd();
231 TH1F *NumTrack =MakeNormTH1F(6,0,6,Analyze,"TauJetPTResol.NumTrack>>NumTrack",1, 0, 1,2,false);
232 NumTrack->Draw();
233 NumTrack->SetTitle("");
234 NumTrack->GetYaxis()->SetTitle("Fraction of events");
235 NumTrack->GetXaxis()->SetTitle("N^{tracks}");
236
237 TPaveText *Delphes = MakeTPave(0.6,0.85,0.85,0.9,"MG/ME + Delphes");
238 Delphes->Draw();
239 TPaveText *numtracks = MakeTPave(0.6,0.75,0.85,0.8,"#DeltaR < 0.4, p_{T}^{track} > 2 GeV");
240 numtracks->Draw();
241 TPaveText *events = MakeTPave(0.6,0.65,0.85,0.7,"Events: WHq'#rightarrow WWWq'#rightarrow lllq', m_{H}=150 GeV ");
242 events->Draw();
243
244}
245
246void ETmisResol()
247{
248
249 setTDRStyle();
250 gROOT->Reset();
251
252 TFile *f1 = new TFile("JET2.root","read");
253 TTree *Analyze = (TTree*)f1->Get("Analysis");
254
255 TF1 *fitfun = new TF1("user","[0]*sqrt(x)",0,600);
256 const Int_t numBin=6;
257 double bins[numBin]={180,260,340,420,500,580};
258 TProfile *ETSoverET = new TProfile("ETSoverET","ETmis resol",(numBin-1),bins,-1000,1000);
259
260
261 double rms[numBin-1];
262
263 TCanvas *c5 = new TCanvas("c5","PTmis resol",0,0,1000,650);
264 c5->cd(); int frame=0;
265 c5->Divide(6,2);
266
267 double x[numBin];
268 double y[numBin];
269
270 for ( int i=0; i<(numBin-1); i++) // premiÚre bin : i ==1 et pas i == 0
271 {
272 TAxis *xaxis = ETSoverET->GetXaxis();
273 float binCenter = xaxis->GetBinCenter(i+1);
274 int binMin = (int)xaxis->GetBinLowEdge(i+1);
275 int binMax = (int)xaxis->GetBinUpEdge(i+1);
276 char tempMin[500];
277 sprintf(tempMin,"ETmisResol.SEt>%d",binMin);
278 string mystringMin(tempMin);
279 char tempMax[500];
280 sprintf(tempMax,"ETmisResol.SEt<%d",binMax);
281 string mystringMax(tempMax);
282
283 char tempName[500];
284 sprintf(tempName,"ETmisResol.ExSmeare>>hETSoverET%d",i);
285 string mystringName(tempName);
286
287 c5->cd(++frame);
288 GaussValuesETmis(Analyze,tempName,rms[i],mystringMin,mystringMax);
289 x[i]=binCenter;
290 y[i]=rms[i];
291
292 }
293
294 TCanvas *c6 = new TCanvas("c6","ETmis resolution",100,100,600,450);
295 c6->cd();
296
297 x[numBin]=0;
298 y[numBin]=0;
299 TGraph *gr11 = new TGraph((numBin),x,rms);
300 gr11->Draw("AP");
301 gr11->GetXaxis()->SetTitle("Offline sum of E_{T} [GeV]");
302 gr11->GetYaxis()->SetTitle("Resolution of x-component of MET [GeV]");
303 gr11->Fit("user","RQ");
304 gr11->Fit("user","RQ");
305 gr11->Fit("user","RQ");
306 gr11->Fit("user","RQ");
307 gr11->GetYaxis()->SetRangeUser(0,60);
308 gr11->GetXaxis()->SetRangeUser(1,600);
309
310 Double_t* params = fitfun->GetParameters();
311
312 char tempResol[500];
313 sprintf(tempResol,"%f * #sqrt{E_{T}}",params[0]);
314
315 TPaveText *leg1 = MakeTPave(0.4,0.6,0.8,0.65,tempResol);
316 leg1->Draw();
317
318 TPaveText *leg2 = MakeTPave(0.2,0.8,0.8,0.85,"WHq'#rightarrow WWWq'#rightarrow lllq', m_{H}=150 GeV ");
319 leg2->Draw();
320
321 TPaveText *Delphes = MakeTPave(0.2,0.15,0.35,0.2,"MG/ME + Delphes");
322 Delphes->Draw();
323
324
325 delete fitfun;
326}
327
328void General()
329{
330 //JetResol();
331 //ElecResol();
332 //ETmisResol();
333 TauJetInfo();
334
335}
336
337
Note: See TracBrowser for help on using the repository browser.