Fork me on GitHub

source: git/examples/Validation.cpp@ 06a83b3

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 06a83b3 was 06a83b3, checked in by Michele Selvaggi <michele.selvaggi@…>, 8 years ago

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