Fork me on GitHub

source: git/validation/Validation.cpp@ d108fdc

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

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