Fork me on GitHub

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

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

added CMS expectation band for etmis

File size: 14.7 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 "TMultiGraph.h"
42#include "TGraphErrors.h"
43#include "TLegend.h"
44
45#include "interface/FuncDef.h"
46
47void JetResol()
48{
49 setTDRStyle();
50 gROOT->Reset();
51
52 TFile *f1 = new TFile("JET2.root","read");
53 if(!f1->IsOpen()) { cout << "could not open "<< f1->GetName() << ". Exiting..." << endl; return;}
54 if(!f1->FindKey("Analysis")) { cout << "Bad input file, could not find the \"Analysis\" tree. Exiting..." << endl; return;}
55 TTree *Analyze = (TTree*)f1->Get("Analysis");
56
57 const Int_t numBin=16;
58 double bins[numBin]={0,10,20,30,40,50,60,70,80,100,120,140,180,220,300,1200};
59 TProfile *ETSoverET = new TProfile("ETSoverET","Jet resolution: E_{T}^{rec}/E_{T}^{mc}",(numBin-1),bins,-10,10);
60
61 double rms[numBin];
62 double mean[numBin];
63
64 TCanvas *c1 = new TCanvas("c1","JET resol",0,0,1000,650);
65 c1->cd(); int frame=0;
66 c1->Divide(6,2);
67
68 float x[numBin-1];
69 float y[numBin-1];
70 float ex[numBin-1];
71 float ey[numBin-1];
72 float erec_over_etruth[numBin-1];
73
74 float finval=0;//valeur a remplir puis a fitter
75
76 for ( int i=0; i<(numBin-1); i++) // premiÚre bin : i ==1 et pas i == 0
77 {
78 TAxis *xaxis = ETSoverET->GetXaxis();
79 float binCenter = xaxis->GetBinCenter(i+1);
80 int binMin = (int)xaxis->GetBinLowEdge(i+1);
81 int binMax = (int)xaxis->GetBinUpEdge(i+1);
82 char tempMin[500];
83 if(i==0)binMin=5;
84 sprintf(tempMin,"JetPTResol.PT > %d",binMin);
85 string mystringMin(tempMin);
86 char tempMax[500];
87 sprintf(tempMax,"JetPTResol.PT < %d",binMax);
88 string mystringMax(tempMax);
89 char tempName[500];
90 sprintf(tempName,"(JetPTResol.SmearedPT)>>hETSoverET%d",i);
91 string mystringName(tempName);
92
93 c1->cd(++frame);
94 GaussValues(Analyze,tempName,rms[i],mean[i],mystringMin,mystringMax);
95 x[i]=binCenter;
96 finval=rms[i]/mean[i];
97 y[i]=(finval*100);
98 ex[i]=0;
99 ey[i]=0;
100 erec_over_etruth[i]=mean[i]+1;
101 }
102
103
104 TCanvas *c3 = new TCanvas("c3","JET linearity",100,100,600,450);
105 c3->cd();
106 //for(int i=0; i<numBin-1; i++) cout << x[i] << "\t\t" << erec_over_etruth[i] << endl;
107 TGraph *linearity = new TGraph(numBin-1,x,erec_over_etruth);
108 linearity->Draw("AP");
109 linearity->SetTitle("");
110 linearity->GetXaxis()->SetTitle("E^{Truth} [GeV]");
111 linearity->GetYaxis()->SetTitle("<E^{Reco}/E^{Truth}>");
112 linearity->GetXaxis()->SetRangeUser(30,1000);
113 linearity->SetMinimum(0.85);
114 linearity->SetMaximum(1.20);
115 gPad->SetLogx();
116
117
118 TCanvas *c2 = new TCanvas("c2","JET resol",100,100,600,450);
119 c2->cd();
120
121 TF1 *fitfun = new TF1("user","sqrt(pow([0]/x,2)+pow([1]/sqrt(x),2)+pow([2],2))",10,800);
122 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);
123 //TF1 *fitfunCMS = new TF1("userCMS","sqrt(pow(1.17*100/sqrt(x),2)+pow(0.039*100,2))",10,800);
124
125 TGraph *gr11 = new TGraph((numBin-1),x,y);
126 gr11->Draw("AP");
127 gr11->SetTitle("");
128 gr11->GetXaxis()->SetTitle("E_{T}^{MC} [GeV]");
129 gr11->GetYaxis()->SetRangeUser(0,50);
130 gr11->GetYaxis()->SetTitle("#sigma(E_{T}^{rec}/E_{T}^{MC})_{fit}/<E_{T}^{rec}/E_{T}^{MC}>_{fit}");
131
132 Double_t* params = fitfun->GetParameters();
133
134 fitfun->SetLineColor(6);
135 gr11->Fit("user","QR");
136 gr11->Fit("user","QR");
137 gr11->Fit("user","QR");
138 gr11->Fit("user","QR");
139 char tempResol[100];
140 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]);
141 char tempResol2[100];
142 sprintf(tempResol2,"sqrt(pow(%f/sqrt(x),2)+pow(%f,2))",params[1],params[2]);
143
144
145 TF1 *fitfunDelphes = new TF1("userDelphes",tempResol2,7,1000);
146 fitfunDelphes->SetLineColor(6);
147 fitfunDelphes->Draw("same");
148
149 //CMS values
150 fitfunCMS->SetLineColor(9);
151 fitfunCMS->Draw("same");
152
153 TPaveText *events = MakeTPave(0.2,0.75,0.35,0.8,"Events: pp #rightarrow gg ");
154 events->Draw();
155
156 TPaveText *Delphes = MakeTPave(0.2,0.15,0.35,0.2,"MG/ME + Delphes");
157 Delphes->Draw();
158
159 TLegend *legend = new TLegend(0.2,0.6,0.9,0.85,NULL,"NDC");
160 legend->AddEntry(fitfunCMS,"CMS resolution","l");
161 legend->AddEntry(fitfunDelphes,tempResol,"l");
162 legend->SetFillColor(10);
163 legend->SetBorderSize(0);
164 legend->Draw();
165
166 delete fitfun;
167}
168
169void ElecResol()
170{
171
172 setTDRStyle();
173 gROOT->Reset();
174
175 TFile *f1 = new TFile("ETMIS2.root","read");
176 if(!f1->IsOpen()) { cout << "could not open "<< f1->GetName() << ". Exiting..." << endl; return;}
177 if(!f1->FindKey("Analysis")) { cout << "Bad input file, could not find the \"Analysis\" tree. Exiting..." << endl; return;}
178 TTree *Analyze = (TTree*)f1->Get("Analysis");
179
180 const Int_t numBin=11;
181 double bins[numBin]={0,10,20,30,40,50,60,70,80,100,120};
182 TProfile *ETSminusET = new TProfile("ETSminusET","Electron resolution: E_{T}^{rec}/E_{T}^{mc}",(numBin-1),bins,-10,10);
183
184 double rms[numBin];
185 double mean[numBin];
186
187 TCanvas *c3 = new TCanvas("c3","ELEC resol",0,0,1000,650);
188 c3->cd(); int frame=0;
189 c3->Divide(6,2);
190
191 float x[numBin-1];
192 float y[numBin-1];
193 float ex[numBin-1];
194 float ey[numBin-1];
195
196 float finval=0;//valeur a remplir puis a fitter
197
198 for ( int i=0; i<(numBin-1); i++) // premiÚre bin : i ==1 et pas i == 0
199 {
200 TAxis *xaxis = ETSminusET->GetXaxis();
201 float binCenter = xaxis->GetBinCenter(i+1);
202 int binMin = (int)xaxis->GetBinLowEdge(i+1);
203 int binMax = (int)xaxis->GetBinUpEdge(i+1);
204 char tempMin[500];
205 if(i==0)binMin=5;
206 sprintf(tempMin,"ElecEResol.E > %d",binMin);
207 string mystringMin(tempMin);
208 char tempMax[500];
209 sprintf(tempMax,"ElecEResol.E < %d",binMax);
210 string mystringMax(tempMax);
211 char tempName[500];
212 sprintf(tempName,"(ElecEResol.E-ElecEResol.SmearedE)>>hETSoverET%d",i);
213 string mystringName(tempName);
214
215 c3->cd(++frame);
216 GaussValuesElec(Analyze,tempName,rms[i],mean[i],mystringMin,mystringMax);
217 //GaussValues(Analyze,tempName,rms[i],mean[i],mystringMin,mystringMax);
218 x[i]=binCenter;
219 finval=rms[i]/binCenter;
220 y[i]=(finval*100);
221 ex[i]=0;
222 ey[i]=0;
223 }
224
225 TCanvas *c4 = new TCanvas("c4","ELEC resol",100,100,600,450);
226 c4->cd();
227
228 TF1 *fitfun = new TF1("user","sqrt(pow([0],2)+pow([1]/sqrt(x),2)+pow([2]/x,2))",1,400);
229
230 TGraphErrors *gr11 = new TGraphErrors((numBin-1),x,y,ex,ey);
231 gr11->Draw("AP");
232 gr11->SetTitle("");
233 gr11->GetXaxis()->SetTitle("E [GeV]");
234 gr11->GetYaxis()->SetRangeUser(0,5);
235 gr11->GetYaxis()->SetTitle("#sigma/E");
236
237 Double_t* params = fitfun->GetParameters();
238
239 gr11->Fit("user","QR");
240 gr11->Fit("user","QRI");
241 gr11->Fit("user","QRI");
242 gr11->Fit("user","QRI");
243 char tempResol[500];
244 sprintf(tempResol,"#frac{#sigma}{E} = #frac{%f}{#sqrt{E}} #oplus #frac{%f}{E} #oplus %f",params[1]/100,params[2]/100,params[0]/100);
245
246 TPaveText *leg1 = MakeTPave(0.4,0.6,0.8,0.65,tempResol);
247 leg1->Draw();
248
249 TPaveText *Delphes = MakeTPave(0.2,0.15,0.35,0.2,"MG/ME + Delphes");
250 Delphes->Draw();
251
252 TPaveText *events = MakeTPave(0.2,0.75,0.35,0.8,"Events: WHq'#rightarrow W#tau#tauq'#rightarrowjjl#tauq', m_{H}=150 GeV ");
253 events->Draw();
254
255 delete fitfun;
256}
257
258void TauJetInfo()
259{
260
261 setTDRStyle();
262 gROOT->Reset();
263
264 TFile *f1 = new TFile("TAUJET2.root","read");
265 if(!f1->IsOpen()) { cout << "could not open "<< f1->GetName() << ". Exiting..." << endl; return;}
266 if(!f1->FindKey("Analysis")) { cout << "Bad input file, could not find the \"Analysis\" tree. Exiting..." << endl; return;}
267 TTree *Analyze = (TTree*)f1->Get("Analysis");
268
269 TCanvas *ct1 = new TCanvas("ct1","Tau information",100,100,600,450);
270 ct1->cd();
271
272 TH1F *tauEnergy =MakeNormTH1F(20,0.8,1,Analyze,"TauJetPTResol.EnergieCen>>tauEnergy",1, 0, 1,2,false);
273 tauEnergy->Draw();
274 tauEnergy->SetTitle("");
275 tauEnergy->GetYaxis()->SetTitle("Fraction of events");
276 tauEnergy->GetXaxis()->SetTitle("C_{#tau}");
277
278
279 TPaveText *Delphes1 = MakeTPave(0.3,0.85,0.45,0.9,"MG/ME + Delphes");
280 Delphes1->Draw();
281 TPaveText *ctau = MakeTPave(0.3,0.75,0.45,0.8,"C_{#tau} = #frac{#sum E^{towers} (#DeltaR = 0.15)}{E^{jet}}");
282 ctau->Draw();
283 TPaveText *events1 = MakeTPave(0.3,0.65,0.45,0.7,"Events: WHq'#rightarrow WWWq'#rightarrow lllq', m_{H}=150 GeV ");
284 events1->Draw();
285
286
287 TCanvas *ct2 = new TCanvas("ct2","Tau information",100,100,600,450);
288 ct2->cd();
289 TH1F *NumTrack =MakeNormTH1F(6,0,6,Analyze,"TauJetPTResol.NumTrack>>NumTrack",1, 0, 1,2,false);
290 NumTrack->Draw();
291 NumTrack->SetTitle("");
292 NumTrack->GetYaxis()->SetTitle("Fraction of events");
293 NumTrack->GetXaxis()->SetTitle("N^{tracks}");
294
295 TPaveText *Delphes = MakeTPave(0.6,0.85,0.85,0.9,"MG/ME + Delphes");
296 Delphes->Draw();
297 TPaveText *numtracks = MakeTPave(0.6,0.75,0.85,0.8,"#DeltaR < 0.4, p_{T}^{track} > 2 GeV");
298 numtracks->Draw();
299 TPaveText *events = MakeTPave(0.6,0.65,0.85,0.7,"Events: WHq'#rightarrow WWWq'#rightarrow lllq', m_{H}=150 GeV ");
300 events->Draw();
301
302}
303
304void ETmisResol()
305{
306
307 setTDRStyle();
308 gROOT->Reset();
309
310 //TFile *f1 = new TFile("JET2.root","read");
311 TFile *f1 = new TFile("JET2_cms_midpoint.root","read");
312
313 if(!f1->IsOpen()) { cout << "could not open "<< f1->GetName() << ". Exiting..." << endl; return;}
314 if(!f1->FindKey("Analysis")) { cout << "Bad input file, could not find the \"Analysis\" tree. Exiting..." << endl; return;}
315 TTree *Analyze = (TTree*)f1->Get("Analysis");
316
317 TF1 *fitfun = new TF1("user","[0]*sqrt(x)",0,600);
318 fitfun->SetLineColor(kBlack);
319 fitfun->SetLineStyle(2);
320 const Int_t numBin=6;
321 double bins[numBin]={180,260,340,420,500,580};
322 TProfile *ETSoverET = new TProfile("ETSoverET","ETmis resol",(numBin-1),bins,-1000,1000);
323
324
325 double rms[numBin-1];
326
327 TCanvas *c5 = new TCanvas("c5","PTmis resol",0,0,1000,650);
328 c5->cd(); int frame=0;
329 c5->Divide(6,2);
330
331 double x[numBin];
332 double y[numBin];
333
334 for ( int i=0; i<(numBin-1); i++) // premiÚre bin : i ==1 et pas i == 0
335 {
336 TAxis *xaxis = ETSoverET->GetXaxis();
337 float binCenter = xaxis->GetBinCenter(i+1);
338 int binMin = (int)xaxis->GetBinLowEdge(i+1);
339 int binMax = (int)xaxis->GetBinUpEdge(i+1);
340 char tempMin[500];
341 sprintf(tempMin,"ETmisResol.SEt>%d",binMin);
342 string mystringMin(tempMin);
343 char tempMax[500];
344 sprintf(tempMax,"ETmisResol.SEt<%d",binMax);
345 string mystringMax(tempMax);
346
347 char tempName[500];
348 sprintf(tempName,"ETmisResol.ExSmeare>>hETSoverET%d",i);
349 string mystringName(tempName);
350
351 c5->cd(++frame);
352 GaussValuesETmis(Analyze,tempName,rms[i],mystringMin,mystringMax);
353 x[i]=binCenter;
354 y[i]=rms[i];
355
356 }
357
358 TCanvas *c6 = new TCanvas("c6","ETmis resolution",100,100,600,450);
359 c6->cd();
360
361 TF1 *fdo = new TF1("fdo","0.6*sqrt(x)",1,600);
362 TF1 *fup = new TF1("fup","0.7*sqrt(x)",1,600);
363 int N = 100;
364 float xx[N], yy[N], eex[N], eey[N];
365 for (int i=0; i<N; i++) {
366 xx[i] = 6*i;
367 yy[i] = (fup->Eval(xx[i]) + fdo->Eval(xx[i]))/2.;
368 eex[i] = 0;
369 eey[i] = (fup->Eval(xx[i]) - fdo->Eval(xx[i]))/2.;
370 }
371 TGraphErrors *gr = new TGraphErrors(N,xx,yy,eex,eey);
372 gr->SetFillColor(18);
373
374
375 x[numBin]=0;
376 y[numBin]=0;
377 TGraph *gr11 = new TGraph(numBin,x,rms);
378 gr11->SetTitle("");
379 gr11->SetMarkerStyle(3);
380 gr11->SetLineStyle(2);
381 gr11->Fit("user","RQ");
382 gr11->Fit("user","RQ");
383 gr11->Fit("user","RQ");
384 gr11->Fit("user","RQ");
385
386
387 TMultiGraph *mgr = new TMultiGraph("mgr","");
388 mgr->Add(gr,"3");
389 mgr->Add(gr11,"P");
390 mgr->Draw("A");
391 mgr->GetXaxis()->SetTitle("#Sigma E_{T} [GeV]");
392 mgr->GetYaxis()->SetTitle("Resolution of x-component of MET [GeV]");
393 mgr->GetYaxis()->SetRangeUser(0,60);
394 mgr->GetXaxis()->SetRangeUser(1,600);
395 mgr->SetMaximum(16);
396
397 Double_t* params = fitfun->GetParameters();
398
399 char tempResol[500];
400 sprintf(tempResol,"E_{x}^{mis} resolution: %.2f #times #sqrt{E_{T}}",params[0]);
401
402 TPaveText *leg1 = MakeTPave(0.26,0.84,0.43,0.89,tempResol);
403 leg1->Draw();
404
405 TPaveText *leg2 = MakeTPave(0.15,0.75,0.45,0.8,"Events: pp #rightarrow ggX");
406 leg2->Draw();
407
408 TPaveText *Delphes = MakeTPave(0.65,0.19,0.8,0.24,"MG/ME + Pythia + Delphes");
409 Delphes->Draw();
410
411 TPaveText *cms = MakeTPave(0.65,0.25,0.8,0.30,"CMS-like detector");
412 cms->Draw();
413
414 TLegend *legend = new TLegend(0.69,0.6,0.94,0.74,NULL,"NDC");
415 legend->AddEntry(gr,"CMS","f");
416 legend->AddEntry(gr11,"Delphes","l");
417 legend->SetFillColor(10);
418 legend->SetBorderSize(0);
419 legend->Draw();
420
421
422 delete fitfun;
423}
424
425void General()
426{
427 JetResol();
428 ElecResol();
429 ETmisResol();
430 //TauJetInfo();
431
432}
433
434
Note: See TracBrowser for help on using the repository browser.