Fork me on GitHub

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

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

update, better&nicer plots

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