Fork me on GitHub

source: svn/trunk/routines/resolutions_atlas.C@ 471

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

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