Fork me on GitHub

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

Last change on this file since 448 was 443, checked in by Xavier Rouby, 15 years ago

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