Fork me on GitHub

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

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

update des labels / securites dans le code

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