Fork me on GitHub

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

Last change on this file since 437 was 264, checked in by Xavier Rouby, 16 years ago

first test 2.0

File size: 9.9 KB
RevLine 
[9]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"
[153]8#include "TLegend.h"
[9]9
10#include "interface/FuncDef.h"
11
[27]12void JetResol()
[9]13{
[39]14 setTDRStyle();
15 gROOT->Reset();
16
[258]17 TFile *f1 = new TFile("JET2.root","read");
[39]18 TTree *Analyze = (TTree*)f1->Get("Analysis");
19
[153]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};
[39]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;
[91]46 sprintf(tempMin,"JetPTResol.PT > %d",binMin);
[39]47 string mystringMin(tempMin);
48 char tempMax[500];
[91]49 sprintf(tempMax,"JetPTResol.PT < %d",binMax);
[39]50 string mystringMax(tempMax);
51 char tempName[500];
[91]52 sprintf(tempName,"(JetPTResol.SmearedPT)>>hETSoverET%d",i);
[39]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
[153]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);
[39]70
[153]71 TGraph *gr11 = new TGraph((numBin-1),x,y);
[39]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();
[153]79
80 fitfun->SetLineColor(6);
[39]81 gr11->Fit("user","QR");
[153]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]);
[27]89
[153]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
[60]99 TPaveText *events = MakeTPave(0.2,0.75,0.35,0.8,"Events: pp #rightarrow gg ");
100 events->Draw();
101
[39]102 TPaveText *Delphes = MakeTPave(0.2,0.15,0.35,0.2,"MG/ME + Delphes");
103 Delphes->Draw();
[153]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
[39]112 delete fitfun;
113}
[27]114
[39]115void ElecResol()
116{
[9]117
[39]118 setTDRStyle();
119 gROOT->Reset();
120
[153]121 TFile *f1 = new TFile("ETMIS2.root","read");
[39]122 TTree *Analyze = (TTree*)f1->Get("Analysis");
123
[153]124 const Int_t numBin=11;
125 double bins[numBin]={0,10,20,30,40,50,60,70,80,100,120};
[39]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;
[89]150 sprintf(tempMin,"ElecEResol.E > %d",binMin);
[39]151 string mystringMin(tempMin);
152 char tempMax[500];
[89]153 sprintf(tempMax,"ElecEResol.E < %d",binMax);
[39]154 string mystringMax(tempMax);
155 char tempName[500];
[89]156 sprintf(tempName,"(ElecEResol.E-ElecEResol.SmearedE)>>hETSoverET%d",i);
[39]157 string mystringName(tempName);
158
159 c3->cd(++frame);
[153]160 GaussValuesElec(Analyze,tempName,rms[i],mean[i],mystringMin,mystringMax);
161 //GaussValues(Analyze,tempName,rms[i],mean[i],mystringMin,mystringMax);
[39]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
[60]196 TPaveText *events = MakeTPave(0.2,0.75,0.35,0.8,"Events: WHq'#rightarrow W#tau#tauq'#rightarrowjjl#tauq', m_{H}=150 GeV ");
[39]197 events->Draw();
198
199 delete fitfun;
200}
[9]201
[39]202void TauJetInfo()
203{
204
205 setTDRStyle();
206 gROOT->Reset();
207
[258]208 TFile *f1 = new TFile("TAUJET2.root","read");
[39]209 TTree *Analyze = (TTree*)f1->Get("Analysis");
210
[60]211 TCanvas *ct1 = new TCanvas("ct1","Tau information",100,100,600,450);
212 ct1->cd();
[39]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();
[153]225 TPaveText *events1 = MakeTPave(0.3,0.65,0.45,0.7,"Events: WHq'#rightarrow WWWq'#rightarrow lllq', m_{H}=150 GeV ");
[39]226 events1->Draw();
227
228
[60]229 TCanvas *ct2 = new TCanvas("ct2","Tau information",100,100,600,450);
230 ct2->cd();
[39]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();
[153]241 TPaveText *events = MakeTPave(0.6,0.65,0.85,0.7,"Events: WHq'#rightarrow WWWq'#rightarrow lllq', m_{H}=150 GeV ");
[39]242 events->Draw();
243
[9]244}
[23]245
[27]246void ETmisResol()
[23]247{
[39]248
249 setTDRStyle();
250 gROOT->Reset();
251
[258]252 TFile *f1 = new TFile("JET2.root","read");
[39]253 TTree *Analyze = (TTree*)f1->Get("Analysis");
254
255 TF1 *fitfun = new TF1("user","[0]*sqrt(x)",0,600);
[153]256 const Int_t numBin=6;
257 double bins[numBin]={180,260,340,420,500,580};
[39]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");
[153]307 gr11->GetYaxis()->SetRangeUser(0,60);
[39]308 gr11->GetXaxis()->SetRangeUser(1,600);
[23]309
[39]310 Double_t* params = fitfun->GetParameters();
311
312 char tempResol[500];
313 sprintf(tempResol,"%f * #sqrt{E_{T}}",params[0]);
[23]314
[39]315 TPaveText *leg1 = MakeTPave(0.4,0.6,0.8,0.65,tempResol);
316 leg1->Draw();
[23]317
[153]318 TPaveText *leg2 = MakeTPave(0.2,0.8,0.8,0.85,"WHq'#rightarrow WWWq'#rightarrow lllq', m_{H}=150 GeV ");
[39]319 leg2->Draw();
[27]320
[39]321 TPaveText *Delphes = MakeTPave(0.2,0.15,0.35,0.2,"MG/ME + Delphes");
322 Delphes->Draw();
323
324
325 delete fitfun;
[23]326}
327
328void General()
329{
[264]330 JetResol();
331 ElecResol();
332 ETmisResol();
333 //TauJetInfo();
[39]334
[27]335}
[23]336
337
Note: See TracBrowser for help on using the repository browser.