Fork me on GitHub

source: git/examples/Validation.C@ 9f8d4e7

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 9f8d4e7 was cf88574, checked in by Pavel Demin <pavel.demin@…>, 8 years ago

update validation code
  • Property mode set to 100644
File size: 27.2 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 = 5;
35
36Color_t objColor = kBlack;
37Color_t trackColor = kBlack;
38Color_t towerColor = kBlack;
39
40double effLegXmin = 0.22;
41double effLegXmax = 0.5;
42double effLegYmin = 0.22;
43double effLegYmax = 0.4;
44
45struct resolPlot
46{
47 TH1 *cenResolHist;
48 TH1 *fwdResolHist;
49 double ptmin;
50 double ptmax;
51 TString obj;
52
53 resolPlot();
54 resolPlot(double ptdown, double ptup, TString object);
55 void set(double ptdown, double ptup, TString object);
56 void print(){std::cout << ptmin << std::endl;}
57};
58
59
60resolPlot::resolPlot()
61{
62}
63
64resolPlot::resolPlot(double ptdown, double ptup, TString object)
65{
66 this->set(ptdown,ptup,object);
67}
68
69void resolPlot::set(double ptdown, double ptup, TString object){
70 ptmin = ptdown;
71 ptmax = ptup;
72 obj = object;
73
74 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", 500, 0, 2);
75 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", 500, 0.4, 0.4);
76
77}
78
79void HistogramsCollection(std::vector<resolPlot> *histos, double ptmin, double ptmax, TString obj)
80{
81 double width;
82 double ptdown;
83 double ptup;
84 resolPlot ptemp;
85
86 for (int i = 0; i < Nbins; i++)
87 {
88 width = (ptmax - ptmin) / Nbins;
89 ptdown = TMath::Power(10,ptmin + i * width );
90 ptup = TMath::Power(10,ptmin + (i+1) * width );
91 ptemp.set(ptdown, ptup, obj);
92 histos->push_back(ptemp);
93 }
94}
95
96//------------------------------------------------------------------------------
97
98class ExRootResult;
99class ExRootTreeReader;
100
101//------------------------------------------------------------------------------
102
103void BinLogX(TH1*h)
104{
105
106 TAxis *axis = h->GetXaxis();
107 int bins = axis->GetNbins();
108
109 Axis_t from = axis->GetXmin();
110 Axis_t to = axis->GetXmax();
111 Axis_t width = (to - from) / bins;
112 Axis_t *new_bins = new Axis_t[bins + 1];
113
114 for (int i = 0; i <= bins; i++) {
115 new_bins[i] = TMath::Power(10, from + i * width);
116
117 }
118 axis->Set(bins, new_bins);
119 delete new_bins;
120}
121
122
123//------------------------------------------------------------------------------
124
125template<typename T>
126std::pair<TH1D*, TH1D*> GetEff(TClonesArray *branchReco, TClonesArray *branchParticle, TString name, int pdgID, ExRootTreeReader *treeReader)
127{
128
129 cout << "** Computing Efficiency of reconstructing "<< branchReco->GetName() << " induced by " << branchParticle->GetName() << " with PID " << pdgID << endl;
130
131 Long64_t allEntries = treeReader->GetEntries();
132
133 GenParticle *particle;
134 T *recoObj;
135
136 TLorentzVector recoMomentum, genMomentum, bestRecoMomentum;
137
138 Float_t deltaR;
139 Float_t pt, eta;
140 Long64_t entry;
141
142 Int_t i, j;
143
144 TH1D *histGenPt = new TH1D(name+" gen spectra Pt",name+" gen spectra Pt", Nbins, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax));
145 TH1D *histRecoPt = new TH1D(name+" reco spectra Pt",name+" reco spectra Pt", Nbins, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax));
146 TH1D *histGenEta = new TH1D(name+" gen spectra Eta",name+" gen spectra Eta", 12, -3, 3);
147 TH1D *histRecoEta = new TH1D(name+" reco spectra Eta",name+" reco spectra Eta", 12, -3, 3);
148
149
150 BinLogX(histGenPt);
151 BinLogX(histRecoPt);
152
153 // Loop over all events
154 for(entry = 0; entry < allEntries; ++entry)
155 {
156 // Load selected branches with data from specified event
157 treeReader->ReadEntry(entry);
158
159 // Loop over all generated particle in event
160 for(i = 0; i < branchParticle->GetEntriesFast(); ++i)
161 {
162
163 particle = (GenParticle*) branchParticle->At(i);
164 genMomentum = particle->P4();
165
166 deltaR = 999;
167
168 if (particle->PID == pdgID && genMomentum.Pt() > ptrangemin && genMomentum.Pt() < ptrangemax )
169 {
170
171 // Loop over all reco object in event
172 for(j = 0; j < branchReco->GetEntriesFast(); ++j)
173 {
174 recoObj = (T*)branchReco->At(j);
175 recoMomentum = recoObj->P4();
176 // this is simply to avoid warnings from initial state particle
177 // having infite rapidity ...
178 //if(Momentum.Px() == 0 && genMomentum.Py() == 0) continue;
179
180 // take the closest parton candidate
181 if(TMath::Abs(pdgID) == 5)
182 {
183 Jet *jet = (Jet *)recoObj;
184 if(jet->BTag != 1) continue;
185 }
186 if(TMath::Abs(pdgID) == 15)
187 {
188 Jet *jet = (Jet *)recoObj;
189 if(jet->TauTag != 1) continue;
190 }
191 if(genMomentum.DeltaR(recoMomentum) < deltaR)
192 {
193 deltaR = genMomentum.DeltaR(recoMomentum);
194 bestRecoMomentum = recoMomentum;
195 }
196 }
197
198 pt = genMomentum.Pt();
199 eta = genMomentum.Eta();
200
201 histGenPt->Fill(pt);
202 histGenEta->Fill(eta);
203
204 if(deltaR < 0.3)
205 {
206 histRecoPt->Fill(pt);
207 histRecoEta->Fill(eta);
208 }
209 }
210 }
211 }
212
213
214 std::pair<TH1D*,TH1D*> histos;
215
216 histRecoPt->Divide(histGenPt);
217 histRecoEta->Divide(histGenEta);
218
219 histos.first = histRecoPt;
220 histos.second = histRecoEta;
221
222 return histos;
223}
224
225template<typename T>
226void GetEres(std::vector<resolPlot> *histos, TClonesArray *branchReco, TClonesArray *branchParticle, int pdgID, ExRootTreeReader *treeReader)
227{
228 Long64_t allEntries = treeReader->GetEntries();
229
230 cout << "** Computing resolution of " << branchReco->GetName() << " induced by " << branchParticle->GetName() << " with PID " << pdgID << endl;
231
232 GenParticle *particle;
233 T* recoObj;
234
235 TLorentzVector recoMomentum, genMomentum, bestGenMomentum;
236
237 Float_t deltaR;
238 Float_t pt, eta;
239 Long64_t entry;
240
241 Int_t i, j, bin;
242
243 // Loop over all events
244 for(entry = 0; entry < allEntries; ++entry)
245 {
246 // Load selected branches with data from specified event
247 treeReader->ReadEntry(entry);
248
249 // Loop over all reconstructed jets in event
250 for(i = 0; i < branchReco->GetEntriesFast(); ++i)
251 {
252 recoObj = (T*) branchReco->At(i);
253 recoMomentum = recoObj->P4();
254
255 deltaR = 999;
256
257 // Loop over all hard partons in event
258 for(j = 0; j < branchParticle->GetEntriesFast(); ++j)
259 {
260 particle = (GenParticle*) branchParticle->At(j);
261 if (particle->PID == pdgID && particle->Status == 1)
262 {
263 genMomentum = particle->P4();
264
265 // this is simply to avoid warnings from initial state particle
266 // having infite rapidity ...
267 if(genMomentum.Px() == 0 && genMomentum.Py() == 0) continue;
268
269 // take the closest parton candidate
270 if(genMomentum.DeltaR(recoMomentum) < deltaR)
271 {
272 deltaR = genMomentum.DeltaR(recoMomentum);
273 bestGenMomentum = genMomentum;
274 }
275 }
276 }
277
278 if(deltaR < 0.3)
279 {
280 pt = bestGenMomentum.Pt();
281 eta = TMath::Abs(bestGenMomentum.Eta());
282
283 for (bin = 0; bin < Nbins; bin++)
284 {
285 if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta > 0.0 && eta < 2.5)
286 {
287 if (eta < 1.5) {histos->at(bin).cenResolHist->Fill(recoMomentum.Pt()/bestGenMomentum.Pt());}
288 else if (eta < 2.5) {histos->at(bin).fwdResolHist->Fill(recoMomentum.Pt()/bestGenMomentum.Pt());}
289 }
290 }
291 }
292 }
293 }
294}
295void GetJetsEres(std::vector<resolPlot> *histos, TClonesArray *branchJet, TClonesArray *branchGenJet, ExRootTreeReader *treeReader)
296{
297
298 Long64_t allEntries = treeReader->GetEntries();
299
300 cout << "** Computing resolution of " << branchJet->GetName() << " induced by " << branchGenJet->GetName() << endl;
301
302 Jet *jet, *genjet;
303
304 TLorentzVector jetMomentum, genJetMomentum, bestGenJetMomentum;
305
306 Float_t deltaR;
307 Float_t pt, eta;
308 Long64_t entry;
309
310 Int_t i, j, bin;
311
312 // Loop over all events
313 for(entry = 0; entry < allEntries; ++entry)
314 {
315 // Load selected branches with data from specified event
316 treeReader->ReadEntry(entry);
317
318 if(entry%10000 == 0) cout << "Event number: "<< entry <<endl;
319
320 // Loop over all reconstructed jets in event
321 for(i = 0; i < TMath::Min(2,branchJet->GetEntriesFast()); ++i) //branchJet->GetEntriesFast(); ++i)
322 {
323
324 jet = (Jet*) branchJet->At(i);
325 jetMomentum = jet->P4();
326
327 deltaR = 999;
328
329 // Loop over all hard partons in event
330 for(j = 0; j < TMath::Min(2,branchGenJet->GetEntriesFast()); ++j)
331 {
332 genjet = (Jet*) branchGenJet->At(j);
333
334 genJetMomentum = genjet->P4();
335
336 // this is simply to avoid warnings from initial state particle
337 // having infite rapidity ...
338 if(genJetMomentum.Px() == 0 && genJetMomentum.Py() == 0) continue;
339
340 // take the closest parton candidate
341 if(genJetMomentum.DeltaR(jetMomentum) < deltaR)
342 {
343 deltaR = genJetMomentum.DeltaR(jetMomentum);
344 bestGenJetMomentum = genJetMomentum;
345 }
346 }
347
348 if(deltaR < 0.25)
349 {
350 pt = genJetMomentum.Pt();
351 eta = TMath::Abs(genJetMomentum.Eta());
352
353 for (bin = 0; bin < Nbins; bin++)
354 {
355 if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta > 0.0 && eta < 1.5)
356 {
357 histos->at(bin).cenResolHist->Fill(jetMomentum.Pt()/bestGenJetMomentum.Pt());
358 }
359 }
360 }
361 }
362 }
363}
364std::pair<Double_t, Double_t> GausFit(TH1* hist)
365{
366 TF1 *f1 = new TF1("f1", "gaus", hist->GetMean()-2*hist->GetRMS(), hist->GetMean()+2*hist->GetRMS());
367 hist->Fit("f1","RQ");
368 Double_t sig = f1->GetParameter(2);
369 Double_t sigErr = f1->GetParError(2);
370 delete f1;
371 return make_pair (sig, sigErr);
372 //return make_pair (hist->GetRMS(), hist->GetRMSError());
373}
374
375
376TGraphErrors EresGraph(std::vector<resolPlot> *histos, bool central)
377{
378 Int_t bin;
379 Int_t count = 0;
380 TGraphErrors gr = TGraphErrors(Nbins/2);
381 Double_t sig = 0;
382 Double_t sigErr = 0;
383 for (bin = 0; bin < Nbins; bin++)
384 {
385 if (central == true && histos->at(bin).cenResolHist->GetEntries() > 100)
386 {
387 std::cout << " pt : " << (histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0;
388 std::cout << " mean : " << histos->at(bin).cenResolHist->GetMean() << " RMS : " << histos->at(bin).cenResolHist->GetRMS();
389 std::cout << " entries : " << histos->at(bin).cenResolHist->GetEntries() << std::endl;
390 std::pair<Double_t, Double_t> sigvalues = GausFit(histos->at(bin).cenResolHist);
391 gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.first);
392 gr.SetPointError(count,0, sigvalues.second);
393 count++;
394 }
395 /*
396 else if (histos->at(bin).cenResolHist->GetEntries() > 10)
397 {
398 histos->at(bin).fwdResolHist->Fit("gaus","","", -2*histos->at(bin).fwdResolHist->GetRMS(), 2*histos->at(bin).fwdResolHist->GetRMS());
399 TF1 *f = histos->at(bin).fwdResolHist->GetFunction("gaus");
400 Double_t sig = f->GetParameter(2);
401 Double_t sigErr = f->GetParError(2);
402 gr.SetPoint(bin,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sig);
403 gr.SetPointError(bin,0, sigErr);
404 }
405 */
406 }
407 return gr;
408}
409
410
411//------------------------------------------------------------------------------
412
413
414// type 1 : object, 2 : track, 3 : tower
415
416void addGraph(TMultiGraph *mg, TGraphErrors *gr, TLegend *leg, int type)
417{
418
419 gr->SetLineWidth(3);
420
421 switch ( type )
422 {
423 case 1:
424 gr->SetLineColor(objColor);
425 gr->SetLineStyle(objStyle);
426 std::cout << "Adding " << gr->GetName() << std::endl;
427 mg->Add(gr);
428 leg->AddEntry(gr,"Reco","l");
429 break;
430
431 case 2:
432 gr->SetLineColor(trackColor);
433 gr->SetLineStyle(trackStyle);
434 mg->Add(gr);
435 leg->AddEntry(gr,"Track","l");
436 break;
437
438 case 3:
439 gr->SetLineColor(towerColor);
440 gr->SetLineStyle(towerStyle);
441 mg->Add(gr);
442 leg->AddEntry(gr,"Tower","l");
443 break;
444
445 default:
446 std::cout << "wrong type, possibles choices are Object, Track and Tower" << std::endl;
447 break;
448 }
449}
450
451void addHist(TH1D *h, TLegend *leg, int type)
452{
453 h->SetLineWidth(3);
454
455 switch ( type )
456 {
457 case 1:
458 h->SetLineColor(objColor);
459 h->SetLineStyle(objStyle);
460 leg->AddEntry(h,"Reco","l");
461 break;
462
463 case 2:
464 h->SetLineColor(trackColor);
465 h->SetLineStyle(trackStyle);
466 leg->AddEntry(h,"Track","l");
467 break;
468
469 case 3:
470 h->SetLineColor(towerColor);
471 h->SetLineStyle(towerStyle);
472 leg->AddEntry(h,"Tower","l");
473 break;
474
475 default:
476 std::cout << "wrong type, possibles choices are Object, Track and Tower" << std::endl;
477 break;
478 }
479}
480
481void DrawAxis(TMultiGraph *mg, TLegend *leg, double max)
482{
483 mg->SetMinimum(0.);
484 mg->SetMaximum(max);
485 mg->GetXaxis()->SetLimits(ptrangemin,ptrangemax);
486 mg->GetYaxis()->SetTitle("#sigma (E) / E_{gen}");
487 mg->GetXaxis()->SetTitle("E_{gen}");
488 mg->GetYaxis()->SetTitleSize(0.07);
489 mg->GetXaxis()->SetTitleSize(0.07);
490 mg->GetYaxis()->SetLabelSize(0.07);
491 mg->GetXaxis()->SetLabelSize(0.07);
492
493 leg->SetLineStyle(0);
494 leg->SetFillStyle(0);
495 leg->SetLineWidth(0);
496 leg->SetLineColor(0);
497
498 gStyle->SetOptTitle(0);
499 gPad->SetLogx();
500 gPad->SetBottomMargin(0.2);
501 gPad->SetLeftMargin(0.2);
502 gPad->Modified();
503 gPad->Update();
504
505}
506
507void DrawAxis(TH1D *h, TLegend *leg, int type)
508{
509
510 h->GetYaxis()->SetRangeUser(0,1.2);
511 if (type == 0) h->GetXaxis()->SetTitle("E_{gen}");
512 else h->GetXaxis()->SetTitle("#eta");
513 h->GetYaxis()->SetTitle("#epsilon");
514 h->GetYaxis()->SetTitleSize(0.07);
515 h->GetXaxis()->SetTitleSize(0.07);
516 h->GetYaxis()->SetLabelSize(0.07);
517 h->GetXaxis()->SetLabelSize(0.07);
518 leg->SetLineStyle(0);
519 leg->SetFillStyle(0);
520 leg->SetLineWidth(0);
521 leg->SetLineColor(0);
522
523 gStyle->SetOptTitle(0);
524 gStyle->SetOptStat(0);
525 gPad->SetBottomMargin(0.2);
526 gPad->SetLeftMargin(0.2);
527
528 gPad->Modified();
529 gPad->Update();
530
531}
532
533
534void Validation(const char *inputFile, const char *outputFile)
535{
536 //gSystem->Load("libDelphes");
537
538 std::cout << "input file : " << inputFile << " " << " , output file : " << outputFile << std::endl;
539
540 TChain *chain = new TChain("Delphes");
541 chain->Add(inputFile);
542
543 ExRootTreeReader *treeReader = new ExRootTreeReader(chain);
544 TClonesArray *branchParticle = treeReader->UseBranch("Particle");
545 TClonesArray *branchElectron = treeReader->UseBranch("Electron");
546 TClonesArray *branchMuon = treeReader->UseBranch("Muon");
547 TClonesArray *branchPhoton = treeReader->UseBranch("Photon");
548 TClonesArray *branchTrack = treeReader->UseBranch("Track");
549 TClonesArray *branchTower = treeReader->UseBranch("Tower");
550 TClonesArray *branchGenJet = treeReader->UseBranch("GenJet");
551 TClonesArray *branchPFJet = treeReader->UseBranch("Jet");
552 TClonesArray *branchCaloJet = treeReader->UseBranch("CaloJet");
553
554#ifdef ELECTRON
555
556 ///////////////
557 // Electrons //
558 ///////////////
559
560 // Reconstruction efficiency
561 TString elecs = "Electron";
562 int elID = 11;
563 std::pair<TH1D*,TH1D*> histos_el = GetEff<Electron>(branchElectron, branchParticle, "Electron", elID, treeReader);
564
565 // tracking reconstruction efficiency
566 std::pair <TH1D*,TH1D*> histos_eltrack = GetEff<Track>(branchTrack, branchParticle, "electronTrack", elID, treeReader);
567
568 // Tower reconstruction efficiency
569 std::pair <TH1D*,TH1D*> histos_eltower = GetEff<Tower>(branchTower, branchParticle, "electronTower", elID, treeReader);
570
571 // Electron Energy Resolution
572 std::vector<resolPlot> plots_el;
573 HistogramsCollection(&plots_el, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "electrons");
574 GetEres<Electron>( &plots_el, branchElectron, branchParticle, elID, treeReader);
575 TGraphErrors gr_el = EresGraph(&plots_el, true);
576 gr_el.SetName("Electron");
577
578 // Electron Track Energy Resolution
579 std::vector<resolPlot> plots_eltrack;
580 HistogramsCollection(&plots_eltrack, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "electronsTracks");
581 GetEres<Track>( &plots_eltrack, branchTrack, branchParticle, elID, treeReader);
582 TGraphErrors gr_eltrack = EresGraph(&plots_eltrack, true);
583 gr_eltrack.SetName("ElectronTracks");
584
585 // Electron Tower Energy Resolution
586 std::vector<resolPlot> plots_eltower;
587 HistogramsCollection(&plots_eltower, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "electronsTower");
588 GetEres<Tower>( &plots_eltower, branchTower, branchParticle, elID, treeReader);
589 TGraphErrors gr_eltower = EresGraph(&plots_eltower, true);
590 gr_eltower.SetName("ElectronTower");
591
592 // Canvases
593 TString elEff = "electronEff";
594 TCanvas *C_el1 = new TCanvas(elEff,elEff, 1000, 500);
595 C_el1->Divide(2);
596 C_el1->cd(1);
597 TLegend *leg_el1 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
598
599 gPad->SetLogx();
600 histos_eltrack.first->Draw();
601 addHist(histos_eltrack.first, leg_el1, 2);
602 histos_eltower.first->Draw("same");
603 addHist(histos_eltower.first, leg_el1, 3);
604 histos_el.first->Draw("same");
605 addHist(histos_el.first, leg_el1, 1);
606
607 DrawAxis(histos_eltrack.first, leg_el1, 0);
608 leg_el1->Draw();
609
610 C_el1->cd(2);
611 TLegend *leg_el2 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
612
613 histos_eltrack.second->Draw();
614 addHist(histos_eltrack.second, leg_el2, 2);
615 histos_eltower.second->Draw("same");
616 addHist(histos_eltower.second, leg_el2, 3);
617 histos_el.second->Draw("same");
618 addHist(histos_el.second, leg_el2, 1);
619
620 DrawAxis(histos_eltrack.second, leg_el2, 1);
621 delete(leg_el2);
622 C_el1->cd(0);
623
624 TString elRes = "electronERes";
625 TCanvas *C_el2 = new TCanvas(elRes,elRes, 1000, 500);
626 TMultiGraph *mg_el = new TMultiGraph(elRes,elRes);
627 TLegend *leg_el = new TLegend(0.52,0.7,0.9,0.9);
628
629 addGraph(mg_el, &gr_eltower, leg_el, 3);
630 addGraph(mg_el, &gr_eltrack, leg_el, 2);
631 addGraph(mg_el, &gr_el, leg_el, 1);
632
633 mg_el->Draw("ACX");
634 leg_el->Draw();
635
636 DrawAxis(mg_el, leg_el, 0.1);
637
638 C_el1->SaveAs(elEff+".eps");
639 C_el2->SaveAs(elRes+".eps");
640
641#endif
642
643 gDirectory->cd(0);
644
645#ifdef MUON
646
647 ///////////
648 // Muons //
649 ///////////
650
651 // Reconstruction efficiency
652 int muID = 13;
653 std::pair<TH1D*,TH1D*> histos_mu = GetEff<Muon>(branchMuon, branchParticle,"Muon", muID, treeReader);
654
655 // muon tracking reconstruction efficiency
656 std::pair <TH1D*,TH1D*> histos_mutrack = GetEff<Track>(branchTrack, branchParticle, "muonTrack", muID, treeReader);
657
658 // Muon Energy Resolution
659 std::vector<resolPlot> plots_mu;
660 HistogramsCollection(&plots_mu, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "muons");
661 GetEres<Muon>( &plots_mu, branchMuon, branchParticle, muID, treeReader);
662 TGraphErrors gr_mu = EresGraph(&plots_mu, true);
663 gr_mu.SetName("Muon");
664
665 // Muon Track Energy Resolution
666 std::vector<resolPlot> plots_mutrack;
667 HistogramsCollection(&plots_mutrack, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "muonsTracks");
668 GetEres<Track>( &plots_mutrack, branchTrack, branchParticle, muID, treeReader);
669 TGraphErrors gr_mutrack = EresGraph(&plots_mutrack, true);
670 gr_mutrack.SetName("MuonTracks");
671
672 // Canvas
673
674 TString muEff = "muonEff";
675 TCanvas *C_mu1 = new TCanvas(muEff,muEff, 1000, 500);
676 C_mu1->Divide(2);
677 C_mu1->cd(1);
678 TLegend *leg_mu1 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
679
680 gPad->SetLogx();
681 histos_mutrack.first->Draw();
682 addHist(histos_mutrack.first, leg_mu1, 2);
683 histos_mu.first->Draw("same");
684 addHist(histos_mu.first, leg_mu1, 1);
685
686 DrawAxis(histos_mutrack.first, leg_mu1, 0);
687 leg_mu1->Draw();
688
689 C_mu1->cd(2);
690 TLegend *leg_mu2 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
691
692 histos_mutrack.second->Draw();
693 addHist(histos_mutrack.second, leg_mu2, 2);
694 histos_mu.second->Draw("same");
695 addHist(histos_mu.second, leg_mu2, 1);
696
697 DrawAxis(histos_mutrack.second, leg_mu2, 1);
698
699 TString muRes = "muonERes";
700 TCanvas *C_mu = new TCanvas(muRes,muRes, 1000, 500);
701 TMultiGraph *mg_mu = new TMultiGraph(muRes,muRes);
702 TLegend *leg_mu = new TLegend(0.52,0.7,0.9,0.9);
703
704 addGraph(mg_mu, &gr_mutrack, leg_mu, 2);
705 addGraph(mg_mu, &gr_mu, leg_mu, 1);
706
707 mg_mu->Draw("ACX");
708 leg_mu->Draw();
709
710 DrawAxis(mg_mu, leg_mu, 0.3);
711
712 C_mu1->SaveAs(muEff+".eps");
713 C_mu->SaveAs(muRes+".eps");
714
715#endif
716
717 gDirectory->cd(0);
718
719#ifdef PHOTON
720
721 /////////////
722 // Photons //
723 /////////////
724
725 // Reconstruction efficiency
726 int phID = 22;
727 std::pair<TH1D*,TH1D*> histos_ph = GetEff<Electron>(branchPhoton, branchParticle, "Photon", phID, treeReader);
728 std::pair<TH1D*,TH1D*> histos_phtower = GetEff<Electron>(branchTower, branchParticle, "Photon", phID, treeReader);
729
730 // Photon Energy Resolution
731 std::vector<resolPlot> plots_ph;
732 HistogramsCollection(&plots_ph, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "photons");
733 GetEres<Photon>( &plots_ph, branchPhoton, branchParticle, phID, treeReader);
734 TGraphErrors gr_ph = EresGraph(&plots_ph, true);
735 gr_ph.SetName("Photon");
736
737 // Photon Tower Energy Resolution
738 std::vector<resolPlot> plots_phtower;
739 HistogramsCollection(&plots_phtower, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "photonsTower");
740 GetEres<Tower>( &plots_phtower, branchTower, branchParticle, phID, treeReader);
741 TGraphErrors gr_phtower = EresGraph(&plots_phtower, true);
742 gr_phtower.SetName("PhotonTower");
743
744 // Canvas
745
746 TString phEff = "photonEff";
747 TCanvas *C_ph1 = new TCanvas(phEff,phEff, 1000, 500);
748 C_ph1->Divide(2);
749 C_ph1->cd(1);
750 TLegend *leg_ph1 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
751
752 gPad->SetLogx();
753 histos_phtower.first->Draw();
754 addHist(histos_phtower.first, leg_ph1, 3);
755 histos_ph.first->Draw("same");
756 addHist(histos_ph.first, leg_ph1, 1);
757
758 DrawAxis(histos_phtower.first, leg_ph1, 0);
759 leg_ph1->Draw();
760
761 C_ph1->cd(2);
762 TLegend *leg_ph2 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
763
764 histos_phtower.second->Draw("same");
765 addHist(histos_phtower.second, leg_ph2, 3);
766 histos_ph.second->Draw("same");
767 addHist(histos_ph.second, leg_ph2, 1);
768
769 DrawAxis(histos_phtower.second, leg_ph2, 1);
770
771 C_ph1->SaveAs(phEff+".eps");
772
773 TString phRes = "phERes";
774 TCanvas *C_ph = new TCanvas(phRes,phRes, 1000, 500);
775 TMultiGraph *mg_ph = new TMultiGraph(phRes,phRes);
776 TLegend *leg_ph = new TLegend(0.52,0.7,0.9,0.9);
777
778 addGraph(mg_ph, &gr_phtower, leg_ph, 3);
779 addGraph(mg_ph, &gr_ph, leg_ph, 1);
780
781 mg_ph->Draw("ACX");
782 leg_ph->Draw();
783
784 DrawAxis(mg_ph, leg_ph, 0.3);
785
786 C_ph->SaveAs(phRes+".eps");
787
788#endif
789
790 gDirectory->cd(0);
791
792#ifdef JET
793
794 //////////
795 // Jets //
796 //////////
797
798 // BJets Reconstruction efficiency
799 int bID = 5;
800 std::pair<TH1D*,TH1D*> histos_btag = GetEff<Jet>(branchPFJet, branchParticle,"BTag", bID, treeReader);
801
802 // TauJets Reconstruction efficiency
803 int tauID = 15;
804 std::pair<TH1D*,TH1D*> histos_tautag = GetEff<Jet>(branchPFJet, branchParticle,"TauTag", tauID, treeReader);
805
806 // PFJets Energy Resolution
807 std::vector<resolPlot> plots_pfjets;
808 HistogramsCollection(&plots_pfjets, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "PFJet");
809 GetJetsEres( &plots_pfjets, branchPFJet, branchGenJet, treeReader);
810 TGraphErrors gr_pfjets = EresGraph(&plots_pfjets, true);
811 gr_pfjets.SetName("pfJet");
812
813
814 // PFJets Energy Resolution
815 std::vector<resolPlot> plots_calojets;
816 HistogramsCollection(&plots_calojets, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "CaloJet");
817 GetJetsEres( &plots_calojets, branchCaloJet, branchGenJet, treeReader);
818 TGraphErrors gr_calojets = EresGraph(&plots_calojets, true);
819 gr_calojets.SetName("caloJet");
820
821 // Canvas
822
823 TString btagEff = "btagEff";
824 TCanvas *C_btag1 = new TCanvas(btagEff,btagEff, 1000, 500);
825 C_btag1->Divide(2);
826 C_btag1->cd(1);
827 TLegend *leg_btag1 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
828
829 gPad->SetLogx();
830 histos_btag.first->Draw();
831 addHist(histos_btag.first, leg_btag1, 1);
832
833 DrawAxis(histos_btag.first, leg_btag1, 0);
834 leg_btag1->Draw();
835
836 C_btag1->cd(2);
837 TLegend *leg_btag2 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
838
839 histos_btag.second->Draw();
840 addHist(histos_btag.second, leg_btag2, 1);
841
842 DrawAxis(histos_btag.second, leg_btag2, 1);
843 delete(leg_btag2);
844 C_btag1->cd(0);
845
846 TString tautagEff = "tautagEff";
847 TCanvas *C_tautag1 = new TCanvas(tautagEff,tautagEff, 1000, 500);
848 C_tautag1->Divide(2);
849 C_tautag1->cd(1);
850 TLegend *leg_tautag1 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
851
852 gPad->SetLogx();
853 histos_tautag.first->Draw();
854 addHist(histos_tautag.first, leg_tautag1, 1);
855
856 DrawAxis(histos_tautag.first, leg_tautag1, 0);
857 leg_tautag1->Draw();
858
859 C_tautag1->cd(2);
860 TLegend *leg_tautag2 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
861
862 histos_tautag.second->Draw();
863 addHist(histos_tautag.second, leg_tautag2, 1);
864
865 DrawAxis(histos_tautag.second, leg_tautag2, 1);
866 delete(leg_tautag2);
867 C_tautag1->cd(0);
868
869 TString jetRes = "jetERes";
870 TCanvas *C_jet = new TCanvas(jetRes,jetRes, 1000, 500);
871 TMultiGraph *mg_jet = new TMultiGraph(jetRes,jetRes);
872 TLegend *leg_jet = new TLegend(0.52,0.7,0.9,0.9);
873
874 addGraph(mg_jet, &gr_calojets, leg_jet, 3);
875 addGraph(mg_jet, &gr_pfjets, leg_jet, 1);
876
877 mg_jet->Draw("ACX");
878 leg_jet->Draw();
879
880 DrawAxis(mg_jet, leg_jet, 0.25);
881
882 C_btag1->SaveAs(btagEff+".eps");
883 C_jet->SaveAs(jetRes+".eps");
884
885#endif
886
887 /*
888 // CaloJets Energy Resolution
889 std::vector<resolPlot> plots_calojets;
890 HistogramsCollection(&plots_calojets, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "caloJet");
891 GetJetsEres( &plots_calojets, branchCaloJet, branchGenJet, treeReader);
892 TGraphErrors gr_calojets = EresGraph(&plots_calojets, true);
893 gr_calojets.SetName("caloJet");
894 */
895
896 TFile *fout = new TFile(outputFile,"recreate");
897
898#ifdef ELECTRON
899 for (int bin = 0; bin < Nbins; bin++)
900 {
901 plots_el.at(bin).cenResolHist->Write();
902 plots_eltrack.at(bin).cenResolHist->Write();
903 plots_eltower.at(bin).cenResolHist->Write();
904 }
905 // gr.Write();
906 histos_el.first->Write();
907 //histos_el.second->Write();
908 histos_eltrack.first->Write();
909 //histos_eltrack.second->Write();
910 histos_eltower.first->Write();
911 C_el1->Write();
912 C_el2->Write();
913#endif
914
915#ifdef MUON
916 histos_mu.first->Write();
917 histos_mu.second->Write();
918 histos_mutrack.first->Write();
919 histos_mutrack.second->Write();
920 C_mu1->Write();
921 C_mu->Write();
922#endif
923
924#ifdef PHOTON
925 histos_ph.first->Write();
926 histos_ph.second->Write();
927 C_ph1->Write();
928 C_ph->Write();
929#endif
930
931#ifdef JET
932 for (int bin = 0; bin < Nbins; bin++)
933 {
934 plots_pfjets.at(bin).cenResolHist->Write();
935 plots_calojets.at(bin).cenResolHist->Write();
936 }
937 histos_btag.first->Write();
938 histos_btag.second->Write();
939 histos_tautag.first->Write();
940 histos_tautag.second->Write();
941 C_btag1->Write();
942 C_tautag1->Write();
943 C_jet->Write();
944#endif
945
946 fout->Write();
947
948 cout << "** Exiting..." << endl;
949
950 //delete plots;
951 //delete result;
952 delete treeReader;
953 delete chain;
954}
955
956//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.