Fork me on GitHub

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

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

update des labels / securites dans le code

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