Fork me on GitHub

source: git/examples/Validation.C@ b25184c

ImprovedOutputFile Timing dual_readout llp
Last change on this file since b25184c was 7bb13cd, checked in by Alexandre Mertens <alexandre.mertens@…>, 8 years ago

validation work
  • Property mode set to 100644
File size: 37.0 KB
Line 
1#ifdef __CLING__
2R__LOAD_LIBRARY(libDelphes)
3#include "classes/DelphesClasses.h"
4#include "external/ExRootAnalysis/ExRootTreeReader.h"
5#include "external/ExRootAnalysis/ExRootResult.h"
6#else
7class ExRootTreeReader;
8class ExRootResult;
9#endif
10
11#include "TCanvas.h"
12#include "TSystem.h"
13#include "TStyle.h"
14#include "TLegend.h"
15#include <TH1.h>
16#include "TString.h"
17#include "vector"
18#include <TMath.h>
19#include <iostream>
20#include "TGraph.h"
21#include "TGraphErrors.h"
22#include "TMultiGraph.h"
23#include <typeinfo>
24#include "TLorentzVector.h"
25
26//------------------------------------------------------------------------------
27
28double ptrangemin = 10;
29double ptrangemax = 10000;
30static const int Nbins = 20;
31
32int objStyle = 1;
33int trackStyle = 7;
34int towerStyle = 3;
35
36Color_t objColor = kBlack;
37Color_t trackColor = kBlack;
38Color_t towerColor = kBlack;
39
40double effLegXmin = 0.22;
41double effLegXmax = 0.7;
42double effLegYmin = 0.22;
43double effLegYmax = 0.5;
44
45double resLegXmin = 0.62;
46double resLegXmax = 0.9;
47double resLegYmin = 0.52;
48double resLegYmax = 0.85;
49
50double topLeftLegXmin = 0.22;
51double topLeftLegXmax = 0.7;
52double topLeftLegYmin = 0.52;
53double topLeftLegYmax = 0.85;
54
55
56struct resolPlot
57{
58 TH1 *cenResolHist;
59 TH1 *fwdResolHist;
60 double ptmin;
61 double ptmax;
62 double xmin;
63 double xmax;
64 TString obj;
65
66 resolPlot();
67 resolPlot(double ptdown, double ptup, TString object);
68 void set(double ptdown, double ptup, TString object, double xmin = 0, double xmax = 2);
69 void print(){std::cout << ptmin << std::endl;}
70};
71
72
73resolPlot::resolPlot()
74{
75}
76
77resolPlot::resolPlot(double ptdown, double ptup, TString object)
78{
79 this->set(ptdown,ptup,object);
80}
81
82void resolPlot::set(double ptdown, double ptup, TString object, double xmin, double xmax){
83 ptmin = ptdown;
84 ptmax = ptup;
85 obj = object;
86
87 cenResolHist = new TH1D(obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_cen", obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_cen", 200, xmin, xmax);
88 fwdResolHist = new TH1D(obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_fwd", obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_fwd", 200, xmin, xmax);
89
90}
91
92void HistogramsCollection(std::vector<resolPlot> *histos, double ptmin, double ptmax, TString obj, double xmin = 0, double xmax = 2)
93{
94 double width;
95 double ptdown;
96 double ptup;
97 resolPlot ptemp;
98
99 for (int i = 0; i < Nbins; i++)
100 {
101 width = (ptmax - ptmin) / Nbins;
102 ptdown = TMath::Power(10,ptmin + i * width );
103 ptup = TMath::Power(10,ptmin + (i+1) * width );
104 ptemp.set(ptdown, ptup, obj, xmin, xmax);
105 histos->push_back(ptemp);
106 }
107}
108
109//------------------------------------------------------------------------------
110
111class ExRootResult;
112class ExRootTreeReader;
113
114//------------------------------------------------------------------------------
115
116void BinLogX(TH1*h)
117{
118
119 TAxis *axis = h->GetXaxis();
120 int bins = axis->GetNbins();
121
122 Axis_t from = axis->GetXmin();
123 Axis_t to = axis->GetXmax();
124 Axis_t width = (to - from) / bins;
125 Axis_t *new_bins = new Axis_t[bins + 1];
126
127 for (int i = 0; i <= bins; i++) {
128 new_bins[i] = TMath::Power(10, from + i * width);
129
130 }
131 axis->Set(bins, new_bins);
132 delete new_bins;
133}
134
135
136//------------------------------------------------------------------------------
137
138template<typename T>
139std::pair<TH1D*, TH1D*> GetEff(TClonesArray *branchReco, TClonesArray *branchParticle, TString name, int pdgID, ExRootTreeReader *treeReader)
140{
141
142 cout << "** Computing Efficiency of reconstructing "<< branchReco->GetName() << " induced by " << branchParticle->GetName() << " with PID " << pdgID << endl;
143
144 Long64_t allEntries = treeReader->GetEntries();
145
146 GenParticle *particle;
147 T *recoObj;
148
149 TLorentzVector recoMomentum, genMomentum, bestRecoMomentum;
150
151 Float_t deltaR;
152 Float_t pt, eta;
153 Long64_t entry;
154
155 Int_t i, j;
156
157 TH1D *histGenPtcen = new TH1D(name+" gen spectra Pt",name+" gen spectra cen", Nbins, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax));
158 TH1D *histRecoPtcen = new TH1D(name+" reco spectra Pt",name+" reco spectra cen", Nbins, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax));
159 TH1D *histGenPtfwd = new TH1D(name+" gen spectra Eta",name+" gen spectra fwd", Nbins, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax));
160 TH1D *histRecoPtfwd = new TH1D(name+" reco spectra Eta",name+" reco spectra fwd", Nbins, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax));
161
162
163 BinLogX(histGenPtcen);
164 BinLogX(histRecoPtcen);
165 BinLogX(histGenPtfwd);
166 BinLogX(histRecoPtfwd);
167
168 // Loop over all events
169 for(entry = 0; entry < allEntries; ++entry)
170 {
171 // Load selected branches with data from specified event
172 treeReader->ReadEntry(entry);
173
174 // Loop over all generated particle in event
175 for(i = 0; i < branchParticle->GetEntriesFast(); ++i)
176 {
177
178 particle = (GenParticle*) branchParticle->At(i);
179 genMomentum = particle->P4();
180
181 deltaR = 999;
182
183 if (particle->PID == pdgID && genMomentum.Pt() > ptrangemin && genMomentum.Pt() < ptrangemax )
184 {
185
186 // Loop over all reco object in event
187 for(j = 0; j < branchReco->GetEntriesFast(); ++j)
188 {
189 recoObj = (T*)branchReco->At(j);
190 recoMomentum = recoObj->P4();
191 // this is simply to avoid warnings from initial state particle
192 // having infite rapidity ...
193 //if(Momentum.Px() == 0 && genMomentum.Py() == 0) continue;
194
195 // take the closest parton candidate
196 if(TMath::Abs(pdgID) == 5)
197 {
198 Jet *jet = (Jet *)recoObj;
199 if(jet->BTag != 1) continue;
200 }
201 if(TMath::Abs(pdgID) == 15)
202 {
203 Jet *jet = (Jet *)recoObj;
204 if(jet->TauTag != 1) continue;
205 }
206 if(genMomentum.DeltaR(recoMomentum) < deltaR)
207 {
208 deltaR = genMomentum.DeltaR(recoMomentum);
209 bestRecoMomentum = recoMomentum;
210 }
211 }
212
213 pt = genMomentum.Pt();
214 eta = genMomentum.Eta();
215
216 if (TMath::Abs(eta) < 1.5)
217 {
218 histGenPtcen->Fill(pt);
219 if(deltaR < 0.3) { histRecoPtcen->Fill(pt); }
220 }
221 else if (TMath::Abs(eta) < 2.5)
222 {
223 histGenPtfwd->Fill(pt);
224 if(deltaR < 0.3) { histRecoPtfwd->Fill(pt); }
225
226 }
227 }
228 }
229 }
230
231
232 std::pair<TH1D*,TH1D*> histos;
233
234 histRecoPtcen->Divide(histGenPtcen);
235 histRecoPtfwd->Divide(histGenPtfwd);
236
237 histos.first = histRecoPtcen;
238 histos.second = histRecoPtfwd;
239
240 return histos;
241}
242
243template<typename T>
244void GetEres(std::vector<resolPlot> *histos, TClonesArray *branchReco, TClonesArray *branchParticle, int pdgID, ExRootTreeReader *treeReader)
245{
246 Long64_t allEntries = treeReader->GetEntries();
247
248 cout << "** Computing resolution of " << branchReco->GetName() << " induced by " << branchParticle->GetName() << " with PID " << pdgID << endl;
249
250 GenParticle *particle;
251 T* recoObj;
252
253 TLorentzVector recoMomentum, genMomentum, bestGenMomentum;
254
255 Float_t deltaR;
256 Float_t pt, eta;
257 Long64_t entry;
258
259 Int_t i, j, bin;
260
261 // Loop over all events
262 for(entry = 0; entry < allEntries; ++entry)
263 {
264 // Load selected branches with data from specified event
265 treeReader->ReadEntry(entry);
266
267 // Loop over all reconstructed jets in event
268 for(i = 0; i < branchReco->GetEntriesFast(); ++i)
269 {
270 recoObj = (T*) branchReco->At(i);
271 recoMomentum = recoObj->P4();
272
273 deltaR = 999;
274
275 // Loop over all hard partons in event
276 for(j = 0; j < branchParticle->GetEntriesFast(); ++j)
277 {
278 particle = (GenParticle*) branchParticle->At(j);
279 if (particle->PID == pdgID && particle->Status == 1)
280 {
281 genMomentum = particle->P4();
282
283 // this is simply to avoid warnings from initial state particle
284 // having infite rapidity ...
285 if(genMomentum.Px() == 0 && genMomentum.Py() == 0) continue;
286
287 // take the closest parton candidate
288 if(genMomentum.DeltaR(recoMomentum) < deltaR)
289 {
290 deltaR = genMomentum.DeltaR(recoMomentum);
291 bestGenMomentum = genMomentum;
292 }
293 }
294 }
295
296 if(deltaR < 0.3)
297 {
298 pt = bestGenMomentum.Pt();
299 eta = TMath::Abs(bestGenMomentum.Eta());
300
301 for (bin = 0; bin < Nbins; bin++)
302 {
303 if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta > 0.0 && eta < 2.5)
304 {
305 if (eta < 1.5) {histos->at(bin).cenResolHist->Fill(recoMomentum.Pt()/bestGenMomentum.Pt());}
306 else if (eta < 2.5) {histos->at(bin).fwdResolHist->Fill(recoMomentum.Pt()/bestGenMomentum.Pt());}
307 }
308 }
309 }
310 }
311 }
312}
313void GetJetsEres(std::vector<resolPlot> *histos, TClonesArray *branchJet, TClonesArray *branchGenJet, ExRootTreeReader *treeReader)
314{
315
316 Long64_t allEntries = treeReader->GetEntries();
317
318 cout << "** Computing resolution of " << branchJet->GetName() << " induced by " << branchGenJet->GetName() << endl;
319
320 Jet *jet, *genjet;
321
322 TLorentzVector jetMomentum, genJetMomentum, bestGenJetMomentum;
323
324 Float_t deltaR;
325 Float_t pt, eta;
326 Long64_t entry;
327
328 Int_t i, j, bin;
329
330 // Loop over all events
331 for(entry = 0; entry < allEntries; ++entry)
332 {
333 // Load selected branches with data from specified event
334 treeReader->ReadEntry(entry);
335
336 if(entry%10000 == 0) cout << "Event number: "<< entry <<endl;
337
338 // Loop over all reconstructed jets in event
339 for(i = 0; i < TMath::Min(2,branchJet->GetEntriesFast()); ++i) //branchJet->GetEntriesFast(); ++i)
340 {
341
342 jet = (Jet*) branchJet->At(i);
343 jetMomentum = jet->P4();
344
345 deltaR = 999;
346
347 // Loop over all hard partons in event
348 for(j = 0; j < TMath::Min(2,branchGenJet->GetEntriesFast()); ++j)
349 {
350 genjet = (Jet*) branchGenJet->At(j);
351
352 genJetMomentum = genjet->P4();
353
354 // this is simply to avoid warnings from initial state particle
355 // having infite rapidity ...
356 if(genJetMomentum.Px() == 0 && genJetMomentum.Py() == 0) continue;
357
358 // take the closest parton candidate
359 if(genJetMomentum.DeltaR(jetMomentum) < deltaR)
360 {
361 deltaR = genJetMomentum.DeltaR(jetMomentum);
362 bestGenJetMomentum = genJetMomentum;
363 }
364 }
365
366 if(deltaR < 0.25)
367 {
368 pt = genJetMomentum.Pt();
369 eta = TMath::Abs(genJetMomentum.Eta());
370
371 for (bin = 0; bin < Nbins; bin++)
372 {
373 if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta < 1.5)
374 {
375 histos->at(bin).cenResolHist->Fill(jetMomentum.Pt()/bestGenJetMomentum.Pt());
376 }
377 else if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta < 2.5)
378 {
379 histos->at(bin).fwdResolHist->Fill(jetMomentum.Pt()/bestGenJetMomentum.Pt());
380 }
381
382 }
383 }
384 }
385 }
386}
387
388void GetMetres(std::vector<resolPlot> *histos, TClonesArray *branchScalarHT, TClonesArray *branchMet, TClonesArray *branchJet, ExRootTreeReader *treeReader)
389{
390
391 Long64_t allEntries = treeReader->GetEntries();
392
393 cout << "** Computing resolution of " << branchMet->GetName() << " vs " << branchScalarHT->GetName() << endl;
394
395 MissingET *met;
396 ScalarHT *scalarHT;
397
398 Long64_t entry;
399
400 Int_t bin;
401 Double_t ht;
402
403 Jet *jet;
404 TLorentzVector p1, p2;
405
406 // Loop over all events
407 for(entry = 0; entry < allEntries; ++entry)
408 {
409 // Load selected branches with data from specified event
410 treeReader->ReadEntry(entry);
411
412 if(entry%10000 == 0) cout << "Event number: "<< entry <<endl;
413
414 if (branchJet->GetEntriesFast() > 1)
415 {
416
417 jet = (Jet*) branchJet->At(0);
418 p1 = jet->P4();
419 jet = (Jet*) branchJet->At(1);
420 p2 = jet->P4();
421
422 met = (MissingET*) branchMet->At(0);
423 scalarHT = (ScalarHT*) branchScalarHT->At(0);
424 ht = scalarHT->HT;
425
426 if(p1.Pt() < 0.75*ht/2) continue;
427 if(p2.Pt() < 0.75*ht/2) continue;
428
429 for (bin = 0; bin < Nbins; bin++)
430 {
431 if(ht > histos->at(bin).ptmin && ht < histos->at(bin).ptmax )
432 {
433 histos->at(bin).cenResolHist->Fill(met->P4().Px());
434 histos->at(bin).fwdResolHist->Fill(met->P4().Py());
435 }
436 }
437 }
438 }
439}
440
441
442std::pair<Double_t, Double_t> GausFit(TH1* hist)
443{
444
445 TF1 *f1 = new TF1("f1", "gaus", hist->GetMean()-2*hist->GetRMS(), hist->GetMean()+2*hist->GetRMS());
446 hist->Fit("f1","RQ");
447 Double_t sig = f1->GetParameter(2);
448 Double_t sigErr = f1->GetParError(2);
449 delete f1;
450 return make_pair (sig, sigErr);
451
452 /*
453 int bin1 = hist->FindFirstBinAbove(hist->GetMaximum()/2);
454 int bin2 = hist->FindLastBinAbove(hist->GetMaximum()/2);
455 double fwhm = hist->GetBinCenter(bin2) - hist->GetBinCenter(bin1);
456
457 return make_pair (fwhm, fwhm);
458 //return make_pair (hist->GetRMS(), hist->GetRMSError());
459 */
460}
461
462
463TGraphErrors EresGraph(std::vector<resolPlot> *histos, bool central, bool rms = false)
464{
465 Int_t bin;
466 Int_t count = 0;
467 TGraphErrors gr = TGraphErrors(Nbins/2);
468 Double_t sig = 0;
469 Double_t sigErr = 0;
470 for (bin = 0; bin < Nbins; bin++)
471 {
472 if (central == true && histos->at(bin).cenResolHist->GetEntries() > 100)
473 {
474 std::pair<Double_t, Double_t> sigvalues = GausFit(histos->at(bin).cenResolHist);
475 if (rms == true)
476 {
477 gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.second);
478 gr.SetPointError(count,0, sigvalues.second); // to correct
479 }
480 else
481 {
482 gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.first);
483 gr.SetPointError(count,0, sigvalues.second);
484 }
485 count++;
486 }
487
488 else if (central == false && histos->at(bin).fwdResolHist->GetEntries() > 10)
489 {
490 std::pair<Double_t, Double_t> sigvalues = GausFit(histos->at(bin).fwdResolHist);
491 if (rms == true)
492 {
493 gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.second);
494 gr.SetPointError(count,0, sigvalues.second); // to correct
495 }
496 else
497 {
498 gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.first);
499 gr.SetPointError(count,0, sigvalues.second);
500 }
501 count++;
502 }
503
504 }
505 return gr;
506}
507
508
509//------------------------------------------------------------------------------
510
511
512// type 1 : object, 2 : track, 3 : tower
513
514void addGraph(TMultiGraph *mg, TGraphErrors *gr, TLegend *leg, int type)
515{
516
517 gr->SetLineWidth(2);
518
519 switch ( type )
520 {
521 case 1:
522 gr->SetLineColor(objColor);
523 gr->SetLineStyle(objStyle);
524 std::cout << "Adding " << gr->GetName() << std::endl;
525 mg->Add(gr);
526 leg->AddEntry(gr,"Reco","l");
527 break;
528
529 case 2:
530 gr->SetLineColor(trackColor);
531 gr->SetLineStyle(trackStyle);
532 mg->Add(gr);
533 leg->AddEntry(gr,"Track","l");
534 break;
535
536 case 3:
537 gr->SetLineColor(towerColor);
538 gr->SetLineStyle(towerStyle);
539 mg->Add(gr);
540 leg->AddEntry(gr,"Tower","l");
541 break;
542
543 case 0:
544 gr->SetLineColor(objColor);
545 gr->SetLineStyle(objStyle);
546 mg->Add(gr);
547 break;
548
549 default:
550 std::cout << "wrong type, possibles choices are Object, Track and Tower" << std::endl;
551 break;
552 }
553}
554
555void addHist(TH1D *h, TLegend *leg, int type)
556{
557 h->SetLineWidth(2);
558
559 switch ( type )
560 {
561 case 1:
562 h->SetLineColor(objColor);
563 h->SetLineStyle(objStyle);
564 leg->AddEntry(h,"Reco","l");
565 break;
566
567 case 2:
568 h->SetLineColor(trackColor);
569 h->SetLineStyle(trackStyle);
570 leg->AddEntry(h,"Track","l");
571 break;
572
573 case 3:
574 h->SetLineColor(towerColor);
575 h->SetLineStyle(towerStyle);
576 leg->AddEntry(h,"Tower","l");
577 break;
578
579 case 0:
580 h->SetLineColor(objColor);
581 h->SetLineStyle(objStyle);
582 break;
583
584 default:
585 std::cout << "wrong type, possibles choices are Object, Track and Tower" << std::endl;
586 break;
587 }
588}
589
590void DrawAxis(TMultiGraph *mg, TLegend *leg, double max)
591{
592 mg->SetMinimum(0.);
593 mg->SetMaximum(max);
594 mg->GetXaxis()->SetLimits(ptrangemin,ptrangemax);
595 mg->GetYaxis()->SetTitle("resolution");
596 mg->GetXaxis()->SetTitle("p_{T} [GeV]");
597 mg->GetYaxis()->SetTitleSize(0.07);
598 mg->GetXaxis()->SetTitleSize(0.07);
599 mg->GetYaxis()->SetLabelSize(0.06);
600 mg->GetXaxis()->SetLabelSize(0.06);
601 mg->GetYaxis()->SetLabelOffset(0.03);
602 mg->GetYaxis()->SetTitleOffset(1.4);
603 mg->GetXaxis()->SetTitleOffset(1.4);
604
605 mg->GetYaxis()->SetNdivisions(505);
606
607 leg->SetBorderSize(0);
608 leg->SetShadowColor(0);
609 leg->SetFillColor(0);
610 leg->SetFillStyle(0);
611
612 gStyle->SetOptTitle(0);
613 gPad->SetLogx();
614 gPad->SetBottomMargin(0.2);
615 gPad->SetLeftMargin(0.2);
616 gPad->Modified();
617 gPad->Update();
618
619}
620
621void DrawAxis(TH1D *h, TLegend *leg, int type)
622{
623
624 h->GetYaxis()->SetRangeUser(0,1.0);
625 if (type == 0) h->GetXaxis()->SetTitle("p_{T} [GeV]");
626 else h->GetXaxis()->SetTitle("#eta");
627 h->GetYaxis()->SetTitle("efficiency");
628 h->GetYaxis()->SetTitleSize(0.07);
629 h->GetXaxis()->SetTitleSize(0.07);
630 h->GetYaxis()->SetLabelSize(0.06);
631 h->GetXaxis()->SetLabelSize(0.06);
632 h->GetYaxis()->SetLabelOffset(0.03);
633 h->GetYaxis()->SetTitleOffset(1.3);
634 h->GetXaxis()->SetTitleOffset(1.4);
635
636 h->GetYaxis()->SetNdivisions(505);
637
638 leg->SetBorderSize(0);
639 leg->SetShadowColor(0);
640 leg->SetFillColor(0);
641 leg->SetFillStyle(0);
642
643 gStyle->SetOptTitle(0);
644 gStyle->SetOptStat(0);
645 gPad->SetBottomMargin(0.2);
646 gPad->SetLeftMargin(0.2);
647
648 gPad->Modified();
649 gPad->Update();
650
651}
652
653
654void Validation(const char *inputFile, const char *outputFile)
655{
656 //gSystem->Load("libDelphes");
657
658 std::cout << "input file : " << inputFile << " " << " , output file : " << outputFile << std::endl;
659
660 TChain *chain = new TChain("Delphes");
661 chain->Add(inputFile);
662
663 ExRootTreeReader *treeReader = new ExRootTreeReader(chain);
664 TClonesArray *branchParticle = treeReader->UseBranch("Particle");
665 TClonesArray *branchElectron = treeReader->UseBranch("Electron");
666 TClonesArray *branchMuon = treeReader->UseBranch("Muon");
667 TClonesArray *branchPhoton = treeReader->UseBranch("Photon");
668 TClonesArray *branchTrack = treeReader->UseBranch("Track");
669 TClonesArray *branchTower = treeReader->UseBranch("Tower");
670 TClonesArray *branchGenJet = treeReader->UseBranch("GenJet");
671 TClonesArray *branchPFJet = treeReader->UseBranch("Jet");
672 TClonesArray *branchCaloJet = treeReader->UseBranch("CaloJet");
673 TClonesArray *branchScalarHT = treeReader->UseBranch("ScalarHT");
674 TClonesArray *branchMet = treeReader->UseBranch("MissingET");
675
676
677#ifdef ELECTRON
678
679 ///////////////
680 // Electrons //
681 ///////////////
682
683 // Reconstruction efficiency
684 TString elecs = "Electron";
685 int elID = 11;
686 std::pair<TH1D*,TH1D*> histos_el = GetEff<Electron>(branchElectron, branchParticle, "Electron", elID, treeReader);
687
688 // tracking reconstruction efficiency
689 std::pair <TH1D*,TH1D*> histos_eltrack = GetEff<Track>(branchTrack, branchParticle, "electronTrack", elID, treeReader);
690
691 // Tower reconstruction efficiency
692 std::pair <TH1D*,TH1D*> histos_eltower = GetEff<Tower>(branchTower, branchParticle, "electronTower", elID, treeReader);
693
694 // Electron Energy Resolution
695 std::vector<resolPlot> plots_el;
696 HistogramsCollection(&plots_el, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "electrons");
697 GetEres<Electron>( &plots_el, branchElectron, branchParticle, elID, treeReader);
698 TGraphErrors gr_el = EresGraph(&plots_el, true);
699 TGraphErrors gr_elFwd = EresGraph(&plots_el, false);
700 gr_el.SetName("Electron");
701 gr_elFwd.SetName("ElectronFwd");
702
703 // Electron Track Energy Resolution
704 std::vector<resolPlot> plots_eltrack;
705 HistogramsCollection(&plots_eltrack, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "electronsTracks");
706 GetEres<Track>( &plots_eltrack, branchTrack, branchParticle, elID, treeReader);
707 TGraphErrors gr_eltrack = EresGraph(&plots_eltrack, true);
708 TGraphErrors gr_eltrackFwd = EresGraph(&plots_eltrack, false);
709 gr_eltrack.SetName("ElectronTracks");
710 gr_eltrackFwd.SetName("ElectronTracksFwd");
711
712 // Electron Tower Energy Resolution
713 std::vector<resolPlot> plots_eltower;
714 HistogramsCollection(&plots_eltower, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "electronsTower");
715 GetEres<Tower>( &plots_eltower, branchTower, branchParticle, elID, treeReader);
716 TGraphErrors gr_eltower = EresGraph(&plots_eltower, true);
717 TGraphErrors gr_eltowerFwd = EresGraph(&plots_eltower, false);
718 gr_eltower.SetName("ElectronTower");
719 gr_eltrackFwd.SetName("ElectronTracksFwd");
720
721 // Canvases
722 TString elEff = "electronEff";
723 TCanvas *C_el1 = new TCanvas(elEff,elEff, 1600, 600);
724 C_el1->Divide(2);
725 C_el1->cd(1);
726 TLegend *leg_el1 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
727 leg_el1->SetHeader("#splitline{electrons}{|#eta| < 1.5}");
728 leg_el1->AddEntry("","","");
729
730 gPad->SetLogx();
731 histos_eltrack.first->Draw("][");
732 addHist(histos_eltrack.first, leg_el1, 2);
733 //histos_eltower.first->Draw("same");
734 //addHist(histos_eltower.first, leg_el1, 3);
735 histos_el.first->Draw("same ][");
736 addHist(histos_el.first, leg_el1, 1);
737 DrawAxis(histos_eltrack.first, leg_el1, 0);
738
739 leg_el1->Draw();
740
741/*
742 TPaveText* txt = new TPaveText(effLegXmin,effLegYmax+0.02,effLegXmax,effLegYmax+0.1,"brNDC") ;
743 txt->AddText("electrons");
744 txt->SetBorderSize(0);
745 txt->SetShadowColor(0);
746 txt->SetFillColor(0);
747 txt->SetFillStyle(0);
748 txt->SetTextAlign(22);
749 txt->Draw();
750*/
751 C_el1->cd(2);
752 TLegend *leg_el2 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
753 leg_el2->SetHeader("#splitline{electrons}{1.5 < |#eta| < 2.5}");
754 leg_el2->AddEntry("","","");
755
756 gPad->SetLogx();
757 histos_eltrack.second->Draw("][");
758 addHist(histos_eltrack.second, leg_el2, 2);
759 //histos_eltower.second->Draw("same");
760 //addHist(histos_eltower.second, leg_el2, 3);
761 histos_el.second->Draw("same ][");
762 addHist(histos_el.second, leg_el2, 1);
763
764 DrawAxis(histos_eltrack.second, leg_el2, 0);
765 leg_el2->Draw();
766
767 //txt->Draw("same");
768
769 C_el1->cd(0);
770
771 TString elRes = "electronERes";
772 TString elResFwd = "electronEResForward";
773 TCanvas *C_el2 = new TCanvas(elRes,elRes, 1600, 600);
774 C_el2->Divide(2);
775 C_el2->cd(1);
776 TMultiGraph *mg_el = new TMultiGraph(elRes,elRes);
777 TLegend *leg_el = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
778 leg_el->SetHeader("#splitline{electrons}{|#eta| < 1.5}");
779 leg_el->AddEntry("","","");
780
781 addGraph(mg_el, &gr_eltower, leg_el, 3);
782 addGraph(mg_el, &gr_eltrack, leg_el, 2);
783 addGraph(mg_el, &gr_el, leg_el, 1);
784
785 mg_el->Draw("ACX");
786 leg_el->Draw();
787
788/*
789 TPaveText* txt2 = new TPaveText(0.72,0.57,0.95,0.65,"brNDC") ;
790
791 txt2->AddText("electrons");
792 txt2->SetBorderSize(0);
793 txt2->SetShadowColor(0);
794 txt2->SetFillColor(0);
795 txt2->SetFillStyle(0);
796 //txt2->SetTextAlign(12);
797 txt2->Draw();
798*/
799
800 DrawAxis(mg_el, leg_el, 0.1);
801
802 C_el2->cd(2);
803 TMultiGraph *mg_elFwd = new TMultiGraph(elResFwd,elResFwd);
804 TLegend *leg_elFwd = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
805 leg_elFwd->SetHeader("#splitline{electrons}{1.5 < |#eta| < 2.5}");
806 leg_elFwd->AddEntry("","","");
807
808
809 addGraph(mg_elFwd, &gr_eltowerFwd, leg_elFwd, 3);
810 addGraph(mg_elFwd, &gr_eltrackFwd, leg_elFwd, 2);
811 addGraph(mg_elFwd, &gr_elFwd, leg_elFwd, 1);
812
813 mg_elFwd->Draw("ACX");
814 leg_elFwd->Draw();
815
816 //txt2->Draw();
817
818 DrawAxis(mg_elFwd, leg_elFwd, 0.2);
819
820 C_el1->Print("electron.pdf(","pdf");
821 C_el2->Print("electron.pdf)","pdf");
822
823 //C_el1->SaveAs(elEff+".eps");
824 //C_el2->SaveAs(elRes+".eps");
825
826#endif
827
828 gDirectory->cd(0);
829
830#ifdef MUON
831
832 ///////////
833 // Muons //
834 ///////////
835
836 // Reconstruction efficiency
837 int muID = 13;
838 std::pair<TH1D*,TH1D*> histos_mu = GetEff<Muon>(branchMuon, branchParticle,"Muon", muID, treeReader);
839
840 // muon tracking reconstruction efficiency
841 std::pair <TH1D*,TH1D*> histos_mutrack = GetEff<Track>(branchTrack, branchParticle, "muonTrack", muID, treeReader);
842
843 // Muon Energy Resolution
844 std::vector<resolPlot> plots_mu;
845 HistogramsCollection(&plots_mu, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "muons");
846 GetEres<Muon>( &plots_mu, branchMuon, branchParticle, muID, treeReader);
847 TGraphErrors gr_mu = EresGraph(&plots_mu, true);
848 TGraphErrors gr_muFwd = EresGraph(&plots_mu, false);
849 gr_mu.SetName("Muon");
850 gr_muFwd.SetName("MuonFwd");
851
852 // Muon Track Energy Resolution
853 std::vector<resolPlot> plots_mutrack;
854 HistogramsCollection(&plots_mutrack, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "muonsTracks");
855 GetEres<Track>( &plots_mutrack, branchTrack, branchParticle, muID, treeReader);
856 TGraphErrors gr_mutrack = EresGraph(&plots_mutrack, true);
857 TGraphErrors gr_mutrackFwd = EresGraph(&plots_mutrack, false);
858 gr_mutrackFwd.SetName("MuonTracksFwd");
859
860 // Canvas
861
862 TString muEff = "muonEff";
863 TCanvas *C_mu1 = new TCanvas(muEff,muEff, 1600, 600);
864 C_mu1->Divide(2);
865 C_mu1->cd(1);
866 TLegend *leg_mu1 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
867 leg_mu1->SetHeader("#splitline{muons}{|#eta| < 1.5}");
868 leg_mu1->AddEntry("","","");
869
870
871 gPad->SetLogx();
872 histos_mutrack.first->Draw("][");
873 addHist(histos_mutrack.first, leg_mu1, 2);
874 histos_mu.first->Draw("same ][");
875 addHist(histos_mu.first, leg_mu1, 1);
876
877 DrawAxis(histos_mutrack.first, leg_mu1, 0);
878
879 leg_mu1->Draw();
880
881 C_mu1->cd(2);
882 TLegend *leg_mu2 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
883 leg_mu2->SetHeader("#splitline{muons}{1.5 < |#eta| < 2.5}");
884 leg_mu2->AddEntry("","","");
885
886 gPad->SetLogx();
887 histos_mutrack.second->Draw("][");
888 addHist(histos_mutrack.second, leg_mu2, 2);
889 histos_mu.second->Draw("same ][");
890 addHist(histos_mu.second, leg_mu2, 1);
891
892 DrawAxis(histos_mutrack.second, leg_mu2, 1);
893 leg_mu2->Draw();
894
895 TString muRes = "muonERes";
896 TString muResFwd = "muonEResFwd";
897
898 TCanvas *C_mu = new TCanvas(muRes,muRes, 1600, 600);
899 C_mu->Divide(2);
900 C_mu->cd(1);
901 TMultiGraph *mg_mu = new TMultiGraph(muRes,muRes);
902 TLegend *leg_mu = new TLegend(topLeftLegXmin,topLeftLegYmin,topLeftLegXmax,topLeftLegYmax);
903 leg_mu->SetHeader("#splitline{muons}{|#eta| < 1.5}");
904 leg_mu->AddEntry("","","");
905
906 addGraph(mg_mu, &gr_mutrack, leg_mu, 2);
907 addGraph(mg_mu, &gr_mu, leg_mu, 1);
908
909 mg_mu->Draw("ACX");
910 leg_mu->Draw();
911
912 DrawAxis(mg_mu, leg_mu, 0.3);
913
914 C_mu->cd(2);
915 TMultiGraph *mg_muFwd = new TMultiGraph(muResFwd,muResFwd);
916 TLegend *leg_muFwd = new TLegend(topLeftLegXmin,topLeftLegYmin,topLeftLegXmax,topLeftLegYmax);
917 leg_muFwd->SetHeader("#splitline{muons}{1.5 < |#eta| < 2.5}");
918 leg_muFwd->AddEntry("","","");
919
920 addGraph(mg_muFwd, &gr_mutrackFwd, leg_muFwd, 2);
921 addGraph(mg_muFwd, &gr_muFwd, leg_muFwd, 1);
922
923 mg_muFwd->Draw("ACX");
924 leg_muFwd->Draw();
925
926 DrawAxis(mg_muFwd, leg_muFwd, 0.3);
927
928 //C_mu1->SaveAs(muEff+".eps");
929 //C_mu->SaveAs(muRes+".eps");
930
931 C_mu1->Print("muon.pdf(","pdf");
932 C_mu->Print("muon.pdf)","pdf");
933
934#endif
935
936 gDirectory->cd(0);
937
938#ifdef PHOTON
939
940 /////////////
941 // Photons //
942 /////////////
943
944 // Reconstruction efficiency
945 int phID = 22;
946 std::pair<TH1D*,TH1D*> histos_ph = GetEff<Electron>(branchPhoton, branchParticle, "Photon", phID, treeReader);
947 std::pair<TH1D*,TH1D*> histos_phtower = GetEff<Electron>(branchTower, branchParticle, "Photon", phID, treeReader);
948
949 // Photon Energy Resolution
950 std::vector<resolPlot> plots_ph;
951 HistogramsCollection(&plots_ph, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "photons");
952 GetEres<Photon>( &plots_ph, branchPhoton, branchParticle, phID, treeReader);
953 TGraphErrors gr_ph = EresGraph(&plots_ph, true);
954 TGraphErrors gr_phFwd = EresGraph(&plots_ph, false);
955 gr_ph.SetName("Photon");
956 gr_phFwd.SetName("PhotonFwd");
957
958
959 // Photon Tower Energy Resolution
960 std::vector<resolPlot> plots_phtower;
961 HistogramsCollection(&plots_phtower, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "photonsTower");
962 GetEres<Tower>( &plots_phtower, branchTower, branchParticle, phID, treeReader);
963 TGraphErrors gr_phtower = EresGraph(&plots_phtower, true);
964 TGraphErrors gr_phtowerFwd = EresGraph(&plots_phtower, false);
965 gr_phtower.SetName("PhotonTower");
966 gr_phtowerFwd.SetName("PhotonTowerFwd");
967
968 // Canvas
969
970 TString phEff = "photonEff";
971 TCanvas *C_ph1 = new TCanvas(phEff,phEff, 1600, 600);
972 C_ph1->Divide(2);
973 C_ph1->cd(1);
974 TLegend *leg_ph1 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
975 leg_ph1->SetHeader("#splitline{photons}{|#eta| < 1.5}");
976 leg_ph1->AddEntry("","","");
977
978
979 gPad->SetLogx();
980 histos_phtower.first->Draw("][");
981 addHist(histos_phtower.first, leg_ph1, 3);
982 histos_ph.first->Draw("same ][");
983 addHist(histos_ph.first, leg_ph1, 1);
984
985 DrawAxis(histos_phtower.first, leg_ph1, 0);
986 leg_ph1->Draw();
987
988 C_ph1->cd(2);
989 TLegend *leg_ph2 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
990 leg_ph2->SetHeader("#splitline{photons}{1.5 < |#eta| < 2.5}");
991 leg_ph2->AddEntry("","","");
992
993
994 gPad->SetLogx();
995 histos_phtower.second->Draw("][");
996 addHist(histos_phtower.second, leg_ph2, 3);
997 histos_ph.second->Draw("same ][");
998 addHist(histos_ph.second, leg_ph2, 1);
999
1000 DrawAxis(histos_phtower.second, leg_ph2, 1);
1001 leg_ph2->Draw();
1002
1003 C_ph1->SaveAs(phEff+".eps");
1004
1005 TString phRes = "phERes";
1006 TString phResFwd = "phEResFwd";
1007
1008 TCanvas *C_ph = new TCanvas(phRes,phRes, 1600, 600);
1009 C_ph->Divide(2);
1010 C_ph->cd(1);
1011 TMultiGraph *mg_ph = new TMultiGraph(phRes,phRes);
1012 TLegend *leg_ph = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
1013 leg_ph->SetHeader("#splitline{photons}{|#eta| < 1.5}");
1014 leg_ph->AddEntry("","","");
1015
1016 addGraph(mg_ph, &gr_phtower, leg_ph, 3);
1017 addGraph(mg_ph, &gr_ph, leg_ph, 1);
1018
1019 mg_ph->Draw("ACX");
1020 leg_ph->Draw();
1021
1022 DrawAxis(mg_ph, leg_ph, 0.3);
1023
1024 C_ph->cd(2);
1025 TMultiGraph *mg_phFwd = new TMultiGraph(phResFwd,phResFwd);
1026 TLegend *leg_phFwd = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
1027 leg_phFwd->SetHeader("#splitline{photons}{1.5 < |#eta| < 2.5}");
1028 leg_phFwd->AddEntry("","","");
1029
1030 addGraph(mg_phFwd, &gr_phtowerFwd, leg_phFwd, 3);
1031 addGraph(mg_phFwd, &gr_phFwd, leg_phFwd, 1);
1032
1033 mg_phFwd->Draw("ACX");
1034 leg_phFwd->Draw();
1035
1036 DrawAxis(mg_phFwd, leg_phFwd, 0.3);
1037
1038 C_ph->SaveAs(phRes+".eps");
1039
1040 C_ph1->Print("photon.pdf(","pdf");
1041 C_ph->Print("photon.pdf)","pdf");
1042
1043
1044#endif
1045
1046 gDirectory->cd(0);
1047
1048#ifdef JET
1049
1050 //////////
1051 // Jets //
1052 //////////
1053
1054 // BJets Reconstruction efficiency
1055 int bID = 5;
1056 std::pair<TH1D*,TH1D*> histos_btag = GetEff<Jet>(branchPFJet, branchParticle,"BTag", bID, treeReader);
1057
1058 // TauJets Reconstruction efficiency
1059 int tauID = 15;
1060 std::pair<TH1D*,TH1D*> histos_tautag = GetEff<Jet>(branchPFJet, branchParticle,"TauTag", tauID, treeReader);
1061
1062 // PFJets Energy Resolution
1063 std::vector<resolPlot> plots_pfjets;
1064 HistogramsCollection(&plots_pfjets, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "PFJet");
1065 GetJetsEres( &plots_pfjets, branchPFJet, branchGenJet, treeReader);
1066 TGraphErrors gr_pfjets = EresGraph(&plots_pfjets, true);
1067 TGraphErrors gr_pfjetsFwd = EresGraph(&plots_pfjets, false);
1068 gr_pfjets.SetName("pfJet");
1069 gr_pfjetsFwd.SetName("pfJetFwd");
1070
1071 // CaloJets Energy Resolution
1072 std::vector<resolPlot> plots_calojets;
1073 HistogramsCollection(&plots_calojets, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "CaloJet");
1074 GetJetsEres( &plots_calojets, branchCaloJet, branchGenJet, treeReader);
1075 TGraphErrors gr_calojets = EresGraph(&plots_calojets, true);
1076 TGraphErrors gr_calojetsFwd = EresGraph(&plots_calojets, false);
1077 gr_calojets.SetName("caloJet");
1078 gr_calojetsFwd.SetName("caloJetFwd");
1079
1080 // MET Resolution vs HT
1081 std::vector<resolPlot> plots_met;
1082 HistogramsCollection(&plots_met, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "MET", -500, 500);
1083 GetMetres( &plots_met, branchScalarHT, branchMet, branchPFJet, treeReader);
1084 TGraphErrors gr_met = EresGraph(&plots_met, true);
1085 gr_calojets.SetName("MET");
1086
1087 // Canvas
1088 TString btagEff = "btagEff";
1089 TCanvas *C_btag1 = new TCanvas(btagEff,btagEff, 1600, 600);
1090 C_btag1->Divide(2);
1091 C_btag1->cd(1);
1092 TLegend *leg_btag1 = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
1093 leg_btag1->SetHeader("#splitline{B-tagging}{|#eta| < 1.5}");
1094 leg_btag1->AddEntry("","","");
1095
1096 gPad->SetLogx();
1097 histos_btag.first->Draw();
1098 addHist(histos_btag.first, leg_btag1, 0);
1099
1100 DrawAxis(histos_btag.first, leg_btag1, 0);
1101 leg_btag1->Draw();
1102
1103 C_btag1->cd(2);
1104 TLegend *leg_btag2 = new TLegend(resLegXmin,resLegYmin,resLegXmax+0.05,resLegYmax);
1105 leg_btag2->SetHeader("#splitline{B-tagging}{1.5 < |#eta| < 2.5}");
1106 leg_btag2->AddEntry("","","");
1107
1108 gPad->SetLogx();
1109 histos_btag.second->Draw();
1110 addHist(histos_btag.second, leg_btag2, 0);
1111
1112 DrawAxis(histos_btag.second, leg_btag2, 0);
1113 leg_btag2->Draw();
1114 C_btag1->cd(0);
1115
1116 TString tautagEff = "tautagEff";
1117 TCanvas *C_tautag1 = new TCanvas(tautagEff,tautagEff, 1600, 600);
1118 C_tautag1->Divide(2);
1119 C_tautag1->cd(1);
1120 TLegend *leg_tautag1 = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
1121 leg_tautag1->SetHeader("#splitline{#tau-tagging}{|#eta| < 1.5}");
1122 leg_tautag1->AddEntry("","","");
1123
1124 gPad->SetLogx();
1125 histos_tautag.first->Draw();
1126 addHist(histos_tautag.first, leg_tautag1, 0);
1127
1128 DrawAxis(histos_tautag.first, leg_tautag1, 0);
1129 leg_tautag1->Draw();
1130
1131 C_tautag1->cd(2);
1132 TLegend *leg_tautag2 = new TLegend(resLegXmin,resLegYmin,resLegXmax+0.05,resLegYmax);
1133 leg_tautag2->SetHeader("#splitline{#tau-tagging}{1.5 < |#eta| < 2.5}");
1134 leg_tautag2->AddEntry("","","");
1135
1136 gPad->SetLogx();
1137 histos_tautag.second->Draw();
1138 addHist(histos_tautag.second, leg_tautag2, 0);
1139
1140 DrawAxis(histos_tautag.second, leg_tautag2, 0);
1141 leg_tautag2->Draw();
1142 C_tautag1->cd(0);
1143
1144 TString jetRes = "jetERes";
1145 TString jetResFwd = "jetEResFwd";
1146 TCanvas *C_jet = new TCanvas(jetRes,jetRes, 1600, 600);
1147 C_jet->Divide(2);
1148
1149 C_jet->cd(1);
1150 TMultiGraph *mg_jet = new TMultiGraph(jetRes,jetRes);
1151 TLegend *leg_jet = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
1152 leg_jet->SetHeader("#splitline{jets}{|#eta| < 1.5}");
1153 leg_jet->AddEntry("","","");
1154
1155 addGraph(mg_jet, &gr_calojets, leg_jet, 3);
1156 addGraph(mg_jet, &gr_pfjets, leg_jet, 1);
1157
1158 mg_jet->Draw("ACX");
1159 leg_jet->Draw();
1160
1161 DrawAxis(mg_jet, leg_jet, 0.25);
1162
1163 C_jet->cd(2);
1164 TMultiGraph *mg_jetFwd = new TMultiGraph(jetResFwd,jetResFwd);
1165 TLegend *leg_jetFwd = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
1166 leg_jetFwd->SetHeader("#splitline{jets}{|#eta| < 1.5}");
1167 leg_jetFwd->AddEntry("","","");
1168
1169 addGraph(mg_jetFwd, &gr_calojetsFwd, leg_jetFwd, 3);
1170 addGraph(mg_jetFwd, &gr_pfjetsFwd, leg_jetFwd, 1);
1171
1172 mg_jetFwd->Draw("ACX");
1173 leg_jetFwd->Draw();
1174
1175 DrawAxis(mg_jetFwd, leg_jetFwd, 0.25);
1176
1177 C_btag1->SaveAs(btagEff+".eps");
1178 C_jet->SaveAs(jetRes+".eps");
1179
1180 TString metRes = "MetRes";
1181 TCanvas *C_met = new TCanvas(metRes,metRes, 800, 600);
1182
1183 TMultiGraph *mg_met = new TMultiGraph(metRes,metRes);
1184 TLegend *leg_met = new TLegend(topLeftLegXmin,topLeftLegYmin+0.2,topLeftLegXmax-0.2,topLeftLegYmax);
1185 leg_met->SetHeader("E_{T}^{miss}");
1186 leg_met->AddEntry("","","");
1187
1188
1189 addGraph(mg_met, &gr_met, leg_met, 0);
1190
1191 mg_met->Draw("ACX");
1192 leg_met->Draw();
1193
1194 DrawAxis(mg_met, leg_met, 300);
1195
1196 mg_met->GetXaxis()->SetTitle("H_{T} [GeV]");
1197 mg_met->GetYaxis()->SetTitle("#sigma(ME_{x}) [GeV]");
1198
1199 C_met->SaveAs(metRes+".eps");
1200
1201 C_jet->Print("jet.pdf(","pdf");
1202 C_btag1->Print("jet.pdf","pdf");
1203 C_tautag1->Print("jet.pdf","pdf");
1204 C_met->Print("jet.pdf)","pdf");
1205
1206
1207#endif
1208
1209 /*
1210 // CaloJets Energy Resolution
1211 std::vector<resolPlot> plots_calojets;
1212 HistogramsCollection(&plots_calojets, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "caloJet");
1213 GetJetsEres( &plots_calojets, branchCaloJet, branchGenJet, treeReader);
1214 TGraphErrors gr_calojets = EresGraph(&plots_calojets, true);
1215 gr_calojets.SetName("caloJet");
1216 */
1217
1218 TFile *fout = new TFile(outputFile,"recreate");
1219
1220#ifdef ELECTRON
1221 for (int bin = 0; bin < Nbins; bin++)
1222 {
1223 plots_el.at(bin).cenResolHist->Write();
1224 plots_eltrack.at(bin).cenResolHist->Write();
1225 plots_eltower.at(bin).cenResolHist->Write();
1226 plots_el.at(bin).fwdResolHist->Write();
1227 plots_eltrack.at(bin).fwdResolHist->Write();
1228 plots_eltower.at(bin).fwdResolHist->Write();
1229 }
1230 // gr.Write();
1231 histos_el.first->Write();
1232 //histos_el.second->Write();
1233 histos_eltrack.first->Write();
1234 //histos_eltrack.second->Write();
1235 histos_eltower.first->Write();
1236 C_el1->Write();
1237 C_el2->Write();
1238#endif
1239
1240#ifdef MUON
1241 histos_mu.first->Write();
1242 histos_mu.second->Write();
1243 histos_mutrack.first->Write();
1244 histos_mutrack.second->Write();
1245 C_mu1->Write();
1246 C_mu->Write();
1247#endif
1248
1249#ifdef PHOTON
1250 histos_ph.first->Write();
1251 histos_ph.second->Write();
1252 C_ph1->Write();
1253 C_ph->Write();
1254#endif
1255
1256#ifdef JET
1257 for (int bin = 0; bin < Nbins; bin++)
1258 {
1259 plots_pfjets.at(bin).cenResolHist->Write();
1260 plots_pfjets.at(bin).fwdResolHist->Write();
1261 plots_calojets.at(bin).cenResolHist->Write();
1262 plots_calojets.at(bin).fwdResolHist->Write();
1263 plots_met.at(bin).cenResolHist->Write();
1264 }
1265 histos_btag.first->Write();
1266 histos_btag.second->Write();
1267 histos_tautag.first->Write();
1268 histos_tautag.second->Write();
1269 C_btag1->Write();
1270 C_tautag1->Write();
1271 C_jet->Write();
1272 C_met->Write();
1273#endif
1274
1275 fout->Write();
1276
1277 cout << "** Exiting..." << endl;
1278
1279 //delete plots;
1280 //delete result;
1281 delete treeReader;
1282 delete chain;
1283}
1284
1285//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.