Fork me on GitHub

source: git/validation/DelphesValidation.cpp@ 78ce8d1

Last change on this file since 78ce8d1 was e09b9bf, checked in by Pavel Demin <pavel.demin@…>, 4 years ago

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