Fork me on GitHub

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

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

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