Fork me on GitHub

source: git/validation/DelphesValidation.cpp@ 505a779

ImprovedOutputFile Timing
Last change on this file since 505a779 was 341014c, checked in by Pavel Demin <pavel-demin@…>, 6 years ago

apply .clang-format to all .h, .cc and .cpp files
  • Property mode set to 100644
File size: 107.5 KB
Line 
1/*
2 * Delphes: a framework for fast simulation of a generic collider experiment
3 * Copyright (C) 2012-2014 Universite catholique de Louvain (UCL), Belgium
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19#include <iostream>
20#include <typeinfo>
21#include <utility>
22#include <vector>
23
24#include "TApplication.h"
25#include "TROOT.h"
26#include "TSystem.h"
27
28#include "TString.h"
29
30#include "TCanvas.h"
31#include "TClonesArray.h"
32#include "TGraph.h"
33#include "TGraphErrors.h"
34#include "TH1.h"
35#include "TH2.h"
36#include "THStack.h"
37#include "TLegend.h"
38#include "TLorentzVector.h"
39#include "TMath.h"
40#include "TMultiGraph.h"
41#include "TPaveText.h"
42#include "TStyle.h"
43
44#include "classes/DelphesClasses.h"
45
46#include "ExRootAnalysis/ExRootResult.h"
47#include "ExRootAnalysis/ExRootTreeBranch.h"
48#include "ExRootAnalysis/ExRootTreeReader.h"
49#include "ExRootAnalysis/ExRootTreeWriter.h"
50#include "ExRootAnalysis/ExRootUtilities.h"
51
52using namespace std;
53
54//------------------------------------------------------------------------------
55
56static const int Nbins = 50;
57
58int objStyle = 1;
59int trackStyle = 7;
60int towerStyle = 3;
61
62Color_t objColor = kBlack;
63Color_t trackColor = kBlack;
64Color_t towerColor = kBlack;
65
66double effLegXmin = 0.22;
67double effLegXmax = 0.7;
68double effLegYmin = 0.22;
69double effLegYmax = 0.5;
70
71double resLegXmin = 0.62;
72double resLegXmax = 0.9;
73double resLegYmin = 0.52;
74double resLegYmax = 0.85;
75
76double topLeftLegXmin = 0.22;
77double topLeftLegXmax = 0.7;
78double topLeftLegYmin = 0.52;
79double topLeftLegYmax = 0.85;
80
81unsigned int k;
82
83struct resolPlot
84{
85 TH1 *resolHist;
86 double ptmin;
87 double ptmax;
88 double etamin;
89 double etamax;
90 double xmin;
91 double xmax;
92 TString obj;
93
94 resolPlot();
95 resolPlot(double ptdown, double ptup, TString object);
96 resolPlot(double etadown, double etaup, double ptdown, double ptup, TString object);
97 void set(double ptdown, double ptup, TString object, double xmin = 0, double xmax = 2);
98 void set(double etadown, double etaup, double ptdown, double ptup, TString object, double xmin = 0, double xmax = 2);
99 void print() { std::cout << ptmin << std::endl; }
100};
101
102resolPlot::resolPlot()
103{
104}
105
106resolPlot::resolPlot(double ptdown, double ptup, TString object)
107{
108 this->set(ptdown, ptup, object);
109}
110
111resolPlot::resolPlot(double etadown, double etaup, double ptdown, double ptup, TString object)
112{
113 this->set(etadown, etaup, ptdown, ptup, object);
114}
115
116void resolPlot::set(double ptdown, double ptup, TString object, double xmin, double xmax)
117{
118 ptmin = ptdown;
119 ptmax = ptup;
120 obj = object;
121
122 resolHist = new TH1D(obj + "_delta_pt_" + Form("%4.2f", ptmin) + "_" + Form("%4.2f", ptmax), obj + "_delta_pt_" + Form("%4.2f", ptmin) + "_" + Form("%4.2f", ptmax), 1000, xmin, xmax);
123}
124
125void resolPlot::set(double etadown, double etaup, double ptdown, double ptup, TString object, double xmin, double xmax)
126{
127 etamin = etadown;
128 etamax = etaup;
129 ptmin = ptdown;
130 ptmax = ptup;
131 obj = object;
132
133 resolHist = new TH1D(obj + "_delta_pt_" + Form("%4.2f", ptmin) + "_" + Form("%4.2f", ptmax) + "_" + Form("%4.2f", etamin) + "_" + Form("%4.2f", etamax), obj + "_delta_pt_" + Form("%4.2f", ptmin) + "_" + Form("%4.2f", ptmax) + "_" + Form("%4.2f", etamin) + "_" + Form("%4.2f", etamax), 1000, xmin, xmax);
134}
135
136void HistogramsCollection(std::vector<resolPlot> *histos, double ptmin, double ptmax, TString obj, double xmin = 0, double xmax = 2)
137{
138 double width;
139 double ptdown;
140 double ptup;
141 resolPlot ptemp;
142
143 for(int i = 0; i < Nbins; i++)
144 {
145 width = (ptmax - ptmin) / Nbins;
146 ptdown = TMath::Power(10, ptmin + i * width);
147 ptup = TMath::Power(10, ptmin + (i + 1) * width);
148 ptemp.set(ptdown, ptup, obj, xmin, xmax);
149 histos->push_back(ptemp);
150 }
151}
152
153void HistogramsCollectionVsEta(std::vector<resolPlot> *histos, double etamin, double etamax, double ptmin, double ptmax, TString obj, double xmin = 0, double xmax = 2)
154{
155 resolPlot ptemp;
156 double width;
157 double etadown;
158 double etaup;
159
160 for(int i = 0; i < Nbins; i++)
161 {
162 width = (etamax - etamin) / Nbins;
163 etadown = etamin + i * width;
164 etaup = etamin + (i + 1) * width;
165
166 ptemp.set(etadown, etaup, ptmin, ptmax, obj, xmin, xmax);
167 histos->push_back(ptemp);
168 }
169}
170
171//------------------------------------------------------------------------------
172
173class ExRootResult;
174class ExRootTreeReader;
175
176//------------------------------------------------------------------------------
177
178void BinLogX(TH1 *h)
179{
180 TAxis *axis = h->GetXaxis();
181 int bins = axis->GetNbins();
182
183 Axis_t from = axis->GetXmin();
184 Axis_t to = axis->GetXmax();
185 Axis_t width = (to - from) / bins;
186 Axis_t *new_bins = new Axis_t[bins + 1];
187
188 for(int i = 0; i <= bins; i++)
189 {
190 new_bins[i] = TMath::Power(10, from + i * width);
191 }
192 axis->Set(bins, new_bins);
193 delete new_bins;
194}
195
196//------------------------------------------------------------------------------
197
198template <typename T>
199TH1D *GetEffPt(TClonesArray *branchReco, TClonesArray *branchParticle, TString name, int pdgID, double ptmin, double ptmax, double etamin, double etamax, ExRootTreeReader *treeReader)
200{
201
202 cout << "** Computing Efficiency of reconstructing " << branchReco->GetName() << " induced by " << branchParticle->GetName() << " with PID " << pdgID << endl;
203
204 Long64_t allEntries = treeReader->GetEntries();
205
206 GenParticle *particle;
207 T *recoObj;
208
209 TLorentzVector recoMomentum, genMomentum, bestRecoMomentum;
210
211 Float_t deltaR;
212 Float_t pt, eta;
213 Long64_t entry;
214
215 Int_t i, j;
216
217 TH1D *histGenPt = new TH1D(name + " gen spectra Pt", name + " gen spectra cen", Nbins, TMath::Log10(ptmin), TMath::Log10(ptmax));
218 TH1D *histRecoPt = new TH1D(name + " reco spectra Pt", name + " reco spectra cen", Nbins, TMath::Log10(ptmin), TMath::Log10(ptmax));
219
220 histGenPt->SetDirectory(0);
221 histRecoPt->SetDirectory(0);
222
223 BinLogX(histGenPt);
224 BinLogX(histRecoPt);
225
226 // Loop over all events
227 for(entry = 0; entry < allEntries; ++entry)
228 {
229 // Load selected branches with data from specified event
230 treeReader->ReadEntry(entry);
231
232 // Loop over all generated particle in event
233 for(i = 0; i < branchParticle->GetEntriesFast(); ++i)
234 {
235
236 particle = (GenParticle *)branchParticle->At(i);
237 genMomentum = particle->P4();
238
239 deltaR = 999;
240
241 pt = genMomentum.Pt();
242 eta = TMath::Abs(genMomentum.Eta());
243
244 if(eta > etamax || eta < etamin) continue;
245
246 if(particle->PID == pdgID && genMomentum.Pt() > ptmin && genMomentum.Pt() < ptmax)
247 //if (TMath::Abs(particle->PID) == pdgID && (particle->Status>20 && particle->Status <30) && genMomentum.Pt() > ptmin && genMomentum.Pt() < ptmax )
248 {
249 // Loop over all reco object in event
250 for(j = 0; j < branchReco->GetEntriesFast(); ++j)
251 {
252 recoObj = (T *)branchReco->At(j);
253 recoMomentum = recoObj->P4();
254 //if(Momentum.Px() == 0 && genMomentum.Py() == 0) continue;
255
256 // take the closest parton candidate
257 if(TMath::Abs(pdgID) == 5)
258 {
259 Jet *jet = (Jet *)recoObj;
260 if(!(jet->BTag & (1 << 0))) continue;
261
262 //if(jet->BTag != ) continue;
263 }
264
265 if(TMath::Abs(pdgID) == 4)
266 {
267 Jet *jet = (Jet *)recoObj;
268 if(!(jet->BTag & (1 << 0))) continue;
269 }
270
271 if(TMath::Abs(pdgID) == 1)
272 {
273 Jet *jet = (Jet *)recoObj;
274 if(!(jet->BTag & (1 << 0))) continue;
275 }
276
277 if(TMath::Abs(pdgID) == 15)
278 {
279 Jet *jet = (Jet *)recoObj;
280 if(jet->TauTag != 1) continue;
281 }
282
283 if(genMomentum.DeltaR(recoMomentum) < deltaR)
284 {
285 deltaR = genMomentum.DeltaR(recoMomentum);
286 bestRecoMomentum = recoMomentum;
287 }
288 }
289 histGenPt->Fill(pt);
290 if(deltaR < 0.3 && bestRecoMomentum.Pt() > 0.20 * pt)
291 {
292 histRecoPt->Fill(pt);
293 }
294 }
295 }
296 }
297
298 histRecoPt->Sumw2();
299 histGenPt->Sumw2();
300
301 histRecoPt->Divide(histGenPt);
302 histRecoPt->Scale(100.);
303
304 return histRecoPt;
305}
306
307template <typename T>
308TH1D *GetEffEta(TClonesArray *branchReco, TClonesArray *branchParticle, TString name, int pdgID, double ptmin, double ptmax, double etamin, double etamax, ExRootTreeReader *treeReader)
309{
310
311 cout << "** Computing Efficiency of reconstructing " << branchReco->GetName() << " induced by " << branchParticle->GetName() << " with PID " << pdgID << endl;
312
313 Long64_t allEntries = treeReader->GetEntries();
314
315 GenParticle *particle;
316 T *recoObj;
317
318 TLorentzVector recoMomentum, genMomentum, bestRecoMomentum;
319
320 Float_t deltaR;
321 Float_t pt, eta;
322 Long64_t entry;
323
324 Int_t i, j;
325
326 TH1D *histGenEta = new TH1D(name + " gen spectra Eta", name + " gen spectra", Nbins, etamin, etamax);
327 TH1D *histRecoEta = new TH1D(name + " reco spectra Eta", name + " reco spectra", Nbins, etamin, etamax);
328
329 histGenEta->SetDirectory(0);
330 histRecoEta->SetDirectory(0);
331
332 // Loop over all events
333 for(entry = 0; entry < allEntries; ++entry)
334 {
335 // Load selected branches with data from specified event
336 treeReader->ReadEntry(entry);
337
338 // Loop over all generated particle in event
339 for(i = 0; i < branchParticle->GetEntriesFast(); ++i)
340 {
341
342 particle = (GenParticle *)branchParticle->At(i);
343 genMomentum = particle->P4();
344
345 deltaR = 999;
346
347 pt = genMomentum.Pt();
348 eta = genMomentum.Eta();
349
350 if(pt > ptmax || pt < ptmin) continue;
351
352 if(particle->PID == pdgID && genMomentum.Pt() > ptmin && genMomentum.Pt() < ptmax)
353 //if (TMath::Abs(particle->PID) == pdgID && (particle->Status>20 && particle->Status <30) && genMomentum.Pt() > ptmin && genMomentum.Pt() < ptmax )
354 {
355 // Loop over all reco object in event
356 for(j = 0; j < branchReco->GetEntriesFast(); ++j)
357 {
358 recoObj = (T *)branchReco->At(j);
359 recoMomentum = recoObj->P4();
360 // this is simply to avoid warnings from initial state particle
361 // having infite rapidity ...
362 //if(Momentum.Px() == 0 && genMomentum.Py() == 0) continue;
363
364 // take the closest parton candidate
365 if(TMath::Abs(pdgID) == 5)
366 {
367 Jet *jet = (Jet *)recoObj;
368 if(!(jet->BTag & (1 << 0))) continue;
369 }
370
371 if(TMath::Abs(pdgID) == 4)
372 {
373 Jet *jet = (Jet *)recoObj;
374 if(!(jet->BTag & (1 << 0))) continue;
375 }
376
377 if(TMath::Abs(pdgID) == 1)
378 {
379 Jet *jet = (Jet *)recoObj;
380 if(!(jet->BTag & (1 << 0))) continue;
381 }
382
383 if(TMath::Abs(pdgID) == 15)
384 {
385 Jet *jet = (Jet *)recoObj;
386 if(jet->TauTag != 1) continue;
387 }
388 if(genMomentum.DeltaR(recoMomentum) < deltaR)
389 {
390 deltaR = genMomentum.DeltaR(recoMomentum);
391 bestRecoMomentum = recoMomentum;
392 }
393 }
394
395 histGenEta->Fill(eta);
396 if(deltaR < 0.3)
397 {
398 histRecoEta->Fill(eta);
399 }
400 }
401 }
402 }
403
404 histRecoEta->Sumw2();
405 histGenEta->Sumw2();
406
407 histRecoEta->Divide(histGenEta);
408 histRecoEta->Scale(100.);
409
410 return histRecoEta;
411}
412
413//------------------------------------------------------------------------------
414
415template <typename T>
416TH1D *GetJetEffPt(TClonesArray *branchJet, TString name, int pdgID, double ptmin, double ptmax, double etamin, double etamax, ExRootTreeReader *treeReader)
417{
418
419 cout << "** Computing Efficiency of reconstructing " << branchJet->GetName() << " with PID " << pdgID << endl;
420
421 Long64_t allEntries = treeReader->GetEntries();
422
423 T *recoObj;
424
425 TLorentzVector recoMomentum, genMomentum, bestRecoMomentum;
426
427 Float_t pt, eta;
428 Long64_t entry;
429
430 Int_t j;
431
432 TH1D *histGenPt = new TH1D(name + " gen spectra Pt", name + " gen spectra cen", Nbins, TMath::Log10(ptmin), TMath::Log10(ptmax));
433 TH1D *histRecoPt = new TH1D(name + " reco spectra Pt", name + " reco spectra cen", Nbins, TMath::Log10(ptmin), TMath::Log10(ptmax));
434
435 histGenPt->SetDirectory(0);
436 histRecoPt->SetDirectory(0);
437
438 BinLogX(histGenPt);
439 BinLogX(histRecoPt);
440
441 // Loop over all events
442 for(entry = 0; entry < allEntries; ++entry)
443 {
444 // Load selected branches with data from specified event
445 treeReader->ReadEntry(entry);
446
447 // Loop over all reco object in event
448 for(j = 0; j < branchJet->GetEntriesFast(); ++j)
449 {
450 recoObj = (T *)branchJet->At(j);
451 recoMomentum = recoObj->P4();
452 pt = recoMomentum.Pt();
453 eta = TMath::Abs(recoMomentum.Eta());
454 Jet *jet = (Jet *)recoObj;
455
456 if(eta > etamax || eta < etamin) continue;
457 if(pt < ptmin || pt > ptmax) continue;
458
459 Int_t flavor = jet->Flavor;
460 if(flavor == 21) flavor = 0;
461
462 if(TMath::Abs(pdgID) == 1)
463 {
464
465 if(flavor < 4)
466 {
467 histGenPt->Fill(pt);
468 if(jet->BTag & (1 << 0)) histRecoPt->Fill(pt);
469 }
470 }
471 if(TMath::Abs(pdgID) == 4)
472 {
473 if(flavor == 4)
474 {
475 histGenPt->Fill(pt);
476 if(jet->BTag & (1 << 0)) histRecoPt->Fill(pt);
477 }
478 }
479 if(TMath::Abs(pdgID) == 5)
480 {
481 if(flavor == 5)
482 {
483 histGenPt->Fill(pt);
484 if(jet->BTag & (1 << 0)) histRecoPt->Fill(pt);
485 }
486 }
487 }
488 }
489
490 histRecoPt->Sumw2();
491 histGenPt->Sumw2();
492
493 histRecoPt->Divide(histGenPt);
494 histRecoPt->Scale(100.);
495
496 return histRecoPt;
497}
498
499// ------------------------------------------------------------------------------------------------------------------------------------------------------
500
501template <typename T>
502TH1D *GetJetEffEta(TClonesArray *branchJet, TString name, int pdgID, double ptmin, double ptmax, double etamin, double etamax, ExRootTreeReader *treeReader)
503{
504
505 cout << "** Computing Efficiency of reconstructing " << branchJet->GetName() << " with PID " << pdgID << endl;
506
507 Long64_t allEntries = treeReader->GetEntries();
508
509 T *recoObj;
510
511 TLorentzVector recoMomentum, genMomentum, bestRecoMomentum;
512
513 Float_t pt, eta;
514 Long64_t entry;
515
516 Int_t j;
517
518 TH1D *histGenEta = new TH1D(name + " gen spectra Eta", name + " gen spectra", Nbins, etamin, etamax);
519 TH1D *histRecoEta = new TH1D(name + " reco spectra Eta", name + " reco spectra", Nbins, etamin, etamax);
520
521 histGenEta->SetDirectory(0);
522 histRecoEta->SetDirectory(0);
523 // Loop over all events
524 for(entry = 0; entry < allEntries; ++entry)
525 {
526 // Load selected branches with data from specified event
527 treeReader->ReadEntry(entry);
528
529 // Loop over all reco object in event
530 for(j = 0; j < branchJet->GetEntriesFast(); ++j)
531 {
532 recoObj = (T *)branchJet->At(j);
533 recoMomentum = recoObj->P4();
534 pt = recoMomentum.Pt();
535 eta = recoMomentum.Eta();
536 Jet *jet = (Jet *)recoObj;
537
538 if(eta > etamax || eta < etamin) continue;
539 if(pt < ptmin || pt > ptmax) continue;
540
541 Int_t flavor = jet->Flavor;
542 if(flavor == 21) flavor = 0;
543
544 if(TMath::Abs(pdgID) == 1)
545 {
546 if(flavor == 1 || flavor == 21)
547 {
548 histGenEta->Fill(eta);
549 if(jet->BTag & (1 << 0)) histRecoEta->Fill(eta);
550 }
551 }
552 if(TMath::Abs(pdgID) == 4)
553 {
554 if(flavor == 4)
555 {
556 histGenEta->Fill(eta);
557 if(jet->BTag & (1 << 0)) histRecoEta->Fill(eta);
558 }
559 }
560 if(TMath::Abs(pdgID) == 5)
561 {
562 if(flavor == 5)
563 {
564 histGenEta->Fill(eta);
565 if(jet->BTag & (1 << 0)) histRecoEta->Fill(eta);
566 }
567 }
568 }
569 }
570
571 histRecoEta->Sumw2();
572 histGenEta->Sumw2();
573
574 histRecoEta->Divide(histGenEta);
575 histRecoEta->Scale(100.);
576
577 return histRecoEta;
578}
579
580// -----------------------------------------------------------------------------------------------------------------------------------------------------
581
582template <typename T>
583TH1D *GetTauEffPt(TClonesArray *branchReco, TClonesArray *branchParticle, TString name, int pdgID, double ptmin, double ptmax, double etamin, double etamax, ExRootTreeReader *treeReader)
584{
585
586 cout << "** Computing Efficiency of reconstructing " << branchReco->GetName() << " induced by " << branchParticle->GetName() << " with PID " << pdgID << endl;
587
588 Long64_t allEntries = treeReader->GetEntries();
589
590 GenParticle *particle;
591 T *recoObj;
592
593 TLorentzVector recoMomentum, genMomentum, bestRecoMomentum;
594
595 Float_t deltaR;
596 Float_t pt, eta;
597 Long64_t entry;
598
599 Int_t i, j;
600
601 TH1D *histGenPt = new TH1D(name + " gen spectra Pt", name + " gen spectra cen", Nbins, TMath::Log10(ptmin), TMath::Log10(ptmax));
602 TH1D *histRecoPt = new TH1D(name + " reco spectra Pt", name + " reco spectra cen", Nbins, TMath::Log10(ptmin), TMath::Log10(ptmax));
603
604 histGenPt->SetDirectory(0);
605 histRecoPt->SetDirectory(0);
606
607 BinLogX(histGenPt);
608 BinLogX(histRecoPt);
609
610 // Loop over all events
611 for(entry = 0; entry < allEntries; ++entry)
612 {
613 // Load selected branches with data from specified event
614 treeReader->ReadEntry(entry);
615
616 // Loop over all generated particle in event
617 for(i = 0; i < branchParticle->GetEntriesFast(); ++i)
618 {
619
620 particle = (GenParticle *)branchParticle->At(i);
621 genMomentum = particle->P4();
622
623 deltaR = 999;
624
625 pt = genMomentum.Pt();
626 eta = TMath::Abs(genMomentum.Eta());
627
628 if(eta > etamax || eta < etamin) continue;
629
630 if(particle->PID == pdgID && genMomentum.Pt() > ptmin && genMomentum.Pt() < ptmax)
631 {
632 // Loop over all reco object in event
633 for(j = 0; j < branchReco->GetEntriesFast(); ++j)
634 {
635 recoObj = (T *)branchReco->At(j);
636 recoMomentum = recoObj->P4();
637 // this is simply to avoid warnings from initial state particle
638 // having infite rapidity ...
639 //if(Momentum.Px() == 0 && genMomentum.Py() == 0) continue;
640
641 if(TMath::Abs(pdgID) == 1)
642 {
643 Jet *jet = (Jet *)recoObj;
644 if(jet->TauTag != 1) continue;
645 }
646
647 if(TMath::Abs(pdgID) == 15)
648 {
649 Jet *jet = (Jet *)recoObj;
650 if(jet->TauTag != 1) continue;
651 }
652
653 if(genMomentum.DeltaR(recoMomentum) < deltaR)
654 {
655 deltaR = genMomentum.DeltaR(recoMomentum);
656 bestRecoMomentum = recoMomentum;
657 }
658 }
659
660 histGenPt->Fill(pt);
661 if(deltaR < 0.3)
662 {
663 histRecoPt->Fill(pt);
664 }
665 }
666 }
667 }
668
669 histRecoPt->Sumw2();
670 histGenPt->Sumw2();
671
672 histRecoPt->Divide(histGenPt);
673 histRecoPt->Scale(100.);
674 if(TMath::Abs(pdgID) == 15) histRecoPt->Scale(1 / 0.648);
675
676 return histRecoPt;
677}
678
679template <typename T>
680TH1D *GetTauEffEta(TClonesArray *branchReco, TClonesArray *branchParticle, TString name, int pdgID, double ptmin, double ptmax, double etamin, double etamax, ExRootTreeReader *treeReader)
681{
682
683 cout << "** Computing Efficiency of reconstructing " << branchReco->GetName() << " induced by " << branchParticle->GetName() << " with PID " << pdgID << endl;
684
685 Long64_t allEntries = treeReader->GetEntries();
686
687 GenParticle *particle;
688 T *recoObj;
689
690 TLorentzVector recoMomentum, genMomentum, bestRecoMomentum;
691
692 Float_t deltaR;
693 Float_t pt, eta;
694 Long64_t entry;
695
696 Int_t i, j;
697
698 TH1D *histGenEta = new TH1D(name + " gen spectra Eta", name + " gen spectra", Nbins, etamin, etamax);
699 TH1D *histRecoEta = new TH1D(name + " reco spectra Eta", name + " reco spectra", Nbins, etamin, etamax);
700
701 histGenEta->SetDirectory(0);
702 histRecoEta->SetDirectory(0);
703
704 // Loop over all events
705 for(entry = 0; entry < allEntries; ++entry)
706 {
707 // Load selected branches with data from specified event
708 treeReader->ReadEntry(entry);
709
710 // Loop over all generated particle in event
711 for(i = 0; i < branchParticle->GetEntriesFast(); ++i)
712 {
713
714 particle = (GenParticle *)branchParticle->At(i);
715 genMomentum = particle->P4();
716
717 deltaR = 999;
718
719 pt = genMomentum.Pt();
720 eta = genMomentum.Eta();
721
722 if(pt > ptmax || pt < ptmin) continue;
723
724 if(particle->PID == pdgID && genMomentum.Pt() > ptmin && genMomentum.Pt() < ptmax)
725 {
726 // Loop over all reco object in event
727 for(j = 0; j < branchReco->GetEntriesFast(); ++j)
728 {
729 recoObj = (T *)branchReco->At(j);
730 recoMomentum = recoObj->P4();
731 // this is simply to avoid warnings from initial state particle
732 // having infite rapidity ...
733 //if(Momentum.Px() == 0 && genMomentum.Py() == 0) continue;
734
735 if(TMath::Abs(pdgID) == 1)
736 {
737 Jet *jet = (Jet *)recoObj;
738 if(jet->TauTag != 1) continue;
739 }
740
741 if(TMath::Abs(pdgID) == 15)
742 {
743 Jet *jet = (Jet *)recoObj;
744 if(jet->TauTag != 1) continue;
745 }
746
747 if(genMomentum.DeltaR(recoMomentum) < deltaR)
748 {
749 deltaR = genMomentum.DeltaR(recoMomentum);
750 bestRecoMomentum = recoMomentum;
751 }
752 }
753
754 histGenEta->Fill(eta);
755 if(deltaR < 0.3)
756 {
757 histRecoEta->Fill(eta);
758 }
759 }
760 }
761 }
762
763 histRecoEta->Sumw2();
764 histGenEta->Sumw2();
765
766 histRecoEta->Divide(histGenEta);
767 histRecoEta->Scale(100.);
768 if(TMath::Abs(pdgID) == 15) histRecoEta->Scale(1 / 0.648);
769
770 return histRecoEta;
771}
772
773template <typename T>
774void GetPtres(std::vector<resolPlot> *histos, TClonesArray *branchReco, TClonesArray *branchParticle, int pdgID, Double_t etaMin, Double_t etaMax, ExRootTreeReader *treeReader)
775{
776 Long64_t allEntries = treeReader->GetEntries();
777
778 cout << "** Computing pt resolution of " << branchReco->GetName() << " induced by " << branchParticle->GetName() << " with PID " << pdgID << endl;
779
780 GenParticle *particle;
781 T *recoObj;
782
783 TLorentzVector recoMomentum, genMomentum, bestGenMomentum;
784
785 Float_t deltaR;
786 Float_t pt, eta;
787 Long64_t entry;
788
789 Int_t i, j, bin;
790
791 // Loop over all events
792 for(entry = 0; entry < allEntries; ++entry)
793 {
794 // Load selected branches with data from specified event
795 treeReader->ReadEntry(entry);
796
797 // Loop over all reconstructed jets in event
798 for(i = 0; i < branchReco->GetEntriesFast(); ++i)
799 {
800 recoObj = (T *)branchReco->At(i);
801 recoMomentum = recoObj->P4();
802
803 deltaR = 999;
804
805 // Loop over all hard partons in event
806 for(j = 0; j < branchParticle->GetEntriesFast(); ++j)
807 {
808 particle = (GenParticle *)branchParticle->At(j);
809 if(particle->PID == pdgID && particle->Status == 1)
810 {
811 genMomentum = particle->P4();
812
813 // this is simply to avoid warnings from initial state particle
814 // having infite rapidity ...
815 if(genMomentum.Px() == 0 && genMomentum.Py() == 0) continue;
816
817 // take the closest parton candidate
818 if(genMomentum.DeltaR(recoMomentum) < deltaR)
819 {
820 deltaR = genMomentum.DeltaR(recoMomentum);
821 bestGenMomentum = genMomentum;
822 }
823 }
824 }
825
826 if(deltaR < 0.3)
827 {
828 pt = bestGenMomentum.Pt();
829 eta = TMath::Abs(bestGenMomentum.Eta());
830
831 for(bin = 0; bin < Nbins; bin++)
832 {
833 if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta > etaMin && eta < etaMax)
834 {
835 histos->at(bin).resolHist->Fill(recoMomentum.Pt() / bestGenMomentum.Pt());
836 }
837 }
838 }
839 }
840 }
841}
842
843template <typename T>
844void GetEres(std::vector<resolPlot> *histos, TClonesArray *branchReco, TClonesArray *branchParticle, int pdgID, Double_t etaMin, Double_t etaMax, ExRootTreeReader *treeReader)
845{
846 Long64_t allEntries = treeReader->GetEntries();
847
848 cout << "** Computing e resolution of " << branchReco->GetName() << " induced by " << branchParticle->GetName() << " with PID " << pdgID << endl;
849
850 GenParticle *particle;
851 T *recoObj;
852
853 TLorentzVector recoMomentum, genMomentum, bestGenMomentum;
854
855 Float_t deltaR;
856 Float_t e, eta;
857 Long64_t entry;
858
859 Int_t i, j, bin;
860
861 // Loop over all events
862 for(entry = 0; entry < allEntries; ++entry)
863 {
864 // Load selected branches with data from specified event
865 treeReader->ReadEntry(entry);
866
867 // Loop over all reconstructed jets in event
868 for(i = 0; i < branchReco->GetEntriesFast(); ++i)
869 {
870 recoObj = (T *)branchReco->At(i);
871 recoMomentum = recoObj->P4();
872
873 deltaR = 999;
874
875 // Loop over all hard partons in event
876 for(j = 0; j < branchParticle->GetEntriesFast(); ++j)
877 {
878 particle = (GenParticle *)branchParticle->At(j);
879 if(particle->PID == pdgID && particle->Status == 1)
880 {
881 genMomentum = particle->P4();
882
883 // this is simply to avoid warnings from initial state particle
884 // having infite rapidity ...
885 if(genMomentum.Px() == 0 && genMomentum.Py() == 0) continue;
886
887 // take the closest parton candidate
888 if(genMomentum.DeltaR(recoMomentum) < deltaR)
889 {
890 deltaR = genMomentum.DeltaR(recoMomentum);
891 bestGenMomentum = genMomentum;
892 }
893 }
894 }
895
896 if(deltaR < 0.3)
897 {
898 e = bestGenMomentum.E();
899 eta = TMath::Abs(bestGenMomentum.Eta());
900
901 for(bin = 0; bin < Nbins; bin++)
902 {
903 if(e > histos->at(bin).ptmin && e < histos->at(bin).ptmax && eta > etaMin && eta < etaMax)
904 {
905 histos->at(bin).resolHist->Fill(recoMomentum.E() / bestGenMomentum.E());
906 }
907 }
908 }
909 }
910 }
911}
912
913template <typename T>
914void GetPtresVsEta(std::vector<resolPlot> *histos, TClonesArray *branchReco, TClonesArray *branchParticle, int pdgID, Double_t ptMin, Double_t ptMax, ExRootTreeReader *treeReader)
915{
916 Long64_t allEntries = treeReader->GetEntries();
917
918 cout << "** Computing pt resolution of " << branchReco->GetName() << " induced by " << branchParticle->GetName() << " with PID " << pdgID << endl;
919
920 GenParticle *particle;
921 T *recoObj;
922
923 TLorentzVector recoMomentum, genMomentum, bestGenMomentum;
924
925 Float_t deltaR;
926 Float_t pt, eta;
927 Long64_t entry;
928
929 Int_t i, j, bin;
930
931 // Loop over all events
932 for(entry = 0; entry < allEntries; ++entry)
933 {
934 // Load selected branches with data from specified event
935 treeReader->ReadEntry(entry);
936
937 // Loop over all reconstructed jets in event
938 for(i = 0; i < branchReco->GetEntriesFast(); ++i)
939 {
940 recoObj = (T *)branchReco->At(i);
941 recoMomentum = recoObj->P4();
942
943 deltaR = 999;
944
945 // Loop over all hard partons in event
946 for(j = 0; j < branchParticle->GetEntriesFast(); ++j)
947 {
948 particle = (GenParticle *)branchParticle->At(j);
949 if(particle->PID == pdgID && particle->Status == 1)
950 {
951 genMomentum = particle->P4();
952
953 // this is simply to avoid warnings from initial state particle
954 // having infite rapidity ...
955 if(genMomentum.Px() == 0 && genMomentum.Py() == 0) continue;
956
957 // take the closest parton candidate
958 if(genMomentum.DeltaR(recoMomentum) < deltaR)
959 {
960 deltaR = genMomentum.DeltaR(recoMomentum);
961 bestGenMomentum = genMomentum;
962 }
963 }
964 }
965
966 if(deltaR < 0.3)
967 {
968 pt = bestGenMomentum.Pt();
969 eta = bestGenMomentum.Eta();
970
971 for(bin = 0; bin < Nbins; bin++)
972 {
973 if(eta > histos->at(bin).etamin && eta < histos->at(bin).etamax && pt > ptMin && pt < ptMax)
974 {
975 histos->at(bin).resolHist->Fill(recoMomentum.Pt() / bestGenMomentum.Pt());
976 }
977 }
978 }
979 }
980 }
981}
982
983template <typename T>
984void GetEresVsEta(std::vector<resolPlot> *histos, TClonesArray *branchReco, TClonesArray *branchParticle, int pdgID, Double_t eMin, Double_t eMax, ExRootTreeReader *treeReader)
985{
986 Long64_t allEntries = treeReader->GetEntries();
987
988 cout << "** Computing E resolution of " << branchReco->GetName() << " induced by " << branchParticle->GetName() << " with PID " << pdgID << endl;
989
990 GenParticle *particle;
991 T *recoObj;
992
993 TLorentzVector recoMomentum, genMomentum, bestGenMomentum;
994
995 Float_t deltaR;
996 Float_t e, eta;
997 Long64_t entry;
998
999 Int_t i, j, bin;
1000
1001 // Loop over all events
1002 for(entry = 0; entry < allEntries; ++entry)
1003 {
1004 // Load selected branches with data from specified event
1005 treeReader->ReadEntry(entry);
1006
1007 // Loop over all reconstructed jets in event
1008 for(i = 0; i < branchReco->GetEntriesFast(); ++i)
1009 {
1010 recoObj = (T *)branchReco->At(i);
1011 recoMomentum = recoObj->P4();
1012
1013 deltaR = 999;
1014
1015 // Loop over all hard partons in event
1016 for(j = 0; j < branchParticle->GetEntriesFast(); ++j)
1017 {
1018 particle = (GenParticle *)branchParticle->At(j);
1019 if(particle->PID == pdgID && particle->Status == 1)
1020 {
1021 genMomentum = particle->P4();
1022
1023 // this is simply to avoid warnings from initial state particle
1024 // having infite rapidity ...
1025 if(genMomentum.Px() == 0 && genMomentum.Py() == 0) continue;
1026
1027 // take the closest parton candidate
1028 if(genMomentum.DeltaR(recoMomentum) < deltaR)
1029 {
1030 deltaR = genMomentum.DeltaR(recoMomentum);
1031 bestGenMomentum = genMomentum;
1032 }
1033 }
1034 }
1035
1036 if(deltaR < 0.3)
1037 {
1038 e = bestGenMomentum.E();
1039 eta = bestGenMomentum.Eta();
1040
1041 for(bin = 0; bin < Nbins; bin++)
1042 {
1043 if(eta > histos->at(bin).etamin && eta < histos->at(bin).etamax && e > eMin && e < eMax)
1044 {
1045 histos->at(bin).resolHist->Fill(recoMomentum.E() / bestGenMomentum.E());
1046 }
1047 }
1048 }
1049 }
1050 }
1051}
1052
1053void GetJetsEres(std::vector<resolPlot> *histos, TClonesArray *branchJet, TClonesArray *branchGenJet, ExRootTreeReader *treeReader, Double_t etaMin, Double_t etaMax)
1054{
1055
1056 Long64_t allEntries = treeReader->GetEntries();
1057
1058 cout << "** Computing resolution of " << branchJet->GetName() << " induced by " << branchGenJet->GetName() << endl;
1059
1060 Jet *jet, *genjet;
1061
1062 TLorentzVector jetMomentum, genJetMomentum, bestGenJetMomentum;
1063
1064 Float_t deltaR;
1065 Float_t pt, eta;
1066 Long64_t entry;
1067
1068 Int_t i, j, bin;
1069
1070 // Loop over all events
1071 for(entry = 0; entry < allEntries; ++entry)
1072 {
1073 // Load selected branches with data from specified event
1074 treeReader->ReadEntry(entry);
1075
1076 if(entry % 10000 == 0) cout << "Event number: " << entry << endl;
1077
1078 // Loop over all reconstructed jets in event
1079 for(i = 0; i < TMath::Min(2, branchJet->GetEntriesFast()); ++i) //branchJet->GetEntriesFast(); ++i)
1080 {
1081
1082 jet = (Jet *)branchJet->At(i);
1083 jetMomentum = jet->P4();
1084
1085 deltaR = 999;
1086
1087 // Loop over all hard partons in event
1088 for(j = 0; j < TMath::Min(2, branchGenJet->GetEntriesFast()); ++j)
1089 {
1090 genjet = (Jet *)branchGenJet->At(j);
1091
1092 genJetMomentum = genjet->P4();
1093
1094 // this is simply to avoid warnings from initial state particle
1095 // having infite rapidity ...
1096 if(genJetMomentum.Px() == 0 && genJetMomentum.Py() == 0) continue;
1097
1098 // take the closest parton candidate
1099 if(genJetMomentum.DeltaR(jetMomentum) < deltaR)
1100 {
1101 deltaR = genJetMomentum.DeltaR(jetMomentum);
1102 bestGenJetMomentum = genJetMomentum;
1103 }
1104 }
1105
1106 if(deltaR < 0.3)
1107 {
1108 pt = genJetMomentum.E();
1109 eta = genJetMomentum.Eta();
1110
1111 for(bin = 0; bin < Nbins; bin++)
1112 {
1113 if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta < etaMax && eta > etaMin)
1114 {
1115 histos->at(bin).resolHist->Fill(jetMomentum.E() / bestGenJetMomentum.E());
1116 }
1117 }
1118 }
1119 }
1120 }
1121}
1122
1123void GetJetsEresVsEta(std::vector<resolPlot> *histos, TClonesArray *branchJet, TClonesArray *branchGenJet, ExRootTreeReader *treeReader, Double_t eMin, Double_t eMax)
1124{
1125
1126 Long64_t allEntries = treeReader->GetEntries();
1127
1128 cout << "** Computing resolution of " << branchJet->GetName() << " induced by " << branchGenJet->GetName() << endl;
1129
1130 Jet *jet, *genjet;
1131
1132 TLorentzVector jetMomentum, genJetMomentum, bestGenJetMomentum;
1133
1134 Float_t deltaR;
1135 Float_t pt, eta;
1136 Long64_t entry;
1137
1138 Int_t i, j, bin;
1139
1140 // Loop over all events
1141 for(entry = 0; entry < allEntries; ++entry)
1142 {
1143 // Load selected branches with data from specified event
1144 treeReader->ReadEntry(entry);
1145
1146 if(entry % 10000 == 0) cout << "Event number: " << entry << endl;
1147
1148 // Loop over all reconstructed jets in event
1149 for(i = 0; i < TMath::Min(2, branchJet->GetEntriesFast()); ++i) //branchJet->GetEntriesFast(); ++i)
1150 {
1151
1152 jet = (Jet *)branchJet->At(i);
1153 jetMomentum = jet->P4();
1154
1155 deltaR = 999;
1156
1157 // Loop over all hard partons in event
1158 for(j = 0; j < TMath::Min(2, branchGenJet->GetEntriesFast()); ++j)
1159 {
1160 genjet = (Jet *)branchGenJet->At(j);
1161
1162 genJetMomentum = genjet->P4();
1163
1164 // this is simply to avoid warnings from initial state particle
1165 // having infite rapidity ...
1166 if(genJetMomentum.Px() == 0 && genJetMomentum.Py() == 0) continue;
1167
1168 // take the closest parton candidate
1169 if(genJetMomentum.DeltaR(jetMomentum) < deltaR)
1170 {
1171 deltaR = genJetMomentum.DeltaR(jetMomentum);
1172 bestGenJetMomentum = genJetMomentum;
1173 }
1174 }
1175
1176 if(deltaR < 0.3)
1177 {
1178
1179 pt = genJetMomentum.E();
1180 eta = genJetMomentum.Eta();
1181
1182 for(bin = 0; bin < Nbins; bin++)
1183 {
1184 if(eta > histos->at(bin).etamin && eta < histos->at(bin).etamax && pt < eMax && pt > eMin)
1185 {
1186 histos->at(bin).resolHist->Fill(jetMomentum.E() / bestGenJetMomentum.E());
1187 }
1188 }
1189 }
1190 }
1191 }
1192}
1193
1194std::pair<Double_t, Double_t> GausFit(TH1 *hist)
1195{
1196 TF1 *f1 = new TF1("f1", "gaus", hist->GetMean() - 2 * hist->GetRMS(), hist->GetMean() + 2 * hist->GetRMS());
1197 hist->Fit("f1", "RQ");
1198
1199 TF1 *f2 = new TF1("f2", "gaus", f1->GetParameter(1) - 2 * f1->GetParameter(2), f1->GetParameter(1) + 2 * f1->GetParameter(2));
1200 hist->Fit("f2", "RQ");
1201
1202 Double_t sig = f2->GetParameter(2);
1203 Double_t sigErr = f2->GetParError(2);
1204
1205 delete f1;
1206 delete f2;
1207 return make_pair(sig, sigErr);
1208}
1209
1210TGraphErrors EresGraph(std::vector<resolPlot> *histos, bool rms = false)
1211{
1212 Int_t bin;
1213 Int_t count = 0;
1214 TGraphErrors gr = TGraphErrors(Nbins / 2);
1215 double val, error;
1216 for(bin = 0; bin < Nbins; bin++)
1217 {
1218 std::pair<Double_t, Double_t> sigvalues = GausFit(histos->at(bin).resolHist);
1219 if(rms == true)
1220 {
1221 gr.SetPoint(count, (histos->at(bin).ptmin + histos->at(bin).ptmax) / 2.0, 100 * histos->at(bin).resolHist->GetRMS());
1222 //gr.SetPointError(count,0, 100*sigvalues.second); // to correct
1223 error = 100 * histos->at(bin).resolHist->GetRMSError();
1224 val = 100 * histos->at(bin).resolHist->GetRMS();
1225 if(error > 0.2 * val) error = 0.2 * val;
1226 gr.SetPointError(count, 0, error); // to correct
1227 }
1228 else
1229 {
1230
1231 gr.SetPoint(count, (histos->at(bin).ptmin + histos->at(bin).ptmax) / 2.0, 100 * sigvalues.first);
1232 error = 100 * sigvalues.second;
1233 val = 100 * sigvalues.first;
1234 if(error > 0.2 * val) error = 0.2 * val;
1235 gr.SetPointError(count, 0, error); // to correct
1236 //gr.SetPointError(count,0, 100*sigvalues.second);
1237 }
1238 count++;
1239 }
1240
1241 return gr;
1242}
1243
1244TGraphErrors MetResGraph(std::vector<resolPlot> *histos, bool rms = false)
1245{
1246 Int_t bin;
1247 Int_t count = 0;
1248 TGraphErrors gr = TGraphErrors(Nbins / 2);
1249 double val, error;
1250 for(bin = 0; bin < Nbins; bin++)
1251 {
1252 std::pair<Double_t, Double_t> sigvalues = GausFit(histos->at(bin).resolHist);
1253 if(rms == true)
1254 {
1255 gr.SetPoint(count, (histos->at(bin).ptmin + histos->at(bin).ptmax) / 2.0, histos->at(bin).resolHist->GetRMS());
1256 error = histos->at(bin).resolHist->GetRMSError();
1257 val = histos->at(bin).resolHist->GetRMS();
1258 if(error > 0.2 * val) error = 0.2 * val;
1259 gr.SetPointError(count, 0, error); // to correct
1260 }
1261 else
1262 {
1263
1264 gr.SetPoint(count, (histos->at(bin).ptmin + histos->at(bin).ptmax) / 2.0, sigvalues.first);
1265 val = sigvalues.first;
1266 error = sigvalues.second;
1267 if(error > 0.2 * val) error = 0.2 * val;
1268 gr.SetPointError(count, 0, error);
1269 //gr.SetPointError(count,0, 100*sigvalues.second);
1270 }
1271 count++;
1272 }
1273
1274 return gr;
1275}
1276
1277TGraphErrors EresGraphVsEta(std::vector<resolPlot> *histos, bool rms = false)
1278{
1279 Int_t bin;
1280 Int_t count = 0;
1281 TGraphErrors gr = TGraphErrors(Nbins / 2);
1282 double val, error;
1283 for(bin = 0; bin < Nbins; bin++)
1284 {
1285
1286 std::pair<Double_t, Double_t> sigvalues = GausFit(histos->at(bin).resolHist);
1287 if(rms == true)
1288 {
1289 gr.SetPoint(count, (histos->at(bin).etamin + histos->at(bin).etamax) / 2.0, histos->at(bin).resolHist->GetRMS());
1290 error = 100 * histos->at(bin).resolHist->GetRMSError();
1291 val = 100 * histos->at(bin).resolHist->GetRMS();
1292 if(error > 0.2 * val) error = 0.2 * val;
1293 gr.SetPointError(count, 0, error); // to correct
1294 }
1295 else
1296 {
1297 gr.SetPoint(count, (histos->at(bin).etamin + histos->at(bin).etamax) / 2.0, 100 * sigvalues.first);
1298 val = 100 * sigvalues.first;
1299 error = 100 * sigvalues.second;
1300 if(error > 0.2 * val) error = 0.2 * val;
1301 gr.SetPointError(count, 0, error);
1302 //gr.SetPointError(count,0, 100*sigvalues.second);
1303 }
1304 count++;
1305 }
1306
1307 return gr;
1308}
1309
1310void GetMetres(std::vector<resolPlot> *histos, TClonesArray *branchScalarHT, TClonesArray *branchMet, TClonesArray *branchJet, ExRootTreeReader *treeReader)
1311{
1312
1313 Long64_t allEntries = treeReader->GetEntries();
1314
1315 cout << "** Computing resolution of " << branchMet->GetName() << " vs " << branchScalarHT->GetName() << endl;
1316
1317 MissingET *met;
1318 ScalarHT *scalarHT;
1319
1320 Long64_t entry;
1321
1322 Int_t bin;
1323 Double_t ht;
1324
1325 Jet *jet;
1326 TLorentzVector p1, p2;
1327
1328 // Loop over all events
1329 for(entry = 0; entry < allEntries; ++entry)
1330 {
1331 // Load selected branches with data from specified event
1332 treeReader->ReadEntry(entry);
1333
1334 if(entry % 10000 == 0) cout << "Event number: " << entry << endl;
1335
1336 if(branchJet->GetEntriesFast() > 1)
1337 {
1338
1339 jet = (Jet *)branchJet->At(0);
1340 p1 = jet->P4();
1341 jet = (Jet *)branchJet->At(1);
1342 p2 = jet->P4();
1343
1344 met = (MissingET *)branchMet->At(0);
1345 scalarHT = (ScalarHT *)branchScalarHT->At(0);
1346 ht = scalarHT->HT;
1347
1348 if(p1.Pt() < 0.75 * ht / 2) continue;
1349 if(p2.Pt() < 0.75 * ht / 2) continue;
1350
1351 for(bin = 0; bin < Nbins; bin++)
1352 {
1353 if(ht > histos->at(bin).ptmin && ht < histos->at(bin).ptmax)
1354 {
1355 histos->at(bin).resolHist->Fill(met->P4().Px());
1356 }
1357 }
1358 }
1359 }
1360}
1361
1362//------------------------------------------------------------------------------
1363
1364void addResoGraph(TMultiGraph *mg, TGraphErrors *gr, TLegend *leg, int style, Color_t color, TString text)
1365{
1366
1367 gr->SetLineWidth(2);
1368 gr->SetLineColor(color);
1369 gr->SetMarkerStyle(style);
1370 gr->SetMarkerColor(color);
1371 gr->SetMarkerSize(1.5);
1372
1373 std::cout << "Adding " << gr->GetName() << std::endl;
1374 mg->Add(gr);
1375 leg->AddEntry(gr, text, "p");
1376}
1377
1378void DrawAxis(TMultiGraph *mg, TLegend *leg, double xmin, double xmax, double ymin, double ymax, TString tx, TString ty, bool logx = 0, bool logy = 0)
1379{
1380 mg->SetMinimum(ymin);
1381 mg->SetMaximum(ymax);
1382 mg->GetXaxis()->SetLimits(xmin, xmax);
1383
1384 mg->GetXaxis()->SetTitle(tx);
1385 mg->GetYaxis()->SetTitle(ty);
1386
1387 mg->GetYaxis()->SetTitleSize(0.07);
1388 mg->GetXaxis()->SetTitleSize(0.07);
1389 mg->GetYaxis()->SetLabelSize(0.06);
1390 mg->GetXaxis()->SetLabelSize(0.06);
1391 mg->GetYaxis()->SetLabelOffset(0.03);
1392 mg->GetYaxis()->SetTitleOffset(1.4);
1393 mg->GetXaxis()->SetTitleOffset(1.4);
1394
1395 mg->GetXaxis()->SetTitleFont(132);
1396 mg->GetYaxis()->SetTitleFont(132);
1397 mg->GetXaxis()->SetLabelFont(132);
1398 mg->GetYaxis()->SetLabelFont(132);
1399
1400 mg->GetYaxis()->SetNdivisions(505);
1401
1402 leg->SetTextFont(132);
1403 leg->SetBorderSize(0);
1404 leg->SetShadowColor(0);
1405 leg->SetFillColor(0);
1406 leg->SetFillStyle(0);
1407
1408 gStyle->SetOptTitle(0);
1409
1410 if(logx) gPad->SetLogx();
1411 if(logy) gPad->SetLogy();
1412
1413 //gPad->SetGridx();
1414 //gPad->SetGridy();
1415 gPad->SetBottomMargin(0.2);
1416 gPad->SetLeftMargin(0.2);
1417 gPad->Modified();
1418 gPad->Update();
1419}
1420
1421void DelphesValidation(
1422 const char *inputFilePion,
1423 const char *inputFileElectron,
1424 const char *inputFileMuon,
1425 const char *inputFilePhoton,
1426 const char *inputFileNeutralHadron,
1427 const char *inputFileJet,
1428 const char *inputFileBJet,
1429 const char *inputFileCJet,
1430 const char *inputFileTauJet,
1431 const char *outputFile,
1432 const char *version)
1433{
1434
1435 TChain *chainPion = new TChain("Delphes");
1436 chainPion->Add(inputFilePion);
1437 ExRootTreeReader *treeReaderPion = new ExRootTreeReader(chainPion);
1438
1439 TChain *chainElectron = new TChain("Delphes");
1440 chainElectron->Add(inputFileElectron);
1441 ExRootTreeReader *treeReaderElectron = new ExRootTreeReader(chainElectron);
1442
1443 TChain *chainMuon = new TChain("Delphes");
1444 chainMuon->Add(inputFileMuon);
1445 ExRootTreeReader *treeReaderMuon = new ExRootTreeReader(chainMuon);
1446
1447 TChain *chainPhoton = new TChain("Delphes");
1448 chainPhoton->Add(inputFilePhoton);
1449 ExRootTreeReader *treeReaderPhoton = new ExRootTreeReader(chainPhoton);
1450
1451 TChain *chainNeutralHadron = new TChain("Delphes");
1452 chainNeutralHadron->Add(inputFileNeutralHadron);
1453 ExRootTreeReader *treeReaderNeutralHadron = new ExRootTreeReader(chainNeutralHadron);
1454
1455 TChain *chainJet = new TChain("Delphes");
1456 chainJet->Add(inputFileJet);
1457 ExRootTreeReader *treeReaderJet = new ExRootTreeReader(chainJet);
1458
1459 TChain *chainBJet = new TChain("Delphes");
1460 chainBJet->Add(inputFileBJet);
1461 ExRootTreeReader *treeReaderBJet = new ExRootTreeReader(chainBJet);
1462
1463 TChain *chainCJet = new TChain("Delphes");
1464 chainCJet->Add(inputFileCJet);
1465 ExRootTreeReader *treeReaderCJet = new ExRootTreeReader(chainCJet);
1466
1467 TChain *chainTauJet = new TChain("Delphes");
1468 chainTauJet->Add(inputFileTauJet);
1469 ExRootTreeReader *treeReaderTauJet = new ExRootTreeReader(chainTauJet);
1470
1471 TClonesArray *branchParticleElectron = treeReaderElectron->UseBranch("Particle");
1472 TClonesArray *branchTrackElectron = treeReaderElectron->UseBranch("Track");
1473 TClonesArray *branchElectron = treeReaderElectron->UseBranch("Electron");
1474 TClonesArray *branchElectronPF = treeReaderElectron->UseBranch("ElectronPF");
1475
1476 TClonesArray *branchParticleMuon = treeReaderMuon->UseBranch("Particle");
1477 TClonesArray *branchTrackMuon = treeReaderMuon->UseBranch("Track");
1478 TClonesArray *branchMuon = treeReaderMuon->UseBranch("Muon");
1479
1480 TClonesArray *branchParticlePion = treeReaderPion->UseBranch("Particle");
1481 TClonesArray *branchTrackPion = treeReaderPion->UseBranch("Track");
1482 TClonesArray *branchPion = treeReaderPion->UseBranch("Pion");
1483
1484 TClonesArray *branchParticlePhoton = treeReaderPhoton->UseBranch("Particle");
1485 TClonesArray *branchTowerPhoton = treeReaderPhoton->UseBranch("Tower");
1486 TClonesArray *branchPhoton = treeReaderPhoton->UseBranch("Photon");
1487
1488 TClonesArray *branchParticleNeutralHadron = treeReaderNeutralHadron->UseBranch("Particle");
1489 TClonesArray *branchTowerNeutralHadron = treeReaderNeutralHadron->UseBranch("Tower");
1490
1491 TClonesArray *branchGenJet = treeReaderJet->UseBranch("GenJet");
1492 TClonesArray *branchParticleJet = treeReaderJet->UseBranch("Particle");
1493 TClonesArray *branchPFJet = treeReaderJet->UseBranch("PFJet");
1494 TClonesArray *branchCaloJet = treeReaderJet->UseBranch("CaloJet");
1495 TClonesArray *branchJet = treeReaderJet->UseBranch("Jet");
1496
1497 TClonesArray *branchParticleBJet = treeReaderBJet->UseBranch("Particle");
1498 TClonesArray *branchPFBJet = treeReaderBJet->UseBranch("Jet");
1499
1500 TClonesArray *branchParticleCJet = treeReaderCJet->UseBranch("Particle");
1501 TClonesArray *branchPFCJet = treeReaderCJet->UseBranch("Jet");
1502
1503 TClonesArray *branchParticleTauJet = treeReaderTauJet->UseBranch("Particle");
1504 TClonesArray *branchPFTauJet = treeReaderTauJet->UseBranch("Jet");
1505
1506 TClonesArray *branchGenScalarHT = treeReaderJet->UseBranch("GenScalarHT");
1507 TClonesArray *branchMet = treeReaderJet->UseBranch("PFMissingET");
1508 TClonesArray *branchCaloMet = treeReaderJet->UseBranch("CaloMissingET");
1509
1510 std::vector<Color_t> colors;
1511
1512 colors.push_back(kBlack);
1513 colors.push_back(kBlue);
1514 colors.push_back(kRed);
1515 colors.push_back(kGreen + 1);
1516 colors.push_back(kMagenta + 1);
1517 colors.push_back(kOrange);
1518
1519 std::vector<Int_t> markerStyles;
1520
1521 markerStyles.push_back(20);
1522 markerStyles.push_back(21);
1523 markerStyles.push_back(22);
1524 markerStyles.push_back(23);
1525 markerStyles.push_back(33);
1526 markerStyles.push_back(34);
1527
1528 TString pdfOutput(outputFile);
1529 pdfOutput.ReplaceAll(".root", ".pdf");
1530
1531 TString figPath = inputFilePion;
1532 figPath.ReplaceAll(".root", "");
1533 figPath.ReplaceAll("root", "www/fig");
1534 Int_t lastSlash = figPath.Last('/');
1535 Int_t sizePath = figPath.Length();
1536 figPath.Remove(lastSlash + 1, sizePath);
1537
1538 TString header = pdfOutput;
1539 header.ReplaceAll(".pdf", "");
1540 header.ReplaceAll("validation_", "");
1541 lastSlash = header.Last('/');
1542 sizePath = header.Length();
1543 header.Remove(0, lastSlash + 1);
1544
1545 TString vrs(version);
1546
1547 TPaveText *pave = new TPaveText(0.0, 0.89, 0.94, 0.94, "NDC");
1548 pave->SetTextAlign(30);
1549 pave->SetTextFont(132);
1550 pave->SetBorderSize(0);
1551 pave->SetShadowColor(0);
1552 pave->SetFillColor(0);
1553 pave->SetFillStyle(0);
1554 pave->AddText("Delphes " + vrs + " - " + header);
1555
1556 TString s_etaMin, s_etaMax, s_eta, s_pt, s_e;
1557
1558 Double_t ptMin = 1.;
1559 Double_t ptMax = 50000.;
1560
1561 Double_t etaMin = -6.;
1562 Double_t etaMax = 6.;
1563
1564 std::vector<Double_t> ptVals;
1565
1566 ptVals.push_back(1.);
1567 ptVals.push_back(10.);
1568 ptVals.push_back(100.);
1569 ptVals.push_back(1000.);
1570 ptVals.push_back(10000.);
1571
1572 std::vector<Double_t> etaVals;
1573
1574 etaVals.push_back(0.);
1575 etaVals.push_back(1.5);
1576 etaVals.push_back(2.5);
1577 etaVals.push_back(4.0);
1578 etaVals.push_back(6.0);
1579
1580 const int n_etabins = etaVals.size() - 1;
1581 const int n_ptbins = ptVals.size();
1582
1583 //////////////////////////
1584 // Tracking performance //
1585 //////////////////////////
1586
1587 // --------- Pion Tracks --------- //
1588
1589 TMultiGraph *mg_trkpi_res_pt = new TMultiGraph("", "");
1590 TMultiGraph *mg_trkpi_eff_pt = new TMultiGraph("", "");
1591 TMultiGraph *mg_trkpi_res_eta = new TMultiGraph("", "");
1592 TMultiGraph *mg_trkpi_eff_eta = new TMultiGraph("", "");
1593
1594 TLegend *leg_trkpi_res_pt = new TLegend(0.55, 0.22, 0.90, 0.48);
1595 TLegend *leg_trkpi_eff_pt = (TLegend *)leg_trkpi_res_pt->Clone();
1596 TLegend *leg_trkpi_res_eta = (TLegend *)leg_trkpi_res_pt->Clone();
1597 TLegend *leg_trkpi_eff_eta = (TLegend *)leg_trkpi_res_eta->Clone();
1598
1599 TGraphErrors *gr_trkpi_res_pt = new TGraphErrors[n_etabins];
1600 TGraphErrors *gr_trkpi_eff_pt = new TGraphErrors[n_etabins];
1601 TGraphErrors *gr_trkpi_res_eta = new TGraphErrors[n_ptbins];
1602 TGraphErrors *gr_trkpi_eff_eta = new TGraphErrors[n_ptbins];
1603 TH1D *h_trkpi_eff_pt, *h_trkpi_eff_eta;
1604
1605 std::vector<resolPlot> *plots_trkpi_res_pt = new std::vector<resolPlot>[n_etabins];
1606 std::vector<resolPlot> *plots_trkpi_res_eta = new std::vector<resolPlot>[n_ptbins];
1607
1608 // loop over eta bins
1609 for(k = 0; k < etaVals.size() - 1; k++)
1610 {
1611 HistogramsCollection(&plots_trkpi_res_pt[k], TMath::Log10(ptMin), TMath::Log10(ptMax), "trkpi");
1612 GetPtres<Track>(&plots_trkpi_res_pt[k], branchTrackPion, branchParticlePion, 211, etaVals.at(k), etaVals.at(k + 1), treeReaderPion);
1613 gr_trkpi_res_pt[k] = EresGraph(&plots_trkpi_res_pt[k]);
1614
1615 h_trkpi_eff_pt = GetEffPt<Track>(branchTrackPion, branchParticlePion, "Pion", 211, ptMin, ptMax, etaVals.at(k), etaVals.at(k + 1), treeReaderPion);
1616 gr_trkpi_eff_pt[k] = TGraphErrors(h_trkpi_eff_pt);
1617
1618 s_etaMin = Form("%.1f", etaVals.at(k));
1619 s_etaMax = Form("%.1f", etaVals.at(k + 1));
1620
1621 s_eta = "#pi^{ #pm} , " + s_etaMin + " < | #eta | < " + s_etaMax;
1622
1623 gr_trkpi_res_pt[k].SetName("trkRes_" + s_etaMin + "_" + s_etaMax);
1624 gr_trkpi_eff_pt[k].SetName("trkEff_" + s_etaMin + "_" + s_etaMax);
1625
1626 addResoGraph(mg_trkpi_res_pt, &gr_trkpi_res_pt[k], leg_trkpi_res_pt, markerStyles.at(k), colors.at(k), s_eta);
1627 addResoGraph(mg_trkpi_eff_pt, &gr_trkpi_eff_pt[k], leg_trkpi_eff_pt, markerStyles.at(k), colors.at(k), s_eta);
1628 }
1629
1630 // loop over pt
1631 for(k = 0; k < ptVals.size(); k++)
1632 {
1633 HistogramsCollectionVsEta(&plots_trkpi_res_eta[k], etaMin, etaMax, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k), "trkpi", 0.0, 2.0);
1634 GetPtresVsEta<Track>(&plots_trkpi_res_eta[k], branchTrackPion, branchParticlePion, 211, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k), treeReaderPion);
1635 gr_trkpi_res_eta[k] = EresGraphVsEta(&plots_trkpi_res_eta[k]);
1636
1637 h_trkpi_eff_eta = GetEffEta<Track>(branchTrackPion, branchParticlePion, "Pion", 211, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k), etaMin, etaMax, treeReaderPion);
1638 gr_trkpi_eff_eta[k] = TGraphErrors(h_trkpi_eff_eta);
1639
1640 s_pt = Form("#pi^{ #pm} , p_{T} = %.0f GeV", ptVals.at(k));
1641 if(ptVals.at(k) >= 1000.) s_pt = Form("#pi^{ #pm} , p_{T} = %.0f TeV", ptVals.at(k) / 1000.);
1642
1643 addResoGraph(mg_trkpi_res_eta, &gr_trkpi_res_eta[k], leg_trkpi_res_eta, markerStyles.at(k), colors.at(k), s_pt);
1644 addResoGraph(mg_trkpi_eff_eta, &gr_trkpi_eff_eta[k], leg_trkpi_eff_eta, markerStyles.at(k), colors.at(k), s_pt);
1645 }
1646
1647 TCanvas *c_trkpi_res_pt = new TCanvas("", "", 800, 600);
1648
1649 mg_trkpi_res_pt->Draw("APE");
1650 DrawAxis(mg_trkpi_res_pt, leg_trkpi_res_pt, ptMin, ptMax, 0.01, 100, "p_{T} [GeV]", "(track resolution in p_{T})/p_{T} (%)", true, true);
1651 leg_trkpi_res_pt->Draw();
1652 pave->Draw();
1653
1654 c_trkpi_res_pt->Print(pdfOutput + "(", "pdf");
1655 c_trkpi_res_pt->Print(figPath + "img_trkpi_res_pt.pdf", "pdf");
1656 c_trkpi_res_pt->Print(figPath + "img_trkpi_res_pt.png", "png");
1657
1658 TCanvas *c_trkpi_res_eta = new TCanvas("", "", 800, 600);
1659
1660 mg_trkpi_res_eta->Draw("APE");
1661 DrawAxis(mg_trkpi_res_eta, leg_trkpi_res_eta, etaMin, etaMax, 0.01, 100, " #eta ", "(track resolution in p_{T})/p_{T} (%)", false, true);
1662 leg_trkpi_res_eta->Draw();
1663 pave->Draw();
1664
1665 c_trkpi_res_eta->Print(pdfOutput, "pdf");
1666 c_trkpi_res_eta->Print(figPath + "img_trkpi_res_eta.pdf", "pdf");
1667 c_trkpi_res_eta->Print(figPath + "img_trkpi_res_eta.png", "png");
1668
1669 TCanvas *c_trkpi_eff_pt = new TCanvas("", "", 800, 600);
1670
1671 mg_trkpi_eff_pt->Draw("APE");
1672 DrawAxis(mg_trkpi_eff_pt, leg_trkpi_eff_pt, ptMin, ptMax, 0.0, 100, "p_{T} [GeV]", "tracking efficiency (%)", true, false);
1673 leg_trkpi_eff_pt->Draw();
1674 pave->Draw();
1675
1676 c_trkpi_eff_pt->Print(pdfOutput, "pdf");
1677 c_trkpi_eff_pt->Print(figPath + "img_trkpi_eff_pt.pdf", "pdf");
1678 c_trkpi_eff_pt->Print(figPath + "img_trkpi_eff_pt.png", "png");
1679
1680 TCanvas *c_trkpi_eff_eta = new TCanvas("", "", 800, 600);
1681
1682 mg_trkpi_eff_eta->Draw("APE");
1683 DrawAxis(mg_trkpi_eff_eta, leg_trkpi_eff_eta, etaMin, etaMax, 0.0, 100, " #eta ", "tracking efficiency (%)", false, false);
1684 leg_trkpi_eff_eta->Draw();
1685 pave->Draw();
1686
1687 c_trkpi_eff_eta->Print(pdfOutput, "pdf");
1688 c_trkpi_eff_eta->Print(figPath + "img_trkpi_eff_eta.pdf", "pdf");
1689 c_trkpi_eff_eta->Print(figPath + "img_trkpi_eff_eta.png", "png");
1690
1691 // --------- Electron Tracks --------- //
1692
1693 TMultiGraph *mg_trkele_res_pt = new TMultiGraph("", "");
1694 TMultiGraph *mg_trkele_eff_pt = new TMultiGraph("", "");
1695 TMultiGraph *mg_trkele_res_eta = new TMultiGraph("", "");
1696 TMultiGraph *mg_trkele_eff_eta = new TMultiGraph("", "");
1697
1698 TLegend *leg_trkele_res_pt = new TLegend(0.55, 0.22, 0.90, 0.48);
1699 TLegend *leg_trkele_eff_pt = (TLegend *)leg_trkele_res_pt->Clone();
1700 TLegend *leg_trkele_res_eta = (TLegend *)leg_trkele_res_pt->Clone();
1701 TLegend *leg_trkele_eff_eta = (TLegend *)leg_trkele_res_eta->Clone();
1702
1703 TGraphErrors *gr_trkele_res_pt = new TGraphErrors[n_etabins];
1704 TGraphErrors *gr_trkele_eff_pt = new TGraphErrors[n_etabins];
1705 TGraphErrors *gr_trkele_res_eta = new TGraphErrors[n_ptbins];
1706 TGraphErrors *gr_trkele_eff_eta = new TGraphErrors[n_ptbins];
1707
1708 TH1D *h_trkele_eff_pt, *h_trkele_eff_eta;
1709
1710 std::vector<resolPlot> *plots_trkele_res_pt = new std::vector<resolPlot>[n_etabins];
1711 std::vector<resolPlot> *plots_trkele_res_eta = new std::vector<resolPlot>[n_ptbins];
1712
1713 // loop over eta bins
1714 for(k = 0; k < etaVals.size() - 1; k++)
1715 {
1716 HistogramsCollection(&plots_trkele_res_pt[k], TMath::Log10(ptMin), TMath::Log10(ptMax), "trkele");
1717 GetPtres<Track>(&plots_trkele_res_pt[k], branchTrackElectron, branchParticleElectron, 11, etaVals.at(k), etaVals.at(k + 1), treeReaderElectron);
1718 gr_trkele_res_pt[k] = EresGraph(&plots_trkele_res_pt[k]);
1719
1720 h_trkele_eff_pt = GetEffPt<Track>(branchTrackElectron, branchParticleElectron, "Electron", 11, ptMin, ptMax, etaVals.at(k), etaVals.at(k + 1), treeReaderElectron);
1721 gr_trkele_eff_pt[k] = TGraphErrors(h_trkele_eff_pt);
1722
1723 s_etaMin = Form("%.1f", etaVals.at(k));
1724 s_etaMax = Form("%.1f", etaVals.at(k + 1));
1725
1726 s_eta = "e^{ #pm} , " + s_etaMin + " < | #eta | < " + s_etaMax;
1727
1728 gr_trkele_res_pt[k].SetName("trkRes_" + s_etaMin + "_" + s_etaMax);
1729 gr_trkele_eff_pt[k].SetName("trkEff_" + s_etaMin + "_" + s_etaMax);
1730
1731 addResoGraph(mg_trkele_res_pt, &gr_trkele_res_pt[k], leg_trkele_res_pt, markerStyles.at(k), colors.at(k), s_eta);
1732 addResoGraph(mg_trkele_eff_pt, &gr_trkele_eff_pt[k], leg_trkele_eff_pt, markerStyles.at(k), colors.at(k), s_eta);
1733 }
1734
1735 // loop over pt
1736 for(k = 0; k < ptVals.size(); k++)
1737 {
1738 HistogramsCollectionVsEta(&plots_trkele_res_eta[k], etaMin, etaMax, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k), "trkele", 0.0, 2.0);
1739 GetPtresVsEta<Track>(&plots_trkele_res_eta[k], branchTrackElectron, branchParticleElectron, 11, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k), treeReaderElectron);
1740 gr_trkele_res_eta[k] = EresGraphVsEta(&plots_trkele_res_eta[k]);
1741
1742 h_trkele_eff_eta = GetEffEta<Track>(branchTrackElectron, branchParticleElectron, "Electron", 11, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k), etaMin, etaMax, treeReaderElectron);
1743 gr_trkele_eff_eta[k] = TGraphErrors(h_trkele_eff_eta);
1744
1745 s_pt = Form("e^{ #pm} , p_{T} = %.0f GeV", ptVals.at(k));
1746 if(ptVals.at(k) >= 1000.) s_pt = Form("e^{ #pm} , p_{T} = %.0f TeV", ptVals.at(k) / 1000.);
1747
1748 addResoGraph(mg_trkele_res_eta, &gr_trkele_res_eta[k], leg_trkele_res_eta, markerStyles.at(k), colors.at(k), s_pt);
1749 addResoGraph(mg_trkele_eff_eta, &gr_trkele_eff_eta[k], leg_trkele_eff_eta, markerStyles.at(k), colors.at(k), s_pt);
1750 }
1751
1752 TCanvas *c_trkele_res_pt = new TCanvas("", "", 800, 600);
1753
1754 mg_trkele_res_pt->Draw("APE");
1755 DrawAxis(mg_trkele_res_pt, leg_trkele_res_pt, ptMin, ptMax, 0.01, 100, "p_{T} [GeV]", "(track resolution in p_{T})/p_{T} (%)", true, true);
1756 leg_trkele_res_pt->Draw();
1757 pave->Draw();
1758
1759 c_trkele_res_pt->Print(pdfOutput, "pdf");
1760 c_trkele_res_pt->Print(figPath + "img_trkele_res_pt.pdf", "pdf");
1761 c_trkele_res_pt->Print(figPath + "img_trkele_res_pt.png", "png");
1762
1763 TCanvas *c_trkele_res_eta = new TCanvas("", "", 800, 600);
1764
1765 mg_trkele_res_eta->Draw("APE");
1766 DrawAxis(mg_trkele_res_eta, leg_trkele_res_eta, etaMin, etaMax, 0.01, 100, " #eta ", "(track resolution in p_{T})/p_{T} (%)", false, true);
1767 leg_trkele_res_eta->Draw();
1768 pave->Draw();
1769
1770 c_trkele_res_eta->Print(pdfOutput, "pdf");
1771 c_trkele_res_eta->Print(figPath + "img_trkele_res_eta.pdf", "pdf");
1772 c_trkele_res_eta->Print(figPath + "img_trkele_res_eta.png", "png");
1773
1774 TCanvas *c_trkele_eff_pt = new TCanvas("", "", 800, 600);
1775
1776 mg_trkele_eff_pt->Draw("APE");
1777 DrawAxis(mg_trkele_eff_pt, leg_trkele_eff_pt, ptMin, ptMax, 0.0, 100, "p_{T} [GeV]", "tracking efficiency (%)", true, false);
1778 leg_trkele_eff_pt->Draw();
1779 pave->Draw();
1780
1781 c_trkele_eff_pt->Print(pdfOutput, "pdf");
1782 c_trkele_eff_pt->Print(figPath + "img_trkele_eff_pt.pdf", "pdf");
1783 c_trkele_eff_pt->Print(figPath + "img_trkele_eff_pt.png", "png");
1784
1785 TCanvas *c_trkele_eff_eta = new TCanvas("", "", 800, 600);
1786
1787 mg_trkele_eff_eta->Draw("APE");
1788 DrawAxis(mg_trkele_eff_eta, leg_trkele_eff_eta, etaMin, etaMax, 0.0, 100, " #eta ", "tracking efficiency (%)", false, false);
1789 leg_trkele_eff_eta->Draw();
1790 pave->Draw();
1791
1792 c_trkele_eff_eta->Print(pdfOutput, "pdf");
1793 c_trkele_eff_eta->Print(figPath + "img_trkele_eff_eta.pdf", "pdf");
1794 c_trkele_eff_eta->Print(figPath + "img_trkele_eff_eta.png", "png");
1795
1796 // --------- Muon Tracks --------- //
1797
1798 TMultiGraph *mg_trkmu_res_pt = new TMultiGraph("", "");
1799 TMultiGraph *mg_trkmu_eff_pt = new TMultiGraph("", "");
1800 TMultiGraph *mg_trkmu_res_eta = new TMultiGraph("", "");
1801 TMultiGraph *mg_trkmu_eff_eta = new TMultiGraph("", "");
1802
1803 TLegend *leg_trkmu_res_pt = new TLegend(0.55, 0.22, 0.90, 0.48);
1804 TLegend *leg_trkmu_eff_pt = (TLegend *)leg_trkmu_res_pt->Clone();
1805
1806 TLegend *leg_trkmu_res_eta = (TLegend *)leg_trkmu_res_pt->Clone();
1807 TLegend *leg_trkmu_eff_eta = (TLegend *)leg_trkmu_res_eta->Clone();
1808
1809 TGraphErrors *gr_trkmu_res_pt = new TGraphErrors[n_etabins];
1810 TGraphErrors *gr_trkmu_eff_pt = new TGraphErrors[n_etabins];
1811 TGraphErrors *gr_trkmu_res_eta = new TGraphErrors[n_ptbins];
1812 TGraphErrors *gr_trkmu_eff_eta = new TGraphErrors[n_ptbins];
1813
1814 TH1D *h_trkmu_eff_pt, *h_trkmu_eff_eta;
1815
1816 std::vector<resolPlot> *plots_trkmu_res_pt = new std::vector<resolPlot>[n_etabins];
1817 std::vector<resolPlot> *plots_trkmu_res_eta = new std::vector<resolPlot>[n_ptbins];
1818
1819 // loop over eta bins
1820 for(k = 0; k < etaVals.size() - 1; k++)
1821 {
1822 HistogramsCollection(&plots_trkmu_res_pt[k], TMath::Log10(ptMin), TMath::Log10(ptMax), "trkmu");
1823 GetPtres<Track>(&plots_trkmu_res_pt[k], branchTrackMuon, branchParticleMuon, 13, etaVals.at(k), etaVals.at(k + 1), treeReaderMuon);
1824 gr_trkmu_res_pt[k] = EresGraph(&plots_trkmu_res_pt[k]);
1825
1826 h_trkmu_eff_pt = GetEffPt<Track>(branchTrackMuon, branchParticleMuon, "Muon", 13, ptMin, ptMax, etaVals.at(k), etaVals.at(k + 1), treeReaderMuon);
1827 gr_trkmu_eff_pt[k] = TGraphErrors(h_trkmu_eff_pt);
1828
1829 s_etaMin = Form("%.1f", etaVals.at(k));
1830 s_etaMax = Form("%.1f", etaVals.at(k + 1));
1831
1832 s_eta = "#mu^{ #pm} , " + s_etaMin + " < | #eta | < " + s_etaMax;
1833
1834 gr_trkmu_res_pt[k].SetName("trkRes_" + s_etaMin + "_" + s_etaMax);
1835 gr_trkmu_eff_pt[k].SetName("trkEff_" + s_etaMin + "_" + s_etaMax);
1836
1837 addResoGraph(mg_trkmu_res_pt, &gr_trkmu_res_pt[k], leg_trkmu_res_pt, markerStyles.at(k), colors.at(k), s_eta);
1838 addResoGraph(mg_trkmu_eff_pt, &gr_trkmu_eff_pt[k], leg_trkmu_eff_pt, markerStyles.at(k), colors.at(k), s_eta);
1839 }
1840
1841 // loop over pt
1842 for(k = 0; k < ptVals.size(); k++)
1843 {
1844 HistogramsCollectionVsEta(&plots_trkmu_res_eta[k], etaMin, etaMax, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k), "trkmu", 0.0, 2.0);
1845 GetPtresVsEta<Track>(&plots_trkmu_res_eta[k], branchTrackMuon, branchParticleMuon, 13, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k), treeReaderMuon);
1846 gr_trkmu_res_eta[k] = EresGraphVsEta(&plots_trkmu_res_eta[k]);
1847
1848 h_trkmu_eff_eta = GetEffEta<Track>(branchTrackMuon, branchParticleMuon, "Muon", 13, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k), etaMin, etaMax, treeReaderMuon);
1849 gr_trkmu_eff_eta[k] = TGraphErrors(h_trkmu_eff_eta);
1850
1851 s_pt = Form("#mu^{ #pm} , p_{T} = %.0f GeV", ptVals.at(k));
1852 if(ptVals.at(k) >= 1000.) s_pt = Form("#mu^{ #pm} , p_{T} = %.0f TeV", ptVals.at(k) / 1000.);
1853
1854 addResoGraph(mg_trkmu_res_eta, &gr_trkmu_res_eta[k], leg_trkmu_res_eta, markerStyles.at(k), colors.at(k), s_pt);
1855 addResoGraph(mg_trkmu_eff_eta, &gr_trkmu_eff_eta[k], leg_trkmu_eff_eta, markerStyles.at(k), colors.at(k), s_pt);
1856 }
1857
1858 TCanvas *c_trkmu_res_pt = new TCanvas("", "", 800, 600);
1859
1860 mg_trkmu_res_pt->Draw("APE");
1861 DrawAxis(mg_trkmu_res_pt, leg_trkmu_res_pt, ptMin, ptMax, 0.01, 100, "p_{T} [GeV]", "(track resolution in p_{T})/p_{T} (%)", true, true);
1862 leg_trkmu_res_pt->Draw();
1863 pave->Draw();
1864
1865 c_trkmu_res_pt->Print(pdfOutput, "pdf");
1866 c_trkmu_res_pt->Print(figPath + "img_trkmu_res_pt.pdf", "pdf");
1867 c_trkmu_res_pt->Print(figPath + "img_trkmu_res_pt.png", "png");
1868
1869 TCanvas *c_trkmu_res_eta = new TCanvas("", "", 800, 600);
1870
1871 mg_trkmu_res_eta->Draw("APE");
1872 DrawAxis(mg_trkmu_res_eta, leg_trkmu_res_eta, etaMin, etaMax, 0.01, 100, " #eta ", "(track resolution in p_{T})/p_{T} (%)", false, true);
1873 leg_trkmu_res_eta->Draw();
1874 pave->Draw();
1875
1876 c_trkmu_res_eta->Print(pdfOutput, "pdf");
1877 c_trkmu_res_eta->Print(figPath + "img_trkmu_res_eta.pdf", "pdf");
1878 c_trkmu_res_eta->Print(figPath + "img_trkmu_res_eta.png", "png");
1879
1880 TCanvas *c_trkmu_eff_pt = new TCanvas("", "", 800, 600);
1881
1882 mg_trkmu_eff_pt->Draw("APE");
1883 DrawAxis(mg_trkmu_eff_pt, leg_trkmu_eff_pt, ptMin, ptMax, 0.0, 100, "p_{T} [GeV]", "tracking efficiency (%)", true, false);
1884 leg_trkmu_eff_pt->Draw();
1885 pave->Draw();
1886
1887 c_trkmu_eff_pt->Print(pdfOutput, "pdf");
1888 c_trkmu_eff_pt->Print(figPath + "img_trkmu_eff_pt.pdf", "pdf");
1889 c_trkmu_eff_pt->Print(figPath + "img_trkmu_eff_pt.png", "png");
1890
1891 TCanvas *c_trkmu_eff_eta = new TCanvas("", "", 800, 600);
1892
1893 mg_trkmu_eff_eta->Draw("APE");
1894 DrawAxis(mg_trkmu_eff_eta, leg_trkmu_eff_eta, etaMin, etaMax, 0.0, 100, " #eta ", "tracking efficiency (%)", false, false);
1895 leg_trkmu_eff_eta->Draw();
1896 pave->Draw();
1897
1898 c_trkmu_eff_eta->Print(pdfOutput, "pdf");
1899 c_trkmu_eff_eta->Print(figPath + "img_trkmu_eff_eta.pdf", "pdf");
1900 c_trkmu_eff_eta->Print(figPath + "img_trkmu_eff_eta.png", "png");
1901
1902 //////////////////////
1903 // Ecal performance //
1904 //////////////////////
1905
1906 TMultiGraph *mg_ecal_res_e = new TMultiGraph("", "");
1907 TMultiGraph *mg_ecal_res_eta = new TMultiGraph("", "");
1908
1909 TLegend *leg_ecal_res_e = new TLegend(0.55, 0.64, 0.90, 0.90);
1910 TLegend *leg_ecal_res_eta = new TLegend(0.60, 0.59, 0.95, 0.90);
1911
1912 TGraphErrors *gr_ecal_res_e = new TGraphErrors[n_etabins];
1913 TGraphErrors *gr_ecal_res_eta = new TGraphErrors[n_ptbins];
1914
1915 std::vector<resolPlot> *plots_ecal_res_e = new std::vector<resolPlot>[n_etabins];
1916 std::vector<resolPlot> *plots_ecal_res_eta = new std::vector<resolPlot>[n_ptbins];
1917
1918 // loop over eta bins
1919 for(k = 0; k < etaVals.size() - 1; k++)
1920 {
1921 HistogramsCollection(&plots_ecal_res_e[k], TMath::Log10(ptMin), TMath::Log10(ptMax), "ecal");
1922 GetEres<Tower>(&plots_ecal_res_e[k], branchTowerPhoton, branchParticlePhoton, 22, etaVals.at(k), etaVals.at(k + 1), treeReaderPhoton);
1923 gr_ecal_res_e[k] = EresGraph(&plots_ecal_res_e[k]);
1924
1925 s_etaMin = Form("%.1f", etaVals.at(k));
1926 s_etaMax = Form("%.1f", etaVals.at(k + 1));
1927
1928 s_eta = "#gamma , " + s_etaMin + " < | #eta | < " + s_etaMax;
1929
1930 gr_ecal_res_e[k].SetName("trkRes_" + s_etaMin + "_" + s_etaMax);
1931
1932 addResoGraph(mg_ecal_res_e, &gr_ecal_res_e[k], leg_ecal_res_e, markerStyles.at(k), colors.at(k), s_eta);
1933 }
1934
1935 // loop over pt
1936 for(k = 0; k < ptVals.size(); k++)
1937 {
1938 HistogramsCollectionVsEta(&plots_ecal_res_eta[k], etaMin, etaMax, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k), "ecal", 0.0, 2.0);
1939 GetEresVsEta<Tower>(&plots_ecal_res_eta[k], branchTowerPhoton, branchParticlePhoton, 22, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k), treeReaderPhoton);
1940 gr_ecal_res_eta[k] = EresGraphVsEta(&plots_ecal_res_eta[k]);
1941
1942 s_e = Form("#gamma , E = %.0f GeV", ptVals.at(k));
1943 if(ptVals.at(k) >= 1000.) s_e = Form("#gamma , E = %.0f TeV", ptVals.at(k) / 1000.);
1944
1945 addResoGraph(mg_ecal_res_eta, &gr_ecal_res_eta[k], leg_ecal_res_eta, markerStyles.at(k), colors.at(k), s_e);
1946 }
1947
1948 TCanvas *c_ecal_res_e = new TCanvas("", "", 800, 600);
1949
1950 mg_ecal_res_e->Draw("APE");
1951 // DrawAxis(mg_ecal_res_e, leg_ecal_res_e, ptMin, ptMax, 0.5, 100, "E [GeV]", "(ECAL resolution in E)/E (%)", true, true);
1952 DrawAxis(mg_ecal_res_e, leg_ecal_res_e, ptMin, ptMax, 0.0, 20, "E [GeV]", "(ECAL resolution in E)/E (%)", true, false);
1953 leg_ecal_res_e->Draw();
1954 pave->Draw();
1955
1956 c_ecal_res_e->Print(pdfOutput, "pdf");
1957 c_ecal_res_e->Print(figPath + "img_ecal_res_e.pdf", "pdf");
1958 c_ecal_res_e->Print(figPath + "img_ecal_res_e.png", "png");
1959
1960 TCanvas *c_ecal_res_eta = new TCanvas("", "", 800, 600);
1961
1962 mg_ecal_res_eta->Draw("APE");
1963 //DrawAxis(mg_ecal_res_eta, leg_ecal_res_eta, etaMin, etaMax, 0.5, 100, " #eta ", "(ECAL resolution in E)/E (%)", false, true);
1964 DrawAxis(mg_ecal_res_eta, leg_ecal_res_eta, etaMin, etaMax, 0.0, 20, " #eta ", "(ECAL resolution in E)/E (%)", false, false);
1965 leg_ecal_res_eta->Draw();
1966 pave->Draw();
1967
1968 c_ecal_res_eta->Print(pdfOutput, "pdf");
1969 c_ecal_res_eta->Print(figPath + "img_ecal_res_eta.pdf", "pdf");
1970 c_ecal_res_eta->Print(figPath + "img_ecal_res_eta.png", "png");
1971
1972 //////////////////////
1973 // Hcal performance //
1974 //////////////////////
1975
1976 TMultiGraph *mg_hcal_res_e = new TMultiGraph("", "");
1977 TMultiGraph *mg_hcal_res_eta = new TMultiGraph("", "");
1978
1979 TLegend *leg_hcal_res_e = new TLegend(0.55, 0.64, 0.90, 0.90);
1980 TLegend *leg_hcal_res_eta = new TLegend(0.60, 0.59, 0.95, 0.90);
1981
1982 TGraphErrors *gr_hcal_res_e = new TGraphErrors[n_etabins];
1983 TGraphErrors *gr_hcal_res_eta = new TGraphErrors[n_ptbins];
1984
1985 std::vector<resolPlot> *plots_hcal_res_e = new std::vector<resolPlot>[n_etabins];
1986 std::vector<resolPlot> *plots_hcal_res_eta = new std::vector<resolPlot>[n_ptbins];
1987
1988 // loop over eta bins
1989 for(k = 0; k < etaVals.size() - 1; k++)
1990 {
1991 HistogramsCollection(&plots_hcal_res_e[k], TMath::Log10(ptMin), TMath::Log10(ptMax), "hcal");
1992 GetEres<Tower>(&plots_hcal_res_e[k], branchTowerNeutralHadron, branchParticleNeutralHadron, 2112, etaVals.at(k), etaVals.at(k + 1), treeReaderNeutralHadron);
1993
1994 gr_hcal_res_e[k] = EresGraph(&plots_hcal_res_e[k]);
1995
1996 s_etaMin = Form("%.1f", etaVals.at(k));
1997 s_etaMax = Form("%.1f", etaVals.at(k + 1));
1998
1999 s_eta = "n , " + s_etaMin + " < | #eta | < " + s_etaMax;
2000
2001 gr_hcal_res_e[k].SetName("trkRes_" + s_etaMin + "_" + s_etaMax);
2002
2003 addResoGraph(mg_hcal_res_e, &gr_hcal_res_e[k], leg_hcal_res_e, markerStyles.at(k), colors.at(k), s_eta);
2004 }
2005
2006 // loop over pt
2007 for(k = 0; k < ptVals.size(); k++)
2008 {
2009 HistogramsCollectionVsEta(&plots_hcal_res_eta[k], etaMin, etaMax, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k), "hcal", 0.0, 2.0);
2010 GetEresVsEta<Tower>(&plots_hcal_res_eta[k], branchTowerNeutralHadron, branchParticleNeutralHadron, 2112, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k), treeReaderNeutralHadron);
2011 gr_hcal_res_eta[k] = EresGraphVsEta(&plots_hcal_res_eta[k]);
2012
2013 s_e = Form("n , E = %.0f GeV", ptVals.at(k));
2014 if(ptVals.at(k) >= 1000.) s_e = Form("n , E = %.0f TeV", ptVals.at(k) / 1000.);
2015
2016 addResoGraph(mg_hcal_res_eta, &gr_hcal_res_eta[k], leg_hcal_res_eta, markerStyles.at(k), colors.at(k), s_e);
2017 }
2018
2019 TCanvas *c_hcal_res_e = new TCanvas("", "", 800, 600);
2020
2021 mg_hcal_res_e->Draw("APE");
2022 //DrawAxis(mg_hcal_res_e, leg_hcal_res_e, ptMin, ptMax, 1, 100, "E [GeV]", "(HCAL resolution in E)/E (%)", true, true);
2023 DrawAxis(mg_hcal_res_e, leg_hcal_res_e, ptMin, ptMax, 0.0, 50, "E [GeV]", "(HCAL resolution in E)/E (%)", true, false);
2024 leg_hcal_res_e->Draw();
2025 pave->Draw();
2026
2027 c_hcal_res_e->Print(pdfOutput, "pdf");
2028 c_hcal_res_e->Print(figPath + "img_hcal_res_e.pdf", "pdf");
2029 c_hcal_res_e->Print(figPath + "img_hcal_res_e.png", "png");
2030
2031 TCanvas *c_hcal_res_eta = new TCanvas("", "", 800, 600);
2032
2033 mg_hcal_res_eta->Draw("APE");
2034 //DrawAxis(mg_hcal_res_eta, leg_hcal_res_eta, etaMin, etaMax, 1, 100, " #eta ", "(HCAL resolution in E)/E (%)", false, true);
2035 DrawAxis(mg_hcal_res_eta, leg_hcal_res_eta, etaMin, etaMax, 0.0, 50, " #eta ", "(HCAL resolution in E)/E (%)", false, false);
2036 leg_hcal_res_eta->Draw();
2037 pave->Draw();
2038
2039 c_hcal_res_eta->Print(pdfOutput, "pdf");
2040 c_hcal_res_eta->Print(figPath + "img_hcal_res_eta.pdf", "pdf");
2041 c_hcal_res_eta->Print(figPath + "img_hcal_res_eta.png", "png");
2042
2043 ////////////////////
2044 // PF - electrons //
2045 ////////////////////
2046
2047 TMultiGraph *mg_pfele_res_e[n_etabins];
2048 TMultiGraph *mg_pfele_res_eta[n_ptbins];
2049
2050 TLegend *leg_pfele_res_e[n_etabins];
2051 TLegend *leg_pfele_res_eta[n_ptbins];
2052
2053 TGraphErrors *gr_pfele_res_e = new TGraphErrors[n_etabins];
2054 TGraphErrors *gr_pfele_res_eta = new TGraphErrors[n_ptbins];
2055 TGraphErrors *gr_trkele_res_e = new TGraphErrors[n_etabins];
2056 TGraphErrors *gr_trkele_res_eeta = new TGraphErrors[n_ptbins];
2057
2058 std::vector<resolPlot> *plots_pfele_res_e = new std::vector<resolPlot>[n_etabins];
2059 std::vector<resolPlot> *plots_pfele_res_eta = new std::vector<resolPlot>[n_ptbins];
2060 std::vector<resolPlot> *plots_trkele_res_e = new std::vector<resolPlot>[n_etabins];
2061 std::vector<resolPlot> *plots_trkele_res_eeta = new std::vector<resolPlot>[n_ptbins];
2062
2063 TCanvas *c_pfele_res_e[n_etabins];
2064 TCanvas *c_pfele_res_eta[n_ptbins];
2065
2066 // loop over eta bins
2067 for(k = 0; k < etaVals.size() - 1; k++)
2068 {
2069 mg_pfele_res_e[k] = new TMultiGraph("", "");
2070 leg_pfele_res_e[k] = new TLegend(0.40, 0.60, 0.75, 0.90);
2071
2072 HistogramsCollection(&plots_pfele_res_e[k], TMath::Log10(ptMin), TMath::Log10(ptMax), "pfele");
2073 GetEres<Electron>(&plots_pfele_res_e[k], branchElectronPF, branchParticleElectron, 11, etaVals.at(k), etaVals.at(k + 1), treeReaderElectron);
2074 gr_pfele_res_e[k] = EresGraph(&plots_pfele_res_e[k]);
2075
2076 HistogramsCollection(&plots_trkele_res_e[k], TMath::Log10(ptMin), TMath::Log10(ptMax), "trkele");
2077 GetEres<Track>(&plots_trkele_res_e[k], branchTrackElectron, branchParticleElectron, 11, etaVals.at(k), etaVals.at(k + 1), treeReaderElectron);
2078 gr_trkele_res_e[k] = EresGraph(&plots_trkele_res_e[k]);
2079
2080 s_etaMin = Form("%.1f", etaVals.at(k));
2081 s_etaMax = Form("%.1f", etaVals.at(k + 1));
2082 s_eta = "e^{ #pm}, " + s_etaMin + " < | #eta | < " + s_etaMax;
2083
2084 leg_pfele_res_e[k]->SetTextFont(132);
2085 leg_pfele_res_e[k]->SetHeader(s_eta);
2086
2087 addResoGraph(mg_pfele_res_e[k], &gr_ecal_res_e[k], leg_pfele_res_e[k], markerStyles.at(0), colors.at(0), "ECAL");
2088 addResoGraph(mg_pfele_res_e[k], &gr_trkele_res_e[k], leg_pfele_res_e[k], markerStyles.at(1), colors.at(1), "Track");
2089 addResoGraph(mg_pfele_res_e[k], &gr_pfele_res_e[k], leg_pfele_res_e[k], markerStyles.at(2), colors.at(2), "Particle-flow");
2090
2091 c_pfele_res_e[k] = new TCanvas("", "", 800, 600);
2092
2093 mg_pfele_res_e[k]->Draw("APE");
2094 //DrawAxis(mg_pfele_res_e[k], leg_pfele_res_e[k], ptMin, ptMax, 0.1, 100, "E [GeV]", "(resolution in E)/E (%)", true, true);
2095 DrawAxis(mg_pfele_res_e[k], leg_pfele_res_e[k], ptMin, ptMax, 0.0, 20, "E [GeV]", "(resolution in E)/E (%)", true, false);
2096 leg_pfele_res_e[k]->Draw();
2097 pave->Draw();
2098
2099 TString s_etarange = "eta_" + s_etaMin + "_" + s_etaMax + "_";
2100
2101 c_pfele_res_e[k]->Print(pdfOutput, "pdf");
2102 c_pfele_res_e[k]->Print(figPath + "img_pfele_res_" + s_etarange + "e.pdf", "pdf");
2103 c_pfele_res_e[k]->Print(figPath + "img_pfele_res_" + s_etarange + "e.png", "png");
2104 }
2105
2106 // loop over eta bins
2107 for(k = 0; k < ptVals.size(); k++)
2108 {
2109
2110 mg_pfele_res_eta[k] = new TMultiGraph("", "");
2111 leg_pfele_res_eta[k] = new TLegend(0.40, 0.60, 0.75, 0.90);
2112
2113 HistogramsCollectionVsEta(&plots_pfele_res_eta[k], etaMin, etaMax, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k), "pfele", 0.0, 2.0);
2114 GetEresVsEta<Electron>(&plots_pfele_res_eta[k], branchElectronPF, branchParticleElectron, 11, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k), treeReaderElectron);
2115 gr_pfele_res_eta[k] = EresGraphVsEta(&plots_pfele_res_eta[k]);
2116
2117 HistogramsCollectionVsEta(&plots_trkele_res_eeta[k], etaMin, etaMax, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k), "trkele", 0.0, 2.0);
2118 GetEresVsEta<Track>(&plots_trkele_res_eeta[k], branchTrackElectron, branchParticleElectron, 11, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k), treeReaderElectron);
2119 gr_trkele_res_eeta[k] = EresGraphVsEta(&plots_trkele_res_eeta[k]);
2120
2121 s_e = Form("e^{ #pm}, E = %.0f GeV", ptVals.at(k));
2122 if(ptVals.at(k) >= 1000.) s_e = Form("e^{ #pm}, E = %.0f TeV", ptVals.at(k) / 1000.);
2123
2124 leg_pfele_res_eta[k]->SetTextFont(132);
2125 leg_pfele_res_eta[k]->SetHeader(s_e);
2126
2127 addResoGraph(mg_pfele_res_eta[k], &gr_ecal_res_eta[k], leg_pfele_res_eta[k], markerStyles.at(0), colors.at(0), "ECAL");
2128 addResoGraph(mg_pfele_res_eta[k], &gr_trkele_res_eeta[k], leg_pfele_res_eta[k], markerStyles.at(1), colors.at(1), "Track");
2129 addResoGraph(mg_pfele_res_eta[k], &gr_pfele_res_eta[k], leg_pfele_res_eta[k], markerStyles.at(2), colors.at(2), "Particle-flow");
2130
2131 c_pfele_res_eta[k] = new TCanvas("", "", 800, 600);
2132
2133 mg_pfele_res_eta[k]->Draw("APE");
2134 //DrawAxis(mg_pfele_res_eta[k], leg_pfele_res_eta[k], etaMin, etaMax, 0.1, 1000, "#eta", "(resolution in E)/E (%)", false, true);
2135 DrawAxis(mg_pfele_res_eta[k], leg_pfele_res_eta[k], etaMin, etaMax, 0.0, 50, "#eta", "(resolution in E)/E (%)", false, false);
2136 leg_pfele_res_eta[k]->Draw();
2137 pave->Draw();
2138
2139 TString s_ptrange = Form("pt_%.0f_", ptVals.at(k));
2140
2141 c_pfele_res_eta[k]->Print(pdfOutput, "pdf");
2142 c_pfele_res_eta[k]->Print(figPath + "img_pfele_res_" + s_ptrange + "eta.pdf", "pdf");
2143 c_pfele_res_eta[k]->Print(figPath + "img_pfele_res_" + s_ptrange + "eta.png", "png");
2144 }
2145
2146 /////////////////
2147 // PF - Pions //
2148 /////////////////
2149
2150 TMultiGraph *mg_pfpi_res_e[n_etabins];
2151 TMultiGraph *mg_pfpi_res_eta[n_ptbins];
2152
2153 TLegend *leg_pfpi_res_e[n_etabins];
2154 TLegend *leg_pfpi_res_eta[n_ptbins];
2155
2156 TGraphErrors *gr_pfpi_res_e = new TGraphErrors[n_etabins];
2157 TGraphErrors *gr_pfpi_res_eta = new TGraphErrors[n_ptbins];
2158
2159 TGraphErrors *gr_trkpi_res_e = new TGraphErrors[n_etabins];
2160 TGraphErrors *gr_trkpi_res_eeta = new TGraphErrors[n_ptbins];
2161
2162 std::vector<resolPlot> *plots_pfpi_res_e = new std::vector<resolPlot>[n_etabins];
2163 std::vector<resolPlot> *plots_pfpi_res_eta = new std::vector<resolPlot>[n_ptbins];
2164 std::vector<resolPlot> *plots_trkpi_res_e = new std::vector<resolPlot>[n_etabins];
2165 std::vector<resolPlot> *plots_trkpi_res_eeta = new std::vector<resolPlot>[n_ptbins];
2166
2167 TCanvas *c_pfpi_res_e[n_etabins];
2168 TCanvas *c_pfpi_res_eta[n_ptbins];
2169
2170 // loop over eta bins
2171 for(k = 0; k < etaVals.size() - 1; k++)
2172 {
2173 mg_pfpi_res_e[k] = new TMultiGraph("", "");
2174 leg_pfpi_res_e[k] = new TLegend(0.40, 0.60, 0.75, 0.90);
2175
2176 HistogramsCollection(&plots_pfpi_res_e[k], TMath::Log10(ptMin), TMath::Log10(ptMax), "pfpi");
2177 GetEres<Track>(&plots_pfpi_res_e[k], branchPion, branchParticlePion, 211, etaVals.at(k), etaVals.at(k + 1), treeReaderPion);
2178 gr_pfpi_res_e[k] = EresGraph(&plots_pfpi_res_e[k]);
2179
2180 HistogramsCollection(&plots_trkpi_res_e[k], TMath::Log10(ptMin), TMath::Log10(ptMax), "trkpi");
2181 GetEres<Track>(&plots_trkpi_res_e[k], branchTrackPion, branchParticlePion, 211, etaVals.at(k), etaVals.at(k + 1), treeReaderPion);
2182 gr_trkpi_res_e[k] = EresGraph(&plots_trkpi_res_e[k]);
2183
2184 s_etaMin = Form("%.1f", etaVals.at(k));
2185 s_etaMax = Form("%.1f", etaVals.at(k + 1));
2186 s_eta = "#pi^{ #pm}, " + s_etaMin + " < | #eta | < " + s_etaMax;
2187
2188 leg_pfpi_res_e[k]->SetTextFont(132);
2189 leg_pfpi_res_e[k]->SetHeader(s_eta);
2190
2191 addResoGraph(mg_pfpi_res_e[k], &gr_hcal_res_e[k], leg_pfpi_res_e[k], markerStyles.at(0), colors.at(0), "HCAL");
2192 addResoGraph(mg_pfpi_res_e[k], &gr_trkpi_res_e[k], leg_pfpi_res_e[k], markerStyles.at(1), colors.at(1), "Track");
2193 addResoGraph(mg_pfpi_res_e[k], &gr_pfpi_res_e[k], leg_pfpi_res_e[k], markerStyles.at(2), colors.at(2), "Particle-flow");
2194
2195 c_pfpi_res_e[k] = new TCanvas("", "", 800, 600);
2196
2197 mg_pfpi_res_e[k]->Draw("APE");
2198 //DrawAxis(mg_pfpi_res_e[k], leg_pfpi_res_e[k], ptMin, ptMax, 0.1, 100, "E [GeV]", "(resolution in E)/E (%)", true, true);
2199 DrawAxis(mg_pfpi_res_e[k], leg_pfpi_res_e[k], ptMin, ptMax, 0.1, 50, "E [GeV]", "(resolution in E)/E (%)", true, false);
2200 leg_pfpi_res_e[k]->Draw();
2201 pave->Draw();
2202
2203 TString s_etarange = "eta_" + s_etaMin + "_" + s_etaMax + "_";
2204
2205 c_pfpi_res_e[k]->Print(pdfOutput, "pdf");
2206 c_pfpi_res_e[k]->Print(figPath + "img_pfpi_res_" + s_etarange + "e.pdf", "pdf");
2207 c_pfpi_res_e[k]->Print(figPath + "img_pfpi_res_" + s_etarange + "e.png", "png");
2208 }
2209
2210 // loop over eta bins
2211 for(k = 0; k < ptVals.size(); k++)
2212 {
2213
2214 mg_pfpi_res_eta[k] = new TMultiGraph("", "");
2215 leg_pfpi_res_eta[k] = new TLegend(0.40, 0.60, 0.75, 0.90);
2216
2217 HistogramsCollectionVsEta(&plots_pfpi_res_eta[k], etaMin, etaMax, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k), "pfpi", 0.0, 2.0);
2218 GetEresVsEta<Track>(&plots_pfpi_res_eta[k], branchPion, branchParticlePion, 211, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k), treeReaderPion);
2219 gr_pfpi_res_eta[k] = EresGraphVsEta(&plots_pfpi_res_eta[k]);
2220
2221 HistogramsCollectionVsEta(&plots_trkpi_res_eeta[k], etaMin, etaMax, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k), "trkpi", 0.0, 2.0);
2222 GetEresVsEta<Track>(&plots_trkpi_res_eeta[k], branchPion, branchParticlePion, 211, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k), treeReaderPion);
2223 gr_trkpi_res_eeta[k] = EresGraphVsEta(&plots_trkpi_res_eeta[k]);
2224
2225 s_e = Form("#pi^{ #pm}, E = %.0f GeV", ptVals.at(k));
2226 if(ptVals.at(k) >= 1000.) s_e = Form("#pi^{ #pm}, E = %.0f TeV", ptVals.at(k) / 1000.);
2227
2228 leg_pfpi_res_eta[k]->SetTextFont(132);
2229 leg_pfpi_res_eta[k]->SetHeader(s_e);
2230
2231 addResoGraph(mg_pfpi_res_eta[k], &gr_hcal_res_eta[k], leg_pfpi_res_eta[k], markerStyles.at(0), colors.at(0), "HCAL");
2232 addResoGraph(mg_pfpi_res_eta[k], &gr_trkpi_res_eeta[k], leg_pfpi_res_eta[k], markerStyles.at(1), colors.at(1), "Track");
2233 addResoGraph(mg_pfpi_res_eta[k], &gr_pfpi_res_eta[k], leg_pfpi_res_eta[k], markerStyles.at(2), colors.at(2), "Particle-flow");
2234
2235 c_pfpi_res_eta[k] = new TCanvas("", "", 800, 600);
2236
2237 mg_pfpi_res_eta[k]->Draw("APE");
2238 //DrawAxis(mg_pfpi_res_eta[k], leg_pfpi_res_eta[k], etaMin, etaMax, 0.1, 1000, "#eta", "(resolution in E)/E (%)", false, true);
2239 DrawAxis(mg_pfpi_res_eta[k], leg_pfpi_res_eta[k], etaMin, etaMax, 0.0, 50, "#eta", "(resolution in E)/E (%)", false, false);
2240 leg_pfpi_res_eta[k]->Draw();
2241 pave->Draw();
2242
2243 TString s_ptrange = Form("pt_%.0f_", ptVals.at(k));
2244
2245 c_pfpi_res_eta[k]->Print(pdfOutput, "pdf");
2246 c_pfpi_res_eta[k]->Print(figPath + "img_pfpi_res_" + s_ptrange + "eta.pdf", "pdf");
2247 c_pfpi_res_eta[k]->Print(figPath + "img_pfpi_res_" + s_ptrange + "eta.png", "png");
2248 }
2249
2250 /////////////////
2251 // PF - jets //
2252 /////////////////
2253
2254 TMultiGraph *mg_pfjet_res_e[n_etabins];
2255 TMultiGraph *mg_pfjet_res_eta[n_ptbins];
2256
2257 TLegend *leg_pfjet_res_e[n_etabins];
2258 TLegend *leg_pfjet_res_eta[n_ptbins];
2259
2260 TGraphErrors *gr_pfjet_res_e = new TGraphErrors[n_etabins];
2261 TGraphErrors *gr_pfjet_res_eta = new TGraphErrors[n_ptbins];
2262
2263 TGraphErrors *gr_cajet_res_e = new TGraphErrors[n_etabins];
2264 TGraphErrors *gr_cajet_res_eta = new TGraphErrors[n_ptbins];
2265
2266 std::vector<resolPlot> *plots_pfjet_res_e = new std::vector<resolPlot>[n_etabins];
2267 std::vector<resolPlot> *plots_pfjet_res_eta = new std::vector<resolPlot>[n_ptbins];
2268 std::vector<resolPlot> *plots_cajet_res_e = new std::vector<resolPlot>[n_etabins];
2269 std::vector<resolPlot> *plots_cajet_res_eta = new std::vector<resolPlot>[n_ptbins];
2270
2271 TCanvas *c_pfjet_res_e[n_etabins];
2272 TCanvas *c_pfjet_res_eta[n_ptbins];
2273
2274 // loop over eta bins
2275 for(k = 0; k < etaVals.size() - 1; k++)
2276 {
2277
2278 mg_pfjet_res_e[k] = new TMultiGraph("", "");
2279 leg_pfjet_res_e[k] = new TLegend(0.40, 0.70, 0.90, 0.90);
2280
2281 HistogramsCollection(&plots_pfjet_res_e[k], TMath::Log10(ptMin), TMath::Log10(ptMax), "pfjet");
2282 HistogramsCollection(&plots_cajet_res_e[k], TMath::Log10(ptMin), TMath::Log10(ptMax), "cajet");
2283
2284 GetJetsEres(&plots_pfjet_res_e[k], branchPFJet, branchGenJet, treeReaderJet, etaVals.at(k), etaVals.at(k + 1));
2285 GetJetsEres(&plots_cajet_res_e[k], branchCaloJet, branchGenJet, treeReaderJet, etaVals.at(k), etaVals.at(k + 1));
2286
2287 gr_pfjet_res_e[k] = EresGraph(&plots_pfjet_res_e[k]);
2288 gr_cajet_res_e[k] = EresGraph(&plots_cajet_res_e[k]);
2289
2290 s_etaMin = Form("%.1f", etaVals.at(k));
2291 s_etaMax = Form("%.1f", etaVals.at(k + 1));
2292 s_eta = "anti-k_{T}, R = 0.4, " + s_etaMin + " < | #eta | < " + s_etaMax;
2293
2294 leg_pfjet_res_e[k]->SetTextFont(132);
2295 leg_pfjet_res_e[k]->SetHeader(s_eta);
2296
2297 addResoGraph(mg_pfjet_res_e[k], &gr_cajet_res_e[k], leg_pfjet_res_e[k], markerStyles.at(0), colors.at(0), "Calorimeter Jets");
2298 addResoGraph(mg_pfjet_res_e[k], &gr_pfjet_res_e[k], leg_pfjet_res_e[k], markerStyles.at(1), colors.at(1), "Particle-flow Jets");
2299
2300 c_pfjet_res_e[k] = new TCanvas("", "", 800, 600);
2301
2302 mg_pfjet_res_e[k]->Draw("APE");
2303 //DrawAxis(mg_pfjet_res_e[k], leg_pfjet_res_e[k], 10, ptMax, 0.5, 100, "E [GeV]", "(resolution in E)/E (%)", true, true);
2304 DrawAxis(mg_pfjet_res_e[k], leg_pfjet_res_e[k], 10, ptMax, 0.0, 30, "E [GeV]", "(resolution in E)/E (%)", true, false);
2305 leg_pfjet_res_e[k]->Draw();
2306 pave->Draw();
2307
2308 TString s_etarange = "eta_" + s_etaMin + "_" + s_etaMax + "_";
2309
2310 c_pfjet_res_e[k]->Print(pdfOutput, "pdf");
2311 c_pfjet_res_e[k]->Print(figPath + "img_pfjet_res_" + s_etarange + "e.pdf", "pdf");
2312 c_pfjet_res_e[k]->Print(figPath + "img_pfjet_res_" + s_etarange + "e.png", "png");
2313 }
2314
2315 // loop over eta bins
2316 for(k = 0; k < ptVals.size(); k++)
2317 {
2318
2319 mg_pfjet_res_eta[k] = new TMultiGraph("", "");
2320 leg_pfjet_res_eta[k] = new TLegend(0.30, 0.70, 0.85, 0.90);
2321
2322 HistogramsCollectionVsEta(&plots_pfjet_res_eta[k], etaMin, etaMax, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k), "pfjet", 0.0, 2.0);
2323 HistogramsCollectionVsEta(&plots_cajet_res_eta[k], etaMin, etaMax, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k), "cajet", 0.0, 2.0);
2324
2325 GetJetsEresVsEta(&plots_pfjet_res_eta[k], branchPFJet, branchGenJet, treeReaderJet, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k));
2326 GetJetsEresVsEta(&plots_cajet_res_eta[k], branchCaloJet, branchGenJet, treeReaderJet, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k));
2327
2328 gr_pfjet_res_eta[k] = EresGraphVsEta(&plots_pfjet_res_eta[k]);
2329 gr_cajet_res_eta[k] = EresGraphVsEta(&plots_cajet_res_eta[k]);
2330
2331 s_e = Form("anti-k_{T}, R = 0.4, jets, E = %.0f GeV", ptVals.at(k));
2332 if(ptVals.at(k) >= 1000.) s_e = Form("anti-k_{T}, R = 0.4, E = %.0f TeV", ptVals.at(k) / 1000.);
2333
2334 leg_pfjet_res_eta[k]->SetTextFont(132);
2335 leg_pfjet_res_eta[k]->SetHeader(s_e);
2336
2337 addResoGraph(mg_pfjet_res_eta[k], &gr_cajet_res_eta[k], leg_pfjet_res_eta[k], markerStyles.at(0), colors.at(0), "Calorimeter Jets");
2338 addResoGraph(mg_pfjet_res_eta[k], &gr_pfjet_res_eta[k], leg_pfjet_res_eta[k], markerStyles.at(1), colors.at(1), "Particle-flow Jets");
2339
2340 c_pfjet_res_eta[k] = new TCanvas("", "", 800, 600);
2341
2342 mg_pfjet_res_eta[k]->Draw("APE");
2343 //DrawAxis(mg_pfjet_res_eta[k], leg_pfjet_res_eta[k], etaMin, etaMax, 0.1, 1000, "#eta", "(resolution in E)/E (%)", false, true);
2344 DrawAxis(mg_pfjet_res_eta[k], leg_pfjet_res_eta[k], etaMin, etaMax, 0.0, 50, "#eta", "(resolution in E)/E (%)", false, false);
2345 leg_pfjet_res_eta[k]->Draw();
2346 pave->Draw();
2347
2348 TString s_ptrange = Form("pt_%.0f_", ptVals.at(k));
2349
2350 c_pfjet_res_eta[k]->Print(pdfOutput, "pdf");
2351 c_pfjet_res_eta[k]->Print(figPath + "img_pfjet_res_" + s_ptrange + "eta.pdf", "pdf");
2352 c_pfjet_res_eta[k]->Print(figPath + "img_pfjet_res_" + s_ptrange + "eta.png", "png");
2353 }
2354
2355 /////////////////////
2356 // PF Missing ET ///
2357 /////////////////////
2358
2359 TMultiGraph *mg_met_res_ht = new TMultiGraph("", "");
2360 TLegend *leg_met_res_ht = new TLegend(0.60, 0.22, 0.90, 0.42);
2361
2362 std::vector<resolPlot> plots_pfmet, plots_camet;
2363
2364 HistogramsCollection(&plots_pfmet, TMath::Log10(ptMin), TMath::Log10(ptMax), "pfMET", -500, 500);
2365 HistogramsCollection(&plots_camet, TMath::Log10(ptMin), TMath::Log10(ptMax), "caMET", -500, 500);
2366
2367 GetMetres(&plots_pfmet, branchGenScalarHT, branchMet, branchPFJet, treeReaderJet);
2368 GetMetres(&plots_camet, branchGenScalarHT, branchCaloMet, branchCaloJet, treeReaderJet);
2369
2370 TGraphErrors gr_pfmet_res_ht = MetResGraph(&plots_pfmet, true);
2371 TGraphErrors gr_camet_res_ht = MetResGraph(&plots_camet, true);
2372
2373 addResoGraph(mg_met_res_ht, &gr_camet_res_ht, leg_met_res_ht, markerStyles.at(0), colors.at(0), "Calorimeter E_{T}^{miss}");
2374 addResoGraph(mg_met_res_ht, &gr_pfmet_res_ht, leg_met_res_ht, markerStyles.at(1), colors.at(1), "Particle-flow E_{T}^{miss}");
2375
2376 TCanvas *c_met_res_ht = new TCanvas("", "", 800, 600);
2377
2378 mg_met_res_ht->Draw("APE");
2379 DrawAxis(mg_met_res_ht, leg_met_res_ht, 1000, 100000, 0.1, 1000, " #sum p_{T} [GeV]", "resolution in E_{x,y}^{miss} [GeV]", true, true);
2380
2381 leg_met_res_ht->Draw();
2382 pave->Draw();
2383 c_met_res_ht->Print(pdfOutput, "pdf");
2384 c_met_res_ht->Print(figPath + "img_met_res_ht.pdf", "pdf");
2385 c_met_res_ht->Print(figPath + "img_met_res_ht.png", "png");
2386
2387 /////////////////////////////////////////
2388 // Electron Reconstruction Efficiency ///
2389 /////////////////////////////////////////
2390
2391 TMultiGraph *mg_recele_eff_pt = new TMultiGraph("", "");
2392 TMultiGraph *mg_recele_eff_eta = new TMultiGraph("", "");
2393
2394 TLegend *leg_recele_eff_pt = new TLegend(0.55, 0.22, 0.90, 0.48);
2395 TLegend *leg_recele_eff_eta = new TLegend(0.55, 0.22, 0.90, 0.48);
2396
2397 TGraphErrors *gr_recele_eff_pt = new TGraphErrors[n_etabins];
2398 TGraphErrors *gr_recele_eff_eta = new TGraphErrors[n_ptbins];
2399 TH1D *h_recele_eff_pt, *h_recele_eff_eta;
2400
2401 // loop over eta bins
2402 for(k = 0; k < etaVals.size() - 1; k++)
2403 {
2404
2405 h_recele_eff_pt = GetEffPt<Electron>(branchElectron, branchParticleElectron, "Electron", 11, ptMin, ptMax, etaVals.at(k), etaVals.at(k + 1), treeReaderElectron);
2406 gr_recele_eff_pt[k] = TGraphErrors(h_recele_eff_pt);
2407
2408 s_etaMin = Form("%.1f", etaVals.at(k));
2409 s_etaMax = Form("%.1f", etaVals.at(k + 1));
2410
2411 s_eta = "e^{ #pm} , " + s_etaMin + " < | #eta | < " + s_etaMax;
2412
2413 gr_recele_eff_pt[k].SetName("recEff_" + s_etaMin + "_" + s_etaMax);
2414
2415 addResoGraph(mg_recele_eff_pt, &gr_recele_eff_pt[k], leg_recele_eff_pt, markerStyles.at(k), colors.at(k), s_eta);
2416 }
2417
2418 // loop over pt
2419 for(k = 0; k < ptVals.size(); k++)
2420 {
2421 h_recele_eff_eta = GetEffEta<Electron>(branchElectron, branchParticleElectron, "Electron", 11, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k), etaMin, etaMax, treeReaderElectron);
2422 gr_recele_eff_eta[k] = TGraphErrors(h_recele_eff_eta);
2423
2424 s_pt = Form("e^{ #pm} , p_{T} = %.0f GeV", ptVals.at(k));
2425 if(ptVals.at(k) >= 1000.) s_pt = Form("e^{ #pm} , p_{T} = %.0f TeV", ptVals.at(k) / 1000.);
2426
2427 addResoGraph(mg_recele_eff_eta, &gr_recele_eff_eta[k], leg_recele_eff_eta, markerStyles.at(k), colors.at(k), s_pt);
2428 }
2429
2430 TCanvas *c_recele_eff_pt = new TCanvas("", "", 800, 600);
2431
2432 mg_recele_eff_pt->Draw("APE");
2433 DrawAxis(mg_recele_eff_pt, leg_recele_eff_pt, ptMin, ptMax, 0.0, 100, "p_{T} [GeV]", "reconstruction efficiency (%)", true, false);
2434 leg_recele_eff_pt->Draw();
2435 pave->Draw();
2436
2437 c_recele_eff_pt->Print(pdfOutput, "pdf");
2438 c_recele_eff_pt->Print(figPath + "img_recele_eff_pt.pdf", "pdf");
2439 c_recele_eff_pt->Print(figPath + "img_recele_eff_pt.png", "png");
2440
2441 TCanvas *c_recele_eff_eta = new TCanvas("", "", 800, 600);
2442
2443 mg_recele_eff_eta->Draw("APE");
2444 DrawAxis(mg_recele_eff_eta, leg_recele_eff_eta, etaMin, etaMax, 0.0, 100, " #eta ", "reconstruction efficiency (%)", false, false);
2445 leg_recele_eff_eta->Draw();
2446 pave->Draw();
2447
2448 c_recele_eff_eta->Print(pdfOutput, "pdf");
2449 c_recele_eff_eta->Print(figPath + "img_recele_eff_eta.pdf", "pdf");
2450 c_recele_eff_eta->Print(figPath + "img_recele_eff_eta.png", "png");
2451
2452 /////////////////////////////////////////
2453 // Muon Reconstruction Efficiency ///
2454 /////////////////////////////////////////
2455
2456 TMultiGraph *mg_recmu_eff_pt = new TMultiGraph("", "");
2457 TMultiGraph *mg_recmu_eff_eta = new TMultiGraph("", "");
2458
2459 TLegend *leg_recmu_eff_pt = new TLegend(0.55, 0.22, 0.90, 0.48);
2460 TLegend *leg_recmu_eff_eta = new TLegend(0.55, 0.22, 0.90, 0.48);
2461
2462 TGraphErrors *gr_recmu_eff_pt = new TGraphErrors[n_etabins];
2463 TGraphErrors *gr_recmu_eff_eta = new TGraphErrors[n_ptbins];
2464 TH1D *h_recmu_eff_pt, *h_recmu_eff_eta;
2465
2466 // loop over eta bins
2467 for(k = 0; k < etaVals.size() - 1; k++)
2468 {
2469
2470 h_recmu_eff_pt = GetEffPt<Muon>(branchMuon, branchParticleMuon, "muon", 13, ptMin, ptMax, etaVals.at(k), etaVals.at(k + 1), treeReaderMuon);
2471 gr_recmu_eff_pt[k] = TGraphErrors(h_recmu_eff_pt);
2472
2473 s_etaMin = Form("%.1f", etaVals.at(k));
2474 s_etaMax = Form("%.1f", etaVals.at(k + 1));
2475
2476 s_eta = "#mu^{ #pm} , " + s_etaMin + " < | #eta | < " + s_etaMax;
2477
2478 gr_recmu_eff_pt[k].SetName("recEff_" + s_etaMin + "_" + s_etaMax);
2479
2480 addResoGraph(mg_recmu_eff_pt, &gr_recmu_eff_pt[k], leg_recmu_eff_pt, markerStyles.at(k), colors.at(k), s_eta);
2481 }
2482
2483 // loop over pt
2484 for(k = 0; k < ptVals.size(); k++)
2485 {
2486 h_recmu_eff_eta = GetEffEta<Muon>(branchMuon, branchParticleMuon, "muon", 13, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k), etaMin, etaMax, treeReaderMuon);
2487 gr_recmu_eff_eta[k] = TGraphErrors(h_recmu_eff_eta);
2488
2489 s_pt = Form("#mu^{ #pm} , p_{T} = %.0f GeV", ptVals.at(k));
2490 if(ptVals.at(k) >= 1000.) s_pt = Form("#mu^{ #pm} , p_{T} = %.0f TeV", ptVals.at(k) / 1000.);
2491
2492 addResoGraph(mg_recmu_eff_eta, &gr_recmu_eff_eta[k], leg_recmu_eff_eta, markerStyles.at(k), colors.at(k), s_pt);
2493 }
2494
2495 TCanvas *c_recmu_eff_pt = new TCanvas("", "", 800, 600);
2496
2497 mg_recmu_eff_pt->Draw("APE");
2498 DrawAxis(mg_recmu_eff_pt, leg_recmu_eff_pt, ptMin, ptMax, 0.0, 100, "p_{T} [GeV]", "reconstruction efficiency (%)", true, false);
2499 leg_recmu_eff_pt->Draw();
2500 pave->Draw();
2501
2502 c_recmu_eff_pt->Print(pdfOutput, "pdf");
2503 c_recmu_eff_pt->Print(figPath + "img_recmu_eff_pt.pdf", "pdf");
2504 c_recmu_eff_pt->Print(figPath + "img_recmu_eff_pt.png", "png");
2505
2506 TCanvas *c_recmu_eff_eta = new TCanvas("", "", 800, 600);
2507
2508 mg_recmu_eff_eta->Draw("APE");
2509 DrawAxis(mg_recmu_eff_eta, leg_recmu_eff_eta, etaMin, etaMax, 0.0, 100, " #eta ", "reconstruction efficiency (%)", false, false);
2510 leg_recmu_eff_eta->Draw();
2511 pave->Draw();
2512
2513 c_recmu_eff_eta->Print(pdfOutput, "pdf");
2514 c_recmu_eff_eta->Print(figPath + "img_recmu_eff_eta.pdf", "pdf");
2515 c_recmu_eff_eta->Print(figPath + "img_recmu_eff_eta.png", "png");
2516
2517 /////////////////////////////////////////
2518 // Photon Reconstruction Efficiency ///
2519 /////////////////////////////////////////
2520
2521 TMultiGraph *mg_recpho_eff_pt = new TMultiGraph("", "");
2522 TMultiGraph *mg_recpho_eff_eta = new TMultiGraph("", "");
2523
2524 TLegend *leg_recpho_eff_pt = new TLegend(0.55, 0.22, 0.90, 0.48);
2525 TLegend *leg_recpho_eff_eta = new TLegend(0.55, 0.22, 0.90, 0.48);
2526
2527 TGraphErrors *gr_recpho_eff_pt = new TGraphErrors[n_etabins];
2528 TGraphErrors *gr_recpho_eff_eta = new TGraphErrors[n_ptbins];
2529 TH1D *h_recpho_eff_pt, *h_recpho_eff_eta;
2530
2531 // loop over eta bins
2532 for(k = 0; k < etaVals.size() - 1; k++)
2533 {
2534
2535 h_recpho_eff_pt = GetEffPt<Photon>(branchPhoton, branchParticlePhoton, "Photon", 22, ptMin, ptMax, etaVals.at(k), etaVals.at(k + 1), treeReaderPhoton);
2536 gr_recpho_eff_pt[k] = TGraphErrors(h_recpho_eff_pt);
2537
2538 s_etaMin = Form("%.1f", etaVals.at(k));
2539 s_etaMax = Form("%.1f", etaVals.at(k + 1));
2540
2541 s_eta = "#gamma , " + s_etaMin + " < | #eta | < " + s_etaMax;
2542
2543 gr_recpho_eff_pt[k].SetName("recEff_" + s_etaMin + "_" + s_etaMax);
2544
2545 addResoGraph(mg_recpho_eff_pt, &gr_recpho_eff_pt[k], leg_recpho_eff_pt, markerStyles.at(k), colors.at(k), s_eta);
2546 }
2547
2548 // loop over pt
2549 for(k = 0; k < ptVals.size(); k++)
2550 {
2551 h_recpho_eff_eta = GetEffEta<Photon>(branchPhoton, branchParticlePhoton, "Photon", 22, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k), etaMin, etaMax, treeReaderPhoton);
2552 gr_recpho_eff_eta[k] = TGraphErrors(h_recpho_eff_eta);
2553
2554 s_pt = Form("#gamma , p_{T} = %.0f GeV", ptVals.at(k));
2555 if(ptVals.at(k) >= 1000.) s_pt = Form("#gamma , p_{T} = %.0f TeV", ptVals.at(k) / 1000.);
2556
2557 addResoGraph(mg_recpho_eff_eta, &gr_recpho_eff_eta[k], leg_recpho_eff_eta, markerStyles.at(k), colors.at(k), s_pt);
2558 }
2559
2560 TCanvas *c_recpho_eff_pt = new TCanvas("", "", 800, 600);
2561
2562 mg_recpho_eff_pt->Draw("APE");
2563 DrawAxis(mg_recpho_eff_pt, leg_recpho_eff_pt, ptMin, ptMax, 0.0, 100, "p_{T} [GeV]", "reconstruction efficiency (%)", true, false);
2564 leg_recpho_eff_pt->Draw();
2565 pave->Draw();
2566
2567 c_recpho_eff_pt->Print(pdfOutput, "pdf");
2568 c_recpho_eff_pt->Print(figPath + "img_recpho_eff_pt.pdf", "pdf");
2569 c_recpho_eff_pt->Print(figPath + "img_recpho_eff_pt.png", "png");
2570
2571 TCanvas *c_recpho_eff_eta = new TCanvas("", "", 800, 600);
2572
2573 mg_recpho_eff_eta->Draw("APE");
2574 DrawAxis(mg_recpho_eff_eta, leg_recpho_eff_eta, etaMin, etaMax, 0.0, 100, " #eta ", "reconstruction efficiency (%)", false, false);
2575 leg_recpho_eff_eta->Draw();
2576 pave->Draw();
2577
2578 c_recpho_eff_eta->Print(pdfOutput, "pdf");
2579 c_recpho_eff_eta->Print(figPath + "img_recpho_eff_eta.pdf", "pdf");
2580 c_recpho_eff_eta->Print(figPath + "img_recpho_eff_eta.png", "png");
2581
2582 /////////////////////////////////////////
2583 // B-jets Efficiency/ mistag rates ///
2584 /////////////////////////////////////////
2585
2586 TMultiGraph *mg_recbjet_eff_pt = new TMultiGraph("", "");
2587 TMultiGraph *mg_recbjet_eff_eta = new TMultiGraph("", "");
2588
2589 TLegend *leg_recbjet_eff_pt = new TLegend(0.50, 0.22, 0.90, 0.48);
2590 TLegend *leg_recbjet_eff_eta = new TLegend(0.50, 0.22, 0.90, 0.48);
2591
2592 TGraphErrors *gr_recbjet_eff_pt = new TGraphErrors[n_etabins];
2593 TGraphErrors *gr_recbjet_eff_eta = new TGraphErrors[n_ptbins];
2594 TH1D *h_recbjet_eff_pt, *h_recbjet_eff_eta;
2595
2596 // loop over eta bins
2597 for(k = 0; k < etaVals.size() - 1; k++)
2598 {
2599
2600 h_recbjet_eff_pt = GetJetEffPt<Jet>(branchPFBJet, "BJet", 5, ptMin, ptMax, etaVals.at(k), etaVals.at(k + 1), treeReaderBJet);
2601 //h_recbjet_eff_pt = GetEffPt<Jet>(branchPFBJet, branchParticleBJet, "BJet", 5, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderBJet);
2602 gr_recbjet_eff_pt[k] = TGraphErrors(h_recbjet_eff_pt);
2603
2604 s_etaMin = Form("%.1f", etaVals.at(k));
2605 s_etaMax = Form("%.1f", etaVals.at(k + 1));
2606
2607 s_eta = "b-jet , " + s_etaMin + " < | #eta | < " + s_etaMax;
2608
2609 gr_recbjet_eff_pt[k].SetName("recEff_" + s_etaMin + "_" + s_etaMax);
2610
2611 addResoGraph(mg_recbjet_eff_pt, &gr_recbjet_eff_pt[k], leg_recbjet_eff_pt, markerStyles.at(k), colors.at(k), s_eta);
2612 }
2613
2614 // loop over pt
2615 for(k = 0; k < ptVals.size(); k++)
2616 {
2617 h_recbjet_eff_eta = GetJetEffEta<Jet>(branchPFBJet, "BJet", 5, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k), etaMin, etaMax, treeReaderBJet);
2618 //h_recbjet_eff_eta = GetEffEta<Jet>(branchPFBJet, branchParticleBJet, "BJet", 5, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderBJet);
2619 gr_recbjet_eff_eta[k] = TGraphErrors(h_recbjet_eff_eta);
2620
2621 s_pt = Form("b-jet , p_{T} = %.0f GeV", ptVals.at(k));
2622 if(ptVals.at(k) >= 1000.) s_pt = Form("b-jet , p_{T} = %.0f TeV", ptVals.at(k) / 1000.);
2623
2624 addResoGraph(mg_recbjet_eff_eta, &gr_recbjet_eff_eta[k], leg_recbjet_eff_eta, markerStyles.at(k), colors.at(k), s_pt);
2625 }
2626
2627 TCanvas *c_recbjet_eff_pt = new TCanvas("", "", 800, 600);
2628
2629 mg_recbjet_eff_pt->Draw("APE");
2630 DrawAxis(mg_recbjet_eff_pt, leg_recbjet_eff_pt, ptMin, ptMax, 0.0, 100, "p_{T} [GeV]", "b - tag efficiency (%)", true, false);
2631 leg_recbjet_eff_pt->Draw();
2632 pave->Draw();
2633
2634 c_recbjet_eff_pt->Print(pdfOutput, "pdf");
2635 c_recbjet_eff_pt->Print(figPath + "img_recbjet_eff_pt.pdf", "pdf");
2636 c_recbjet_eff_pt->Print(figPath + "img_recbjet_eff_pt.png", "png");
2637
2638 TCanvas *c_recbjet_eff_eta = new TCanvas("", "", 800, 600);
2639
2640 mg_recbjet_eff_eta->Draw("APE");
2641 DrawAxis(mg_recbjet_eff_eta, leg_recbjet_eff_eta, etaMin, etaMax, 0.0, 100, " #eta ", "b - tag efficiency (%)", false, false);
2642 leg_recbjet_eff_eta->Draw();
2643 pave->Draw();
2644
2645 c_recbjet_eff_eta->Print(pdfOutput, "pdf");
2646 c_recbjet_eff_eta->Print(figPath + "img_recbjet_eff_eta.pdf", "pdf");
2647 c_recbjet_eff_eta->Print(figPath + "img_recbjet_eff_eta.png", "png");
2648
2649 // ------ c - mistag ------
2650
2651 TMultiGraph *mg_recbjet_cmis_pt = new TMultiGraph("", "");
2652 TMultiGraph *mg_recbjet_cmis_eta = new TMultiGraph("", "");
2653
2654 TLegend *leg_recbjet_cmis_pt = new TLegend(0.50, 0.64, 0.90, 0.90);
2655 TLegend *leg_recbjet_cmis_eta = new TLegend(0.50, 0.64, 0.90, 0.90);
2656
2657 TGraphErrors *gr_recbjet_cmis_pt = new TGraphErrors[n_etabins];
2658 TGraphErrors *gr_recbjet_cmis_eta = new TGraphErrors[n_ptbins];
2659 TH1D *h_recbjet_cmis_pt, *h_recbjet_cmis_eta;
2660
2661 // loop over eta bins
2662 for(k = 0; k < etaVals.size() - 1; k++)
2663 {
2664
2665 h_recbjet_cmis_pt = GetJetEffPt<Jet>(branchPFCJet, "CJet", 4, ptMin, ptMax, etaVals.at(k), etaVals.at(k + 1), treeReaderCJet);
2666 //h_recbjet_cmis_pt = GetEffPt<Jet>(branchPFCJet, branchParticleCJet, "CJet", 4, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderCJet);
2667 gr_recbjet_cmis_pt[k] = TGraphErrors(h_recbjet_cmis_pt);
2668
2669 s_etaMin = Form("%.1f", etaVals.at(k));
2670 s_etaMax = Form("%.1f", etaVals.at(k + 1));
2671
2672 s_eta = "c-jet , " + s_etaMin + " < | #eta | < " + s_etaMax;
2673
2674 gr_recbjet_cmis_pt[k].SetName("recEff_" + s_etaMin + "_" + s_etaMax);
2675
2676 addResoGraph(mg_recbjet_cmis_pt, &gr_recbjet_cmis_pt[k], leg_recbjet_cmis_pt, markerStyles.at(k), colors.at(k), s_eta);
2677 }
2678
2679 // loop over pt
2680 for(k = 0; k < ptVals.size(); k++)
2681 {
2682 h_recbjet_cmis_eta = GetJetEffEta<Jet>(branchPFCJet, "CJet", 4, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k), etaMin, etaMax, treeReaderCJet);
2683 //h_recbjet_cmis_eta = GetEffEta<Jet>(branchPFCJet, branchParticleCJet, "CJet", 4, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderCJet);
2684 gr_recbjet_cmis_eta[k] = TGraphErrors(h_recbjet_cmis_eta);
2685
2686 s_pt = Form("c-jet , p_{T} = %.0f GeV", ptVals.at(k));
2687 if(ptVals.at(k) >= 1000.) s_pt = Form("c-jet , p_{T} = %.0f TeV", ptVals.at(k) / 1000.);
2688
2689 addResoGraph(mg_recbjet_cmis_eta, &gr_recbjet_cmis_eta[k], leg_recbjet_cmis_eta, markerStyles.at(k), colors.at(k), s_pt);
2690 }
2691
2692 TCanvas *c_recbjet_cmis_pt = new TCanvas("", "", 800, 600);
2693
2694 mg_recbjet_cmis_pt->Draw("APE");
2695 DrawAxis(mg_recbjet_cmis_pt, leg_recbjet_cmis_pt, ptMin, ptMax, 0.0, 20, "p_{T} [GeV]", "c - mistag rate (%)", true, false);
2696 leg_recbjet_cmis_pt->Draw();
2697 pave->Draw();
2698
2699 c_recbjet_cmis_pt->Print(pdfOutput, "pdf");
2700 c_recbjet_cmis_pt->Print(figPath + "img_recbjet_cmis_pt.pdf", "pdf");
2701 c_recbjet_cmis_pt->Print(figPath + "img_recbjet_cmis_pt.png", "png");
2702
2703 TCanvas *c_recbjet_cmis_eta = new TCanvas("", "", 800, 600);
2704
2705 mg_recbjet_cmis_eta->Draw("APE");
2706 DrawAxis(mg_recbjet_cmis_eta, leg_recbjet_cmis_eta, etaMin, etaMax, 0.0, 20, " #eta ", "c - mistag rate (%)", false, false);
2707 leg_recbjet_cmis_eta->Draw();
2708 pave->Draw();
2709
2710 c_recbjet_cmis_eta->Print(pdfOutput, "pdf");
2711 c_recbjet_cmis_eta->Print(figPath + "img_recbjet_cmis_eta.pdf", "pdf");
2712 c_recbjet_cmis_eta->Print(figPath + "img_recbjet_cmis_eta.png", "png");
2713
2714 // ------ light - mistag ------
2715
2716 TMultiGraph *mg_recbjet_lmis_pt = new TMultiGraph("", "");
2717 TMultiGraph *mg_recbjet_lmis_eta = new TMultiGraph("", "");
2718
2719 TLegend *leg_recbjet_lmis_pt = new TLegend(0.50, 0.64, 0.90, 0.90);
2720 TLegend *leg_recbjet_lmis_eta = new TLegend(0.50, 0.64, 0.90, 0.90);
2721
2722 TGraphErrors *gr_recbjet_lmis_pt = new TGraphErrors[n_etabins];
2723 TGraphErrors *gr_recbjet_lmis_eta = new TGraphErrors[n_ptbins];
2724 TH1D *h_recbjet_lmis_pt, *h_recbjet_lmis_eta;
2725
2726 // loop over eta bins
2727 for(k = 0; k < etaVals.size() - 1; k++)
2728 {
2729
2730 h_recbjet_lmis_pt = GetJetEffPt<Jet>(branchJet, "Jet", 1, ptMin, ptMax, etaVals.at(k), etaVals.at(k + 1), treeReaderJet);
2731 //h_recbjet_lmis_pt = GetEffPt<Jet>(branchJet, branchParticleJet, "Jet", 1, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderJet);
2732 gr_recbjet_lmis_pt[k] = TGraphErrors(h_recbjet_lmis_pt);
2733
2734 s_etaMin = Form("%.1f", etaVals.at(k));
2735 s_etaMax = Form("%.1f", etaVals.at(k + 1));
2736
2737 s_eta = "uds-jet , " + s_etaMin + " < | #eta | < " + s_etaMax;
2738
2739 gr_recbjet_lmis_pt[k].SetName("recEff_" + s_etaMin + "_" + s_etaMax);
2740
2741 addResoGraph(mg_recbjet_lmis_pt, &gr_recbjet_lmis_pt[k], leg_recbjet_lmis_pt, markerStyles.at(k), colors.at(k), s_eta);
2742 }
2743
2744 // loop over pt
2745 for(k = 0; k < ptVals.size(); k++)
2746 {
2747 h_recbjet_lmis_eta = GetJetEffEta<Jet>(branchJet, "Jet", 1, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k), etaMin, etaMax, treeReaderJet);
2748 //h_recbjet_lmis_eta = GetEffEta<Jet>(branchJet, branchParticleJet, "Jet", 1, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderJet);
2749 gr_recbjet_lmis_eta[k] = TGraphErrors(h_recbjet_lmis_eta);
2750
2751 s_pt = Form("uds-jet , p_{T} = %.0f GeV", ptVals.at(k));
2752 if(ptVals.at(k) >= 1000.) s_pt = Form("uds-jet , p_{T} = %.0f TeV", ptVals.at(k) / 1000.);
2753
2754 addResoGraph(mg_recbjet_lmis_eta, &gr_recbjet_lmis_eta[k], leg_recbjet_lmis_eta, markerStyles.at(k), colors.at(k), s_pt);
2755 }
2756
2757 TCanvas *c_recbjet_lmis_pt = new TCanvas("", "", 800, 600);
2758
2759 mg_recbjet_lmis_pt->Draw("APE");
2760
2761 DrawAxis(mg_recbjet_lmis_pt, leg_recbjet_lmis_pt, ptMin, ptMax, 0.0, 1.0, "p_{T} [GeV]", "light - mistag rate (%)", true, false);
2762
2763 leg_recbjet_lmis_pt->Draw();
2764 pave->Draw();
2765
2766 c_recbjet_lmis_pt->Print(pdfOutput, "pdf");
2767 c_recbjet_lmis_pt->Print(figPath + "img_recbjet_lmis_pt.pdf", "pdf");
2768 c_recbjet_lmis_pt->Print(figPath + "img_recbjet_lmis_pt.png", "png");
2769
2770 TCanvas *c_recbjet_lmis_eta = new TCanvas("", "", 800, 600);
2771
2772 mg_recbjet_lmis_eta->Draw("APE");
2773 DrawAxis(mg_recbjet_lmis_eta, leg_recbjet_lmis_eta, etaMin, etaMax, 0.0, 1.0, " #eta ", "light - mistag rate (%)", false, false);
2774 leg_recbjet_lmis_eta->Draw();
2775 pave->Draw();
2776
2777 c_recbjet_lmis_eta->Print(pdfOutput, "pdf");
2778 c_recbjet_lmis_eta->Print(figPath + "img_recbjet_lmis_eta.pdf", "pdf");
2779 c_recbjet_lmis_eta->Print(figPath + "img_recbjet_lmis_eta.png", "png");
2780
2781 ///////////////////////////////////////////
2782 // tau-jets Efficiency/ mistag rates ///
2783 ///////////////////////////////////////////
2784
2785 TMultiGraph *mg_rectaujet_eff_pt = new TMultiGraph("", "");
2786 TMultiGraph *mg_rectaujet_eff_eta = new TMultiGraph("", "");
2787
2788 TLegend *leg_rectaujet_eff_pt = new TLegend(0.50, 0.22, 0.90, 0.48);
2789 TLegend *leg_rectaujet_eff_eta = new TLegend(0.50, 0.22, 0.90, 0.48);
2790
2791 TGraphErrors *gr_rectaujet_eff_pt = new TGraphErrors[n_etabins];
2792 TGraphErrors *gr_rectaujet_eff_eta = new TGraphErrors[n_ptbins];
2793 TH1D *h_rectaujet_eff_pt, *h_rectaujet_eff_eta;
2794
2795 // loop over eta bins
2796 for(k = 0; k < etaVals.size() - 1; k++)
2797 {
2798
2799 h_rectaujet_eff_pt = GetTauEffPt<Jet>(branchPFTauJet, branchParticleTauJet, "TauJet", 15, ptMin, ptMax, etaVals.at(k), etaVals.at(k + 1), treeReaderTauJet);
2800 gr_rectaujet_eff_pt[k] = TGraphErrors(h_rectaujet_eff_pt);
2801
2802 s_etaMin = Form("%.1f", etaVals.at(k));
2803 s_etaMax = Form("%.1f", etaVals.at(k + 1));
2804
2805 s_eta = "#tau-jet , " + s_etaMin + " < | #eta | < " + s_etaMax;
2806
2807 gr_rectaujet_eff_pt[k].SetName("recEff_" + s_etaMin + "_" + s_etaMax);
2808
2809 addResoGraph(mg_rectaujet_eff_pt, &gr_rectaujet_eff_pt[k], leg_rectaujet_eff_pt, markerStyles.at(k), colors.at(k), s_eta);
2810 }
2811
2812 // loop over pt
2813 for(k = 0; k < ptVals.size(); k++)
2814 {
2815 h_rectaujet_eff_eta = GetTauEffEta<Jet>(branchPFTauJet, branchParticleTauJet, "TauJet", 15, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k), etaMin, etaMax, treeReaderTauJet);
2816 gr_rectaujet_eff_eta[k] = TGraphErrors(h_rectaujet_eff_eta);
2817
2818 s_pt = Form("#tau-jet , p_{T} = %.0f GeV", ptVals.at(k));
2819 if(ptVals.at(k) >= 1000.) s_pt = Form("#tau-jet , p_{T} = %.0f TeV", ptVals.at(k) / 1000.);
2820
2821 addResoGraph(mg_rectaujet_eff_eta, &gr_rectaujet_eff_eta[k], leg_rectaujet_eff_eta, markerStyles.at(k), colors.at(k), s_pt);
2822 }
2823
2824 TCanvas *c_rectaujet_eff_pt = new TCanvas("", "", 800, 600);
2825
2826 mg_rectaujet_eff_pt->Draw("APE");
2827 DrawAxis(mg_rectaujet_eff_pt, leg_rectaujet_eff_pt, ptMin, ptMax, 0.0, 100, "p_{T} [GeV]", "#tau - tag efficiency (%)", true, false);
2828 leg_rectaujet_eff_pt->Draw();
2829 pave->Draw();
2830
2831 c_rectaujet_eff_pt->Print(pdfOutput, "pdf");
2832 c_rectaujet_eff_pt->Print(figPath + "img_rectaujet_eff_pt.pdf", "pdf");
2833 c_rectaujet_eff_pt->Print(figPath + "img_rectaujet_eff_pt.png", "png");
2834
2835 TCanvas *c_rectaujet_eff_eta = new TCanvas("", "", 800, 600);
2836
2837 mg_rectaujet_eff_eta->Draw("APE");
2838 DrawAxis(mg_rectaujet_eff_eta, leg_rectaujet_eff_eta, etaMin, etaMax, 0.0, 100., " #eta ", "#tau - tag efficiency (%)", false, false);
2839 leg_rectaujet_eff_eta->Draw();
2840 pave->Draw();
2841
2842 c_rectaujet_eff_eta->Print(pdfOutput, "pdf");
2843 c_rectaujet_eff_eta->Print(figPath + "img_rectaujet_eff_eta.pdf", "pdf");
2844 c_rectaujet_eff_eta->Print(figPath + "img_rectaujet_eff_eta.png", "png");
2845
2846 //--------------- tau mistag rate ----------
2847
2848 TMultiGraph *mg_rectaujet_mis_pt = new TMultiGraph("", "");
2849 TMultiGraph *mg_rectaujet_mis_eta = new TMultiGraph("", "");
2850
2851 TLegend *leg_rectaujet_mis_pt = new TLegend(0.50, 0.64, 0.90, 0.90);
2852 TLegend *leg_rectaujet_mis_eta = new TLegend(0.50, 0.64, 0.90, 0.90);
2853
2854 TGraphErrors *gr_rectaujet_mis_pt = new TGraphErrors[n_etabins];
2855 TGraphErrors *gr_rectaujet_mis_eta = new TGraphErrors[n_ptbins];
2856 TH1D *h_rectaujet_mis_pt, *h_rectaujet_mis_eta;
2857
2858 // loop over eta bins
2859 for(k = 0; k < etaVals.size() - 1; k++)
2860 {
2861
2862 h_rectaujet_mis_pt = GetTauEffPt<Jet>(branchJet, branchParticleJet, "TauJet", 1, ptMin, ptMax, etaVals.at(k), etaVals.at(k + 1), treeReaderJet);
2863 gr_rectaujet_mis_pt[k] = TGraphErrors(h_rectaujet_mis_pt);
2864
2865 s_etaMin = Form("%.1f", etaVals.at(k));
2866 s_etaMax = Form("%.1f", etaVals.at(k + 1));
2867
2868 s_eta = "uds-jet , " + s_etaMin + " < | #eta | < " + s_etaMax;
2869
2870 gr_rectaujet_mis_pt[k].SetName("recEff_" + s_etaMin + "_" + s_etaMax);
2871
2872 addResoGraph(mg_rectaujet_mis_pt, &gr_rectaujet_mis_pt[k], leg_rectaujet_mis_pt, markerStyles.at(k), colors.at(k), s_eta);
2873 }
2874
2875 // loop over pt
2876 for(k = 0; k < ptVals.size(); k++)
2877 {
2878 h_rectaujet_mis_eta = GetTauEffEta<Jet>(branchJet, branchParticleJet, "TauJet", 1, 0.5 * ptVals.at(k), 2.0 * ptVals.at(k), etaMin, etaMax, treeReaderJet);
2879 gr_rectaujet_mis_eta[k] = TGraphErrors(h_rectaujet_mis_eta);
2880
2881 s_pt = Form("uds-jet , p_{T} = %.0f GeV", ptVals.at(k));
2882 if(ptVals.at(k) >= 1000.) s_pt = Form("uds-jet , p_{T} = %.0f TeV", ptVals.at(k) / 1000.);
2883
2884 addResoGraph(mg_rectaujet_mis_eta, &gr_rectaujet_mis_eta[k], leg_rectaujet_mis_eta, markerStyles.at(k), colors.at(k), s_pt);
2885 }
2886
2887 TCanvas *c_rectaujet_mis_pt = new TCanvas("", "", 800, 600);
2888
2889 mg_rectaujet_mis_pt->Draw("APE");
2890 DrawAxis(mg_rectaujet_mis_pt, leg_rectaujet_mis_pt, ptMin, ptMax, 0.0, 5., "p_{T} [GeV]", "#tau - mistag(%)", true, false);
2891 leg_rectaujet_mis_pt->Draw();
2892 pave->Draw();
2893
2894 c_rectaujet_mis_pt->Print(pdfOutput, "pdf");
2895 c_rectaujet_mis_pt->Print(figPath + "img_rectaujet_mis_pt.pdf", "pdf");
2896 c_rectaujet_mis_pt->Print(figPath + "img_rectaujet_mis_pt.png", "png");
2897
2898 TCanvas *c_rectaujet_mis_eta = new TCanvas("", "", 800, 600);
2899
2900 mg_rectaujet_mis_eta->Draw("APE");
2901 DrawAxis(mg_rectaujet_mis_eta, leg_rectaujet_mis_eta, etaMin, etaMax, 0.0, 5., " #eta ", "#tau - mistag (%)", false, false);
2902 leg_rectaujet_mis_eta->Draw();
2903 pave->Draw();
2904
2905 c_rectaujet_mis_eta->Print(pdfOutput + ")", "pdf");
2906 c_rectaujet_mis_eta->Print(figPath + "img_rectaujet_mis_eta.pdf", "pdf");
2907 c_rectaujet_mis_eta->Print(figPath + "img_rectaujet_mis_eta.png", "png");
2908
2909 // ---- store resolution histograms in the output (for leave efficiencies out) ---
2910
2911 TFile *fout = new TFile(outputFile, "recreate");
2912
2913 for(int bin = 0; bin < Nbins; bin++)
2914 {
2915
2916 for(k = 0; k < etaVals.size() - 1; k++)
2917 {
2918 plots_trkpi_res_pt[k].at(bin).resolHist->Write();
2919 plots_trkele_res_pt[k].at(bin).resolHist->Write();
2920 plots_trkmu_res_pt[k].at(bin).resolHist->Write();
2921 plots_ecal_res_e[k].at(bin).resolHist->Write();
2922 plots_hcal_res_e[k].at(bin).resolHist->Write();
2923 plots_pfele_res_e[k].at(bin).resolHist->Write();
2924 plots_pfpi_res_e[k].at(bin).resolHist->Write();
2925 plots_pfjet_res_e[k].at(bin).resolHist->Write();
2926 plots_cajet_res_e[k].at(bin).resolHist->Write();
2927 }
2928
2929 for(k = 0; k < ptVals.size(); k++)
2930 {
2931 plots_trkpi_res_eta[k].at(bin).resolHist->Write();
2932 plots_trkele_res_eta[k].at(bin).resolHist->Write();
2933 plots_trkmu_res_eta[k].at(bin).resolHist->Write();
2934 plots_ecal_res_eta[k].at(bin).resolHist->Write();
2935 plots_hcal_res_eta[k].at(bin).resolHist->Write();
2936 plots_pfele_res_eta[k].at(bin).resolHist->Write();
2937 plots_pfpi_res_eta[k].at(bin).resolHist->Write();
2938 plots_pfjet_res_eta[k].at(bin).resolHist->Write();
2939 plots_cajet_res_eta[k].at(bin).resolHist->Write();
2940 }
2941
2942 plots_pfmet.at(bin).resolHist->Write();
2943 plots_camet.at(bin).resolHist->Write();
2944 }
2945
2946 fout->Write();
2947
2948 cout << "** Exiting..." << endl;
2949
2950 delete treeReaderElectron;
2951 delete treeReaderMuon;
2952 delete treeReaderPhoton;
2953 delete treeReaderJet;
2954 delete treeReaderBJet;
2955 delete treeReaderTauJet;
2956 delete chainElectron;
2957 delete chainMuon;
2958 delete chainPhoton;
2959 delete chainJet;
2960 delete chainBJet;
2961 delete chainTauJet;
2962}
2963
2964//------------------------------------------------------------------------------
2965
2966int main(int argc, char *argv[])
2967{
2968 char *appName = "DelphesValidation";
2969
2970 if(argc != 12)
2971 {
2972 cout << " Usage: " << appName << " input_file_electron input_file_muon input_file_photon input_file_jet input_file_bjet input_file_taujet output_file version" << endl;
2973 cout << " input_file_pion - input file in ROOT format ('Delphes' tree)," << endl;
2974 cout << " input_file_electron - input file in ROOT format ('Delphes' tree)," << endl;
2975 cout << " input_file_muon - input file in ROOT format ('Delphes' tree)," << endl;
2976 cout << " input_file_photon - input file in ROOT format ('Delphes' tree)," << endl;
2977 cout << " input_file_neutralhadron - input file in ROOT format ('Delphes' tree)," << endl;
2978 cout << " input_file_jet - input file in ROOT format ('Delphes' tree)," << endl;
2979 cout << " input_file_bjet - input file in ROOT format ('Delphes' tree)," << endl;
2980 cout << " input_file_cjet - input file in ROOT format ('Delphes' tree)," << endl;
2981 cout << " input_file_taujet - input file in ROOT format ('Delphes' tree)," << endl;
2982 cout << " output_file - output file in ROOT format," << endl;
2983 cout << " version - Delphes version" << endl;
2984
2985 return 1;
2986 }
2987
2988 gROOT->SetBatch();
2989
2990 int appargc = 1;
2991 char *appargv[] = {appName};
2992 TApplication app(appName, &appargc, appargv);
2993
2994 DelphesValidation(argv[1], argv[2], argv[3], argv[4], argv[5], argv[6], argv[7], argv[8], argv[9], argv[10], argv[11]);
2995}
Note: See TracBrowser for help on using the repository browser.