Fork me on GitHub

source: git/validation/DelphesValidation.cpp@ bdf9e2e

Last change on this file since bdf9e2e was 181c061, checked in by Pavel Demin <pavel-demin@…>, 8 years ago

rename Validation.cpp to DelphesValidation.cpp
  • Property mode set to 100644
File size: 107.7 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 DelphesValidation(
1454 const char *inputFilePion,
1455 const char *inputFileElectron,
1456 const char *inputFileMuon,
1457 const char *inputFilePhoton,
1458 const char *inputFileNeutralHadron,
1459 const char *inputFileJet,
1460 const char *inputFileBJet,
1461 const char *inputFileCJet,
1462 const char *inputFileTauJet,
1463 const char *outputFile,
1464 const char *version)
1465{
1466
1467 TChain *chainPion = new TChain("Delphes");
1468 chainPion->Add(inputFilePion);
1469 ExRootTreeReader *treeReaderPion = new ExRootTreeReader(chainPion);
1470
1471 TChain *chainElectron = new TChain("Delphes");
1472 chainElectron->Add(inputFileElectron);
1473 ExRootTreeReader *treeReaderElectron = new ExRootTreeReader(chainElectron);
1474
1475 TChain *chainMuon = new TChain("Delphes");
1476 chainMuon->Add(inputFileMuon);
1477 ExRootTreeReader *treeReaderMuon = new ExRootTreeReader(chainMuon);
1478
1479 TChain *chainPhoton = new TChain("Delphes");
1480 chainPhoton->Add(inputFilePhoton);
1481 ExRootTreeReader *treeReaderPhoton = new ExRootTreeReader(chainPhoton);
1482
1483 TChain *chainNeutralHadron = new TChain("Delphes");
1484 chainNeutralHadron->Add(inputFileNeutralHadron);
1485 ExRootTreeReader *treeReaderNeutralHadron = new ExRootTreeReader(chainNeutralHadron);
1486
1487 TChain *chainJet = new TChain("Delphes");
1488 chainJet->Add(inputFileJet);
1489 ExRootTreeReader *treeReaderJet = new ExRootTreeReader(chainJet);
1490
1491 TChain *chainBJet = new TChain("Delphes");
1492 chainBJet->Add(inputFileBJet);
1493 ExRootTreeReader *treeReaderBJet = new ExRootTreeReader(chainBJet);
1494
1495 TChain *chainCJet = new TChain("Delphes");
1496 chainCJet->Add(inputFileCJet);
1497 ExRootTreeReader *treeReaderCJet = new ExRootTreeReader(chainCJet);
1498
1499 TChain *chainTauJet = new TChain("Delphes");
1500 chainTauJet->Add(inputFileTauJet);
1501 ExRootTreeReader *treeReaderTauJet = new ExRootTreeReader(chainTauJet);
1502
1503 TClonesArray *branchParticleElectron = treeReaderElectron->UseBranch("Particle");
1504 TClonesArray *branchTrackElectron = treeReaderElectron->UseBranch("Track");
1505 TClonesArray *branchElectron = treeReaderElectron->UseBranch("Electron");
1506 TClonesArray *branchElectronPF = treeReaderElectron->UseBranch("ElectronPF");
1507
1508 TClonesArray *branchParticleMuon = treeReaderMuon->UseBranch("Particle");
1509 TClonesArray *branchTrackMuon = treeReaderMuon->UseBranch("Track");
1510 TClonesArray *branchMuon = treeReaderMuon->UseBranch("Muon");
1511
1512 TClonesArray *branchParticlePion = treeReaderPion->UseBranch("Particle");
1513 TClonesArray *branchTrackPion = treeReaderPion->UseBranch("Track");
1514 TClonesArray *branchPion = treeReaderPion->UseBranch("Pion");
1515
1516 TClonesArray *branchParticlePhoton = treeReaderPhoton->UseBranch("Particle");
1517 TClonesArray *branchTowerPhoton = treeReaderPhoton->UseBranch("Tower");
1518 TClonesArray *branchPhoton = treeReaderPhoton->UseBranch("Photon");
1519
1520 TClonesArray *branchParticleNeutralHadron = treeReaderNeutralHadron->UseBranch("Particle");
1521 TClonesArray *branchTowerNeutralHadron = treeReaderNeutralHadron->UseBranch("Tower");
1522
1523 TClonesArray *branchGenJet = treeReaderJet->UseBranch("GenJet");
1524 TClonesArray *branchParticleJet = treeReaderJet->UseBranch("Particle");
1525 TClonesArray *branchPFJet = treeReaderJet->UseBranch("PFJet");
1526 TClonesArray *branchCaloJet = treeReaderJet->UseBranch("CaloJet");
1527 TClonesArray *branchJet = treeReaderJet->UseBranch("Jet");
1528
1529 TClonesArray *branchParticleBJet = treeReaderBJet->UseBranch("Particle");
1530 TClonesArray *branchPFBJet = treeReaderBJet->UseBranch("Jet");
1531
1532 TClonesArray *branchParticleCJet = treeReaderCJet->UseBranch("Particle");
1533 TClonesArray *branchPFCJet = treeReaderCJet->UseBranch("Jet");
1534
1535 TClonesArray *branchParticleTauJet = treeReaderTauJet->UseBranch("Particle");
1536 TClonesArray *branchPFTauJet = treeReaderTauJet->UseBranch("Jet");
1537
1538 TClonesArray *branchGenScalarHT = treeReaderJet->UseBranch("GenScalarHT");
1539 TClonesArray *branchMet = treeReaderJet->UseBranch("PFMissingET");
1540 TClonesArray *branchCaloMet = treeReaderJet->UseBranch("CaloMissingET");
1541
1542 std::vector<Color_t> colors;
1543
1544 colors.push_back(kBlack);
1545 colors.push_back(kBlue);
1546 colors.push_back(kRed);
1547 colors.push_back(kGreen+1);
1548 colors.push_back(kMagenta+1);
1549 colors.push_back(kOrange);
1550
1551 std::vector<Int_t> markerStyles;
1552
1553 markerStyles.push_back(20);
1554 markerStyles.push_back(21);
1555 markerStyles.push_back(22);
1556 markerStyles.push_back(23);
1557 markerStyles.push_back(33);
1558 markerStyles.push_back(34);
1559
1560 TString pdfOutput(outputFile);
1561 pdfOutput.ReplaceAll(".root", ".pdf");
1562
1563 TString figPath = inputFilePion;
1564 figPath.ReplaceAll(".root", "");
1565 figPath.ReplaceAll("root", "www/fig");
1566 Int_t lastSlash = figPath.Last('/');
1567 Int_t sizePath = figPath.Length();
1568 figPath.Remove(lastSlash+1,sizePath);
1569
1570 TString header = pdfOutput;
1571 header.ReplaceAll(".pdf", "");
1572 header.ReplaceAll("validation_", "");
1573 lastSlash = header.Last('/');
1574 sizePath = header.Length();
1575 header.Remove(0,lastSlash+1);
1576
1577 TString vrs(version);
1578
1579 TPaveText *pave = new TPaveText(0.0, 0.89, 0.94, 0.94,"NDC");
1580 pave->SetTextAlign(30);
1581 pave->SetTextFont(132);
1582 pave->SetBorderSize(0);
1583 pave->SetShadowColor(0);
1584 pave->SetFillColor(0);
1585 pave->SetFillStyle(0);
1586 pave->AddText("Delphes "+vrs+" - "+header);
1587
1588 TString s_etaMin, s_etaMax, s_eta, s_pt, s_e;
1589
1590 Double_t ptMin = 1.;
1591 Double_t ptMax = 50000.;
1592
1593 Double_t etaMin = -6.;
1594 Double_t etaMax = 6.;
1595
1596 std::vector<Double_t> ptVals;
1597
1598 ptVals.push_back(1.);
1599 ptVals.push_back(10.);
1600 ptVals.push_back(100.);
1601 ptVals.push_back(1000.);
1602 ptVals.push_back(10000.);
1603
1604 std::vector<Double_t> etaVals;
1605
1606 etaVals.push_back(0.);
1607 etaVals.push_back(1.5);
1608 etaVals.push_back(2.5);
1609 etaVals.push_back(4.0);
1610 etaVals.push_back(6.0);
1611
1612 const int n_etabins = etaVals.size()-1;
1613 const int n_ptbins = ptVals.size();
1614
1615 //////////////////////////
1616 // Tracking performance //
1617 //////////////////////////
1618
1619 // --------- Pion Tracks --------- //
1620
1621 TMultiGraph *mg_trkpi_res_pt = new TMultiGraph("","");
1622 TMultiGraph *mg_trkpi_eff_pt = new TMultiGraph("","");
1623 TMultiGraph *mg_trkpi_res_eta = new TMultiGraph("","");
1624 TMultiGraph *mg_trkpi_eff_eta = new TMultiGraph("","");
1625
1626 TLegend *leg_trkpi_res_pt = new TLegend(0.55,0.22,0.90,0.48);
1627 TLegend *leg_trkpi_eff_pt = (TLegend*)leg_trkpi_res_pt->Clone();
1628 TLegend *leg_trkpi_res_eta = (TLegend*)leg_trkpi_res_pt->Clone();
1629 TLegend *leg_trkpi_eff_eta = (TLegend*)leg_trkpi_res_eta->Clone();
1630
1631 TGraphErrors *gr_trkpi_res_pt = new TGraphErrors[n_etabins];
1632 TGraphErrors *gr_trkpi_eff_pt = new TGraphErrors[n_etabins];
1633 TGraphErrors *gr_trkpi_res_eta = new TGraphErrors[n_ptbins];
1634 TGraphErrors *gr_trkpi_eff_eta = new TGraphErrors[n_ptbins];
1635 TH1D* h_trkpi_eff_pt, *h_trkpi_eff_eta;
1636
1637 std::vector<resolPlot> *plots_trkpi_res_pt = new std::vector<resolPlot>[n_etabins];
1638 std::vector<resolPlot> *plots_trkpi_res_eta = new std::vector<resolPlot>[n_ptbins];
1639
1640 // loop over eta bins
1641 for (k = 0; k < etaVals.size()-1; k++)
1642 {
1643 HistogramsCollection(&plots_trkpi_res_pt[k], TMath::Log10(ptMin), TMath::Log10(ptMax), "trkpi");
1644 GetPtres<Track>(&plots_trkpi_res_pt[k], branchTrackPion, branchParticlePion, 211, etaVals.at(k), etaVals.at(k+1), treeReaderPion);
1645 gr_trkpi_res_pt[k] = EresGraph(&plots_trkpi_res_pt[k]);
1646
1647 h_trkpi_eff_pt = GetEffPt<Track>(branchTrackPion, branchParticlePion, "Pion", 211, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderPion);
1648 gr_trkpi_eff_pt[k] = TGraphErrors(h_trkpi_eff_pt);
1649
1650 s_etaMin = Form("%.1f",etaVals.at(k));
1651 s_etaMax = Form("%.1f",etaVals.at(k+1));
1652
1653 s_eta = "#pi^{ #pm} , " + s_etaMin + " < | #eta | < "+s_etaMax;
1654
1655 gr_trkpi_res_pt[k].SetName("trkRes_"+s_etaMin+"_"+s_etaMax);
1656 gr_trkpi_eff_pt[k].SetName("trkEff_"+s_etaMin+"_"+s_etaMax);
1657
1658 addResoGraph(mg_trkpi_res_pt, &gr_trkpi_res_pt[k], leg_trkpi_res_pt, markerStyles.at(k), colors.at(k), s_eta);
1659 addResoGraph(mg_trkpi_eff_pt, &gr_trkpi_eff_pt[k], leg_trkpi_eff_pt, markerStyles.at(k), colors.at(k), s_eta);
1660 }
1661
1662 // loop over pt
1663 for (k = 0; k < ptVals.size(); k++)
1664 {
1665 HistogramsCollectionVsEta(&plots_trkpi_res_eta[k], etaMin, etaMax, 0.5*ptVals.at(k), 2.0*ptVals.at(k), "trkpi", 0.0, 2.0);
1666 GetPtresVsEta<Track>(&plots_trkpi_res_eta[k], branchTrackPion, branchParticlePion, 211, 0.5*ptVals.at(k), 2.0*ptVals.at(k), treeReaderPion);
1667 gr_trkpi_res_eta[k] = EresGraphVsEta(&plots_trkpi_res_eta[k]);
1668
1669 h_trkpi_eff_eta = GetEffEta<Track>(branchTrackPion, branchParticlePion, "Pion", 211, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderPion);
1670 gr_trkpi_eff_eta[k] = TGraphErrors(h_trkpi_eff_eta);
1671
1672 s_pt = Form("#pi^{ #pm} , p_{T} = %.0f GeV",ptVals.at(k));
1673 if(ptVals.at(k) >= 1000.) s_pt = Form("#pi^{ #pm} , p_{T} = %.0f TeV",ptVals.at(k)/1000.);
1674
1675 addResoGraph(mg_trkpi_res_eta, &gr_trkpi_res_eta[k], leg_trkpi_res_eta, markerStyles.at(k), colors.at(k), s_pt );
1676 addResoGraph(mg_trkpi_eff_eta, &gr_trkpi_eff_eta[k], leg_trkpi_eff_eta, markerStyles.at(k), colors.at(k), s_pt );
1677
1678 }
1679
1680 TCanvas *c_trkpi_res_pt = new TCanvas("","", 800, 600);
1681
1682 mg_trkpi_res_pt->Draw("APE");
1683 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);
1684 leg_trkpi_res_pt->Draw();
1685 pave->Draw();
1686
1687 c_trkpi_res_pt->Print(pdfOutput+"(","pdf");
1688 c_trkpi_res_pt->Print(figPath+"img_trkpi_res_pt.pdf","pdf");
1689 c_trkpi_res_pt->Print(figPath+"img_trkpi_res_pt.png","png");
1690
1691 TCanvas *c_trkpi_res_eta = new TCanvas("","", 800, 600);
1692
1693 mg_trkpi_res_eta->Draw("APE");
1694 DrawAxis(mg_trkpi_res_eta, leg_trkpi_res_eta, etaMin, etaMax, 0.01, 100, " #eta ", "(track resolution in p_{T})/p_{T} (%)", false, true);
1695 leg_trkpi_res_eta->Draw();
1696 pave->Draw();
1697
1698 c_trkpi_res_eta->Print(pdfOutput,"pdf");
1699 c_trkpi_res_eta->Print(figPath+"img_trkpi_res_eta.pdf","pdf");
1700 c_trkpi_res_eta->Print(figPath+"img_trkpi_res_eta.png","png");
1701
1702 TCanvas *c_trkpi_eff_pt = new TCanvas("","", 800, 600);
1703
1704 mg_trkpi_eff_pt->Draw("APE");
1705 DrawAxis(mg_trkpi_eff_pt, leg_trkpi_eff_pt, ptMin, ptMax, 0.0, 100, "p_{T} [GeV]", "tracking efficiency (%)", true, false);
1706 leg_trkpi_eff_pt->Draw();
1707 pave->Draw();
1708
1709 c_trkpi_eff_pt->Print(pdfOutput,"pdf");
1710 c_trkpi_eff_pt->Print(figPath+"img_trkpi_eff_pt.pdf","pdf");
1711 c_trkpi_eff_pt->Print(figPath+"img_trkpi_eff_pt.png","png");
1712
1713 TCanvas *c_trkpi_eff_eta = new TCanvas("","", 800, 600);
1714
1715 mg_trkpi_eff_eta->Draw("APE");
1716 DrawAxis(mg_trkpi_eff_eta, leg_trkpi_eff_eta, etaMin, etaMax, 0.0, 100, " #eta ", "tracking efficiency (%)", false, false);
1717 leg_trkpi_eff_eta->Draw();
1718 pave->Draw();
1719
1720
1721 c_trkpi_eff_eta->Print(pdfOutput,"pdf");
1722 c_trkpi_eff_eta->Print(figPath+"img_trkpi_eff_eta.pdf","pdf");
1723 c_trkpi_eff_eta->Print(figPath+"img_trkpi_eff_eta.png","png");
1724
1725 // --------- Electron Tracks --------- //
1726
1727 TMultiGraph *mg_trkele_res_pt = new TMultiGraph("","");
1728 TMultiGraph *mg_trkele_eff_pt = new TMultiGraph("","");
1729 TMultiGraph *mg_trkele_res_eta = new TMultiGraph("","");
1730 TMultiGraph *mg_trkele_eff_eta = new TMultiGraph("","");
1731
1732 TLegend *leg_trkele_res_pt = new TLegend(0.55,0.22,0.90,0.48);
1733 TLegend *leg_trkele_eff_pt = (TLegend*)leg_trkele_res_pt->Clone();
1734 TLegend *leg_trkele_res_eta = (TLegend*)leg_trkele_res_pt->Clone();
1735 TLegend *leg_trkele_eff_eta = (TLegend*)leg_trkele_res_eta->Clone();
1736
1737 TGraphErrors *gr_trkele_res_pt = new TGraphErrors[n_etabins];
1738 TGraphErrors *gr_trkele_eff_pt = new TGraphErrors[n_etabins];
1739 TGraphErrors *gr_trkele_res_eta = new TGraphErrors[n_ptbins];
1740 TGraphErrors *gr_trkele_eff_eta = new TGraphErrors[n_ptbins];
1741
1742 TH1D* h_trkele_eff_pt, *h_trkele_eff_eta;
1743
1744 std::vector<resolPlot> *plots_trkele_res_pt = new std::vector<resolPlot>[n_etabins];
1745 std::vector<resolPlot> *plots_trkele_res_eta = new std::vector<resolPlot>[n_ptbins];
1746
1747 // loop over eta bins
1748 for (k = 0; k < etaVals.size()-1; k++)
1749 {
1750 HistogramsCollection(&plots_trkele_res_pt[k], TMath::Log10(ptMin), TMath::Log10(ptMax), "trkele");
1751 GetPtres<Track>(&plots_trkele_res_pt[k], branchTrackElectron, branchParticleElectron, 11, etaVals.at(k), etaVals.at(k+1), treeReaderElectron);
1752 gr_trkele_res_pt[k] = EresGraph(&plots_trkele_res_pt[k]);
1753
1754 h_trkele_eff_pt = GetEffPt<Track>(branchTrackElectron, branchParticleElectron, "Electron", 11, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderElectron);
1755 gr_trkele_eff_pt[k] = TGraphErrors(h_trkele_eff_pt);
1756
1757 s_etaMin = Form("%.1f",etaVals.at(k));
1758 s_etaMax = Form("%.1f",etaVals.at(k+1));
1759
1760 s_eta = "e^{ #pm} , " + s_etaMin + " < | #eta | < "+s_etaMax;
1761
1762 gr_trkele_res_pt[k].SetName("trkRes_"+s_etaMin+"_"+s_etaMax);
1763 gr_trkele_eff_pt[k].SetName("trkEff_"+s_etaMin+"_"+s_etaMax);
1764
1765 addResoGraph(mg_trkele_res_pt, &gr_trkele_res_pt[k], leg_trkele_res_pt, markerStyles.at(k), colors.at(k), s_eta);
1766 addResoGraph(mg_trkele_eff_pt, &gr_trkele_eff_pt[k], leg_trkele_eff_pt, markerStyles.at(k), colors.at(k), s_eta);
1767 }
1768
1769 // loop over pt
1770 for (k = 0; k < ptVals.size(); k++)
1771 {
1772 HistogramsCollectionVsEta(&plots_trkele_res_eta[k], etaMin, etaMax, 0.5*ptVals.at(k), 2.0*ptVals.at(k), "trkele", 0.0, 2.0);
1773 GetPtresVsEta<Track>(&plots_trkele_res_eta[k], branchTrackElectron, branchParticleElectron, 11, 0.5*ptVals.at(k), 2.0*ptVals.at(k), treeReaderElectron);
1774 gr_trkele_res_eta[k] = EresGraphVsEta(&plots_trkele_res_eta[k]);
1775
1776 h_trkele_eff_eta = GetEffEta<Track>(branchTrackElectron, branchParticleElectron, "Electron", 11, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderElectron);
1777 gr_trkele_eff_eta[k] = TGraphErrors(h_trkele_eff_eta);
1778
1779 s_pt = Form("e^{ #pm} , p_{T} = %.0f GeV",ptVals.at(k));
1780 if(ptVals.at(k) >= 1000.) s_pt = Form("e^{ #pm} , p_{T} = %.0f TeV",ptVals.at(k)/1000.);
1781
1782 addResoGraph(mg_trkele_res_eta, &gr_trkele_res_eta[k], leg_trkele_res_eta, markerStyles.at(k), colors.at(k), s_pt );
1783 addResoGraph(mg_trkele_eff_eta, &gr_trkele_eff_eta[k], leg_trkele_eff_eta, markerStyles.at(k), colors.at(k), s_pt );
1784
1785 }
1786
1787 TCanvas *c_trkele_res_pt = new TCanvas("","", 800, 600);
1788
1789 mg_trkele_res_pt->Draw("APE");
1790 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);
1791 leg_trkele_res_pt->Draw();
1792 pave->Draw();
1793
1794 c_trkele_res_pt->Print(pdfOutput,"pdf");
1795 c_trkele_res_pt->Print(figPath+"img_trkele_res_pt.pdf","pdf");
1796 c_trkele_res_pt->Print(figPath+"img_trkele_res_pt.png","png");
1797
1798 TCanvas *c_trkele_res_eta = new TCanvas("","", 800, 600);
1799
1800 mg_trkele_res_eta->Draw("APE");
1801 DrawAxis(mg_trkele_res_eta, leg_trkele_res_eta, etaMin, etaMax, 0.01, 100, " #eta ", "(track resolution in p_{T})/p_{T} (%)", false, true);
1802 leg_trkele_res_eta->Draw();
1803 pave->Draw();
1804
1805 c_trkele_res_eta->Print(pdfOutput,"pdf");
1806 c_trkele_res_eta->Print(figPath+"img_trkele_res_eta.pdf","pdf");
1807 c_trkele_res_eta->Print(figPath+"img_trkele_res_eta.png","png");
1808
1809 TCanvas *c_trkele_eff_pt = new TCanvas("","", 800, 600);
1810
1811 mg_trkele_eff_pt->Draw("APE");
1812 DrawAxis(mg_trkele_eff_pt, leg_trkele_eff_pt, ptMin, ptMax, 0.0, 100, "p_{T} [GeV]", "tracking efficiency (%)", true, false);
1813 leg_trkele_eff_pt->Draw();
1814 pave->Draw();
1815
1816 c_trkele_eff_pt->Print(pdfOutput,"pdf");
1817 c_trkele_eff_pt->Print(figPath+"img_trkele_eff_pt.pdf","pdf");
1818 c_trkele_eff_pt->Print(figPath+"img_trkele_eff_pt.png","png");
1819
1820 TCanvas *c_trkele_eff_eta = new TCanvas("","", 800, 600);
1821
1822 mg_trkele_eff_eta->Draw("APE");
1823 DrawAxis(mg_trkele_eff_eta, leg_trkele_eff_eta, etaMin, etaMax, 0.0, 100, " #eta ", "tracking efficiency (%)", false, false);
1824 leg_trkele_eff_eta->Draw();
1825 pave->Draw();
1826
1827 c_trkele_eff_eta->Print(pdfOutput,"pdf");
1828 c_trkele_eff_eta->Print(figPath+"img_trkele_eff_eta.pdf","pdf");
1829 c_trkele_eff_eta->Print(figPath+"img_trkele_eff_eta.png","png");
1830
1831
1832 // --------- Muon Tracks --------- //
1833
1834 TMultiGraph *mg_trkmu_res_pt = new TMultiGraph("","");
1835 TMultiGraph *mg_trkmu_eff_pt = new TMultiGraph("","");
1836 TMultiGraph *mg_trkmu_res_eta = new TMultiGraph("","");
1837 TMultiGraph *mg_trkmu_eff_eta = new TMultiGraph("","");
1838
1839 TLegend *leg_trkmu_res_pt = new TLegend(0.55,0.22,0.90,0.48);
1840 TLegend *leg_trkmu_eff_pt = (TLegend*)leg_trkmu_res_pt->Clone();
1841
1842 TLegend *leg_trkmu_res_eta = (TLegend*)leg_trkmu_res_pt->Clone();
1843 TLegend *leg_trkmu_eff_eta = (TLegend*)leg_trkmu_res_eta->Clone();
1844
1845
1846 TGraphErrors *gr_trkmu_res_pt = new TGraphErrors[n_etabins];
1847 TGraphErrors *gr_trkmu_eff_pt = new TGraphErrors[n_etabins];
1848 TGraphErrors *gr_trkmu_res_eta = new TGraphErrors[n_ptbins];
1849 TGraphErrors *gr_trkmu_eff_eta = new TGraphErrors[n_ptbins];
1850
1851 TH1D* h_trkmu_eff_pt, *h_trkmu_eff_eta;
1852
1853 std::vector<resolPlot> *plots_trkmu_res_pt = new std::vector<resolPlot>[n_etabins];
1854 std::vector<resolPlot> *plots_trkmu_res_eta = new std::vector<resolPlot>[n_ptbins];
1855
1856 // loop over eta bins
1857 for (k = 0; k < etaVals.size()-1; k++)
1858 {
1859 HistogramsCollection(&plots_trkmu_res_pt[k], TMath::Log10(ptMin), TMath::Log10(ptMax), "trkmu");
1860 GetPtres<Track>(&plots_trkmu_res_pt[k], branchTrackMuon, branchParticleMuon, 13, etaVals.at(k), etaVals.at(k+1), treeReaderMuon);
1861 gr_trkmu_res_pt[k] = EresGraph(&plots_trkmu_res_pt[k]);
1862
1863 h_trkmu_eff_pt = GetEffPt<Track>(branchTrackMuon, branchParticleMuon, "Muon", 13, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderMuon);
1864 gr_trkmu_eff_pt[k] = TGraphErrors(h_trkmu_eff_pt);
1865
1866 s_etaMin = Form("%.1f",etaVals.at(k));
1867 s_etaMax = Form("%.1f",etaVals.at(k+1));
1868
1869 s_eta = "#mu^{ #pm} , " + s_etaMin + " < | #eta | < "+s_etaMax;
1870
1871 gr_trkmu_res_pt[k].SetName("trkRes_"+s_etaMin+"_"+s_etaMax);
1872 gr_trkmu_eff_pt[k].SetName("trkEff_"+s_etaMin+"_"+s_etaMax);
1873
1874 addResoGraph(mg_trkmu_res_pt, &gr_trkmu_res_pt[k], leg_trkmu_res_pt, markerStyles.at(k), colors.at(k), s_eta);
1875 addResoGraph(mg_trkmu_eff_pt, &gr_trkmu_eff_pt[k], leg_trkmu_eff_pt, markerStyles.at(k), colors.at(k), s_eta);
1876 }
1877
1878 // loop over pt
1879 for (k = 0; k < ptVals.size(); k++)
1880 {
1881 HistogramsCollectionVsEta(&plots_trkmu_res_eta[k], etaMin, etaMax, 0.5*ptVals.at(k), 2.0*ptVals.at(k), "trkmu", 0.0, 2.0);
1882 GetPtresVsEta<Track>(&plots_trkmu_res_eta[k], branchTrackMuon, branchParticleMuon, 13, 0.5*ptVals.at(k), 2.0*ptVals.at(k), treeReaderMuon);
1883 gr_trkmu_res_eta[k] = EresGraphVsEta(&plots_trkmu_res_eta[k]);
1884
1885 h_trkmu_eff_eta = GetEffEta<Track>(branchTrackMuon, branchParticleMuon, "Muon", 13, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderMuon);
1886 gr_trkmu_eff_eta[k] = TGraphErrors(h_trkmu_eff_eta);
1887
1888 s_pt = Form("#mu^{ #pm} , p_{T} = %.0f GeV",ptVals.at(k));
1889 if(ptVals.at(k) >= 1000.) s_pt = Form("#mu^{ #pm} , p_{T} = %.0f TeV",ptVals.at(k)/1000.);
1890
1891 addResoGraph(mg_trkmu_res_eta, &gr_trkmu_res_eta[k], leg_trkmu_res_eta, markerStyles.at(k), colors.at(k), s_pt );
1892 addResoGraph(mg_trkmu_eff_eta, &gr_trkmu_eff_eta[k], leg_trkmu_eff_eta, markerStyles.at(k), colors.at(k), s_pt );
1893
1894 }
1895
1896 TCanvas *c_trkmu_res_pt = new TCanvas("","", 800, 600);
1897
1898 mg_trkmu_res_pt->Draw("APE");
1899 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);
1900 leg_trkmu_res_pt->Draw();
1901 pave->Draw();
1902
1903 c_trkmu_res_pt->Print(pdfOutput,"pdf");
1904 c_trkmu_res_pt->Print(figPath+"img_trkmu_res_pt.pdf","pdf");
1905 c_trkmu_res_pt->Print(figPath+"img_trkmu_res_pt.png","png");
1906
1907 TCanvas *c_trkmu_res_eta = new TCanvas("","", 800, 600);
1908
1909 mg_trkmu_res_eta->Draw("APE");
1910 DrawAxis(mg_trkmu_res_eta, leg_trkmu_res_eta, etaMin, etaMax, 0.01, 100, " #eta ", "(track resolution in p_{T})/p_{T} (%)", false, true);
1911 leg_trkmu_res_eta->Draw();
1912 pave->Draw();
1913
1914 c_trkmu_res_eta->Print(pdfOutput,"pdf");
1915 c_trkmu_res_eta->Print(figPath+"img_trkmu_res_eta.pdf","pdf");
1916 c_trkmu_res_eta->Print(figPath+"img_trkmu_res_eta.png","png");
1917
1918 TCanvas *c_trkmu_eff_pt = new TCanvas("","", 800, 600);
1919
1920 mg_trkmu_eff_pt->Draw("APE");
1921 DrawAxis(mg_trkmu_eff_pt, leg_trkmu_eff_pt, ptMin, ptMax, 0.0, 100, "p_{T} [GeV]", "tracking efficiency (%)", true, false);
1922 leg_trkmu_eff_pt->Draw();
1923 pave->Draw();
1924
1925 c_trkmu_eff_pt->Print(pdfOutput,"pdf");
1926 c_trkmu_eff_pt->Print(figPath+"img_trkmu_eff_pt.pdf","pdf");
1927 c_trkmu_eff_pt->Print(figPath+"img_trkmu_eff_pt.png","png");
1928
1929 TCanvas *c_trkmu_eff_eta = new TCanvas("","", 800, 600);
1930
1931 mg_trkmu_eff_eta->Draw("APE");
1932 DrawAxis(mg_trkmu_eff_eta, leg_trkmu_eff_eta, etaMin, etaMax, 0.0, 100, " #eta ", "tracking efficiency (%)", false, false);
1933 leg_trkmu_eff_eta->Draw();
1934 pave->Draw();
1935
1936 c_trkmu_eff_eta->Print(pdfOutput,"pdf");
1937 c_trkmu_eff_eta->Print(figPath+"img_trkmu_eff_eta.pdf","pdf");
1938 c_trkmu_eff_eta->Print(figPath+"img_trkmu_eff_eta.png","png");
1939
1940
1941 //////////////////////
1942 // Ecal performance //
1943 //////////////////////
1944
1945
1946 TMultiGraph *mg_ecal_res_e = new TMultiGraph("","");
1947 TMultiGraph *mg_ecal_res_eta = new TMultiGraph("","");
1948
1949 TLegend *leg_ecal_res_e = new TLegend(0.55,0.64,0.90,0.90);
1950 TLegend *leg_ecal_res_eta = new TLegend(0.60,0.59,0.95,0.90);
1951
1952 TGraphErrors *gr_ecal_res_e = new TGraphErrors[n_etabins];
1953 TGraphErrors *gr_ecal_res_eta = new TGraphErrors[n_ptbins];
1954
1955 std::vector<resolPlot> *plots_ecal_res_e = new std::vector<resolPlot>[n_etabins];
1956 std::vector<resolPlot> *plots_ecal_res_eta = new std::vector<resolPlot>[n_ptbins];
1957
1958 // loop over eta bins
1959 for (k = 0; k < etaVals.size()-1; k++)
1960 {
1961 HistogramsCollection(&plots_ecal_res_e[k], TMath::Log10(ptMin), TMath::Log10(ptMax), "ecal");
1962 GetEres<Tower>(&plots_ecal_res_e[k], branchTowerPhoton, branchParticlePhoton, 22, etaVals.at(k), etaVals.at(k+1), treeReaderPhoton);
1963 gr_ecal_res_e[k] = EresGraph(&plots_ecal_res_e[k]);
1964
1965 s_etaMin = Form("%.1f",etaVals.at(k));
1966 s_etaMax = Form("%.1f",etaVals.at(k+1));
1967
1968 s_eta = "#gamma , " + s_etaMin + " < | #eta | < "+s_etaMax;
1969
1970 gr_ecal_res_e[k].SetName("trkRes_"+s_etaMin+"_"+s_etaMax);
1971
1972 addResoGraph(mg_ecal_res_e, &gr_ecal_res_e[k], leg_ecal_res_e, markerStyles.at(k), colors.at(k), s_eta);
1973 }
1974
1975 // loop over pt
1976 for (k = 0; k < ptVals.size(); k++)
1977 {
1978 HistogramsCollectionVsEta(&plots_ecal_res_eta[k], etaMin, etaMax, 0.5*ptVals.at(k), 2.0*ptVals.at(k), "ecal", 0.0, 2.0);
1979 GetEresVsEta<Tower>(&plots_ecal_res_eta[k], branchTowerPhoton, branchParticlePhoton, 22, 0.5*ptVals.at(k), 2.0*ptVals.at(k), treeReaderPhoton);
1980 gr_ecal_res_eta[k] = EresGraphVsEta(&plots_ecal_res_eta[k]);
1981
1982 s_e = Form("#gamma , E = %.0f GeV",ptVals.at(k));
1983 if(ptVals.at(k) >= 1000.) s_e = Form("#gamma , E = %.0f TeV",ptVals.at(k)/1000.);
1984
1985 addResoGraph(mg_ecal_res_eta, &gr_ecal_res_eta[k], leg_ecal_res_eta, markerStyles.at(k), colors.at(k), s_e );
1986 }
1987
1988 TCanvas *c_ecal_res_e = new TCanvas("","", 800, 600);
1989
1990 mg_ecal_res_e->Draw("APE");
1991 // DrawAxis(mg_ecal_res_e, leg_ecal_res_e, ptMin, ptMax, 0.5, 100, "E [GeV]", "(ECAL resolution in E)/E (%)", true, true);
1992 DrawAxis(mg_ecal_res_e, leg_ecal_res_e, ptMin, ptMax, 0.0, 20, "E [GeV]", "(ECAL resolution in E)/E (%)", true, false);
1993 leg_ecal_res_e->Draw();
1994 pave->Draw();
1995
1996 c_ecal_res_e->Print(pdfOutput,"pdf");
1997 c_ecal_res_e->Print(figPath+"img_ecal_res_e.pdf","pdf");
1998 c_ecal_res_e->Print(figPath+"img_ecal_res_e.png","png");
1999
2000 TCanvas *c_ecal_res_eta = new TCanvas("","", 800, 600);
2001
2002 mg_ecal_res_eta->Draw("APE");
2003 //DrawAxis(mg_ecal_res_eta, leg_ecal_res_eta, etaMin, etaMax, 0.5, 100, " #eta ", "(ECAL resolution in E)/E (%)", false, true);
2004 DrawAxis(mg_ecal_res_eta, leg_ecal_res_eta, etaMin, etaMax, 0.0, 20, " #eta ", "(ECAL resolution in E)/E (%)", false, false);
2005 leg_ecal_res_eta->Draw();
2006 pave->Draw();
2007
2008 c_ecal_res_eta->Print(pdfOutput,"pdf");
2009 c_ecal_res_eta->Print(figPath+"img_ecal_res_eta.pdf","pdf");
2010 c_ecal_res_eta->Print(figPath+"img_ecal_res_eta.png","png");
2011
2012 //////////////////////
2013 // Hcal performance //
2014 //////////////////////
2015
2016
2017 TMultiGraph *mg_hcal_res_e = new TMultiGraph("","");
2018 TMultiGraph *mg_hcal_res_eta = new TMultiGraph("","");
2019
2020 TLegend *leg_hcal_res_e = new TLegend(0.55,0.64,0.90,0.90);
2021 TLegend *leg_hcal_res_eta = new TLegend(0.60,0.59,0.95,0.90);
2022
2023 TGraphErrors *gr_hcal_res_e = new TGraphErrors[n_etabins];
2024 TGraphErrors *gr_hcal_res_eta = new TGraphErrors[n_ptbins];
2025
2026 std::vector<resolPlot> *plots_hcal_res_e = new std::vector<resolPlot>[n_etabins];
2027 std::vector<resolPlot> *plots_hcal_res_eta = new std::vector<resolPlot>[n_ptbins];
2028
2029 // loop over eta bins
2030 for (k = 0; k < etaVals.size()-1; k++)
2031 {
2032 HistogramsCollection(&plots_hcal_res_e[k], TMath::Log10(ptMin), TMath::Log10(ptMax), "hcal");
2033 GetEres<Tower>(&plots_hcal_res_e[k], branchTowerNeutralHadron, branchParticleNeutralHadron, 2112, etaVals.at(k), etaVals.at(k+1), treeReaderNeutralHadron);
2034
2035 gr_hcal_res_e[k] = EresGraph(&plots_hcal_res_e[k]);
2036
2037 s_etaMin = Form("%.1f",etaVals.at(k));
2038 s_etaMax = Form("%.1f",etaVals.at(k+1));
2039
2040 s_eta = "n , " + s_etaMin + " < | #eta | < "+s_etaMax;
2041
2042 gr_hcal_res_e[k].SetName("trkRes_"+s_etaMin+"_"+s_etaMax);
2043
2044 addResoGraph(mg_hcal_res_e, &gr_hcal_res_e[k], leg_hcal_res_e, markerStyles.at(k), colors.at(k), s_eta);
2045 }
2046
2047 // loop over pt
2048 for (k = 0; k < ptVals.size(); k++)
2049 {
2050 HistogramsCollectionVsEta(&plots_hcal_res_eta[k], etaMin, etaMax, 0.5*ptVals.at(k), 2.0*ptVals.at(k), "hcal", 0.0, 2.0);
2051 GetEresVsEta<Tower>(&plots_hcal_res_eta[k], branchTowerNeutralHadron, branchParticleNeutralHadron, 2112, 0.5*ptVals.at(k), 2.0*ptVals.at(k), treeReaderNeutralHadron);
2052 gr_hcal_res_eta[k] = EresGraphVsEta(&plots_hcal_res_eta[k]);
2053
2054 s_e = Form("n , E = %.0f GeV",ptVals.at(k));
2055 if(ptVals.at(k) >= 1000.) s_e = Form("n , E = %.0f TeV",ptVals.at(k)/1000.);
2056
2057 addResoGraph(mg_hcal_res_eta, &gr_hcal_res_eta[k], leg_hcal_res_eta, markerStyles.at(k), colors.at(k), s_e );
2058 }
2059
2060
2061 TCanvas *c_hcal_res_e = new TCanvas("","", 800, 600);
2062
2063 mg_hcal_res_e->Draw("APE");
2064 //DrawAxis(mg_hcal_res_e, leg_hcal_res_e, ptMin, ptMax, 1, 100, "E [GeV]", "(HCAL resolution in E)/E (%)", true, true);
2065 DrawAxis(mg_hcal_res_e, leg_hcal_res_e, ptMin, ptMax, 0.0, 50, "E [GeV]", "(HCAL resolution in E)/E (%)", true, false);
2066 leg_hcal_res_e->Draw();
2067 pave->Draw();
2068
2069 c_hcal_res_e->Print(pdfOutput,"pdf");
2070 c_hcal_res_e->Print(figPath+"img_hcal_res_e.pdf","pdf");
2071 c_hcal_res_e->Print(figPath+"img_hcal_res_e.png","png");
2072
2073 TCanvas *c_hcal_res_eta = new TCanvas("","", 800, 600);
2074
2075 mg_hcal_res_eta->Draw("APE");
2076 //DrawAxis(mg_hcal_res_eta, leg_hcal_res_eta, etaMin, etaMax, 1, 100, " #eta ", "(HCAL resolution in E)/E (%)", false, true);
2077 DrawAxis(mg_hcal_res_eta, leg_hcal_res_eta, etaMin, etaMax, 0.0, 50, " #eta ", "(HCAL resolution in E)/E (%)", false, false);
2078 leg_hcal_res_eta->Draw();
2079 pave->Draw();
2080
2081 c_hcal_res_eta->Print(pdfOutput,"pdf");
2082 c_hcal_res_eta->Print(figPath+"img_hcal_res_eta.pdf","pdf");
2083 c_hcal_res_eta->Print(figPath+"img_hcal_res_eta.png","png");
2084
2085 ////////////////////
2086 // PF - electrons //
2087 ////////////////////
2088
2089 TMultiGraph *mg_pfele_res_e[n_etabins];
2090 TMultiGraph *mg_pfele_res_eta[n_ptbins];
2091
2092 TLegend *leg_pfele_res_e[n_etabins];
2093 TLegend *leg_pfele_res_eta[n_ptbins];
2094
2095 TGraphErrors *gr_pfele_res_e = new TGraphErrors[n_etabins];
2096 TGraphErrors *gr_pfele_res_eta = new TGraphErrors[n_ptbins];
2097 TGraphErrors *gr_trkele_res_e = new TGraphErrors[n_etabins];
2098 TGraphErrors *gr_trkele_res_eeta = new TGraphErrors[n_ptbins];
2099
2100 std::vector<resolPlot> *plots_pfele_res_e = new std::vector<resolPlot>[n_etabins];
2101 std::vector<resolPlot> *plots_pfele_res_eta = new std::vector<resolPlot>[n_ptbins];
2102 std::vector<resolPlot> *plots_trkele_res_e = new std::vector<resolPlot>[n_etabins];
2103 std::vector<resolPlot> *plots_trkele_res_eeta = new std::vector<resolPlot>[n_ptbins];
2104
2105 TCanvas *c_pfele_res_e[n_etabins];
2106 TCanvas *c_pfele_res_eta[n_ptbins];
2107
2108
2109 // loop over eta bins
2110 for (k = 0; k < etaVals.size()-1; k++)
2111 {
2112 mg_pfele_res_e[k] = new TMultiGraph("","");
2113 leg_pfele_res_e[k] = new TLegend(0.40,0.60,0.75,0.90);
2114
2115 HistogramsCollection(&plots_pfele_res_e[k], TMath::Log10(ptMin), TMath::Log10(ptMax), "pfele");
2116 GetEres<Electron>(&plots_pfele_res_e[k], branchElectronPF, branchParticleElectron, 11, etaVals.at(k), etaVals.at(k+1), treeReaderElectron);
2117 gr_pfele_res_e[k] = EresGraph(&plots_pfele_res_e[k]);
2118
2119 HistogramsCollection(&plots_trkele_res_e[k], TMath::Log10(ptMin), TMath::Log10(ptMax), "trkele");
2120 GetEres<Track>(&plots_trkele_res_e[k], branchTrackElectron, branchParticleElectron, 11, etaVals.at(k), etaVals.at(k+1), treeReaderElectron);
2121 gr_trkele_res_e[k] = EresGraph(&plots_trkele_res_e[k]);
2122
2123 s_etaMin = Form("%.1f",etaVals.at(k));
2124 s_etaMax = Form("%.1f",etaVals.at(k+1));
2125 s_eta = "e^{ #pm}, "+ s_etaMin + " < | #eta | < " + s_etaMax;
2126
2127 leg_pfele_res_e[k]->SetTextFont(132);
2128 leg_pfele_res_e[k]->SetHeader(s_eta);
2129
2130 addResoGraph(mg_pfele_res_e[k], &gr_ecal_res_e[k], leg_pfele_res_e[k], markerStyles.at(0), colors.at(0), "ECAL");
2131 addResoGraph(mg_pfele_res_e[k], &gr_trkele_res_e[k], leg_pfele_res_e[k], markerStyles.at(1), colors.at(1), "Track");
2132 addResoGraph(mg_pfele_res_e[k], &gr_pfele_res_e[k], leg_pfele_res_e[k], markerStyles.at(2), colors.at(2), "Particle-flow");
2133
2134 c_pfele_res_e[k] = new TCanvas("","", 800, 600);
2135
2136 mg_pfele_res_e[k]->Draw("APE");
2137 //DrawAxis(mg_pfele_res_e[k], leg_pfele_res_e[k], ptMin, ptMax, 0.1, 100, "E [GeV]", "(resolution in E)/E (%)", true, true);
2138 DrawAxis(mg_pfele_res_e[k], leg_pfele_res_e[k], ptMin, ptMax, 0.0, 20, "E [GeV]", "(resolution in E)/E (%)", true, false);
2139 leg_pfele_res_e[k]->Draw();
2140 pave->Draw();
2141
2142 TString s_etarange = "eta_"+s_etaMin+"_"+s_etaMax+"_";
2143
2144 c_pfele_res_e[k]->Print(pdfOutput,"pdf");
2145 c_pfele_res_e[k]->Print(figPath+"img_pfele_res_"+s_etarange+"e.pdf","pdf");
2146 c_pfele_res_e[k]->Print(figPath+"img_pfele_res_"+s_etarange+"e.png","png");
2147
2148 }
2149
2150
2151 // loop over eta bins
2152 for (k = 0; k < ptVals.size(); k++)
2153 {
2154
2155 mg_pfele_res_eta[k] = new TMultiGraph("","");
2156 leg_pfele_res_eta[k] = new TLegend(0.40,0.60,0.75,0.90);
2157
2158 HistogramsCollectionVsEta(&plots_pfele_res_eta[k], etaMin, etaMax, 0.5*ptVals.at(k), 2.0*ptVals.at(k), "pfele", 0.0, 2.0);
2159 GetEresVsEta<Electron>(&plots_pfele_res_eta[k], branchElectronPF, branchParticleElectron, 11, 0.5*ptVals.at(k), 2.0*ptVals.at(k), treeReaderElectron);
2160 gr_pfele_res_eta[k] = EresGraphVsEta(&plots_pfele_res_eta[k]);
2161
2162 HistogramsCollectionVsEta(&plots_trkele_res_eeta[k], etaMin, etaMax, 0.5*ptVals.at(k), 2.0*ptVals.at(k), "trkele", 0.0, 2.0);
2163 GetEresVsEta<Track>(&plots_trkele_res_eeta[k], branchTrackElectron, branchParticleElectron, 11, 0.5*ptVals.at(k), 2.0*ptVals.at(k), treeReaderElectron);
2164 gr_trkele_res_eeta[k] = EresGraphVsEta(&plots_trkele_res_eeta[k]);
2165
2166 s_e = Form("e^{ #pm}, E = %.0f GeV",ptVals.at(k));
2167 if(ptVals.at(k) >= 1000.) s_e = Form("e^{ #pm}, E = %.0f TeV",ptVals.at(k)/1000.);
2168
2169
2170 leg_pfele_res_eta[k]->SetTextFont(132);
2171 leg_pfele_res_eta[k]->SetHeader(s_e);
2172
2173 addResoGraph(mg_pfele_res_eta[k], &gr_ecal_res_eta[k], leg_pfele_res_eta[k], markerStyles.at(0), colors.at(0), "ECAL");
2174 addResoGraph(mg_pfele_res_eta[k], &gr_trkele_res_eeta[k], leg_pfele_res_eta[k], markerStyles.at(1), colors.at(1), "Track");
2175 addResoGraph(mg_pfele_res_eta[k], &gr_pfele_res_eta[k], leg_pfele_res_eta[k], markerStyles.at(2), colors.at(2), "Particle-flow");
2176
2177 c_pfele_res_eta[k] = new TCanvas("","", 800, 600);
2178
2179 mg_pfele_res_eta[k]->Draw("APE");
2180 //DrawAxis(mg_pfele_res_eta[k], leg_pfele_res_eta[k], etaMin, etaMax, 0.1, 1000, "#eta", "(resolution in E)/E (%)", false, true);
2181 DrawAxis(mg_pfele_res_eta[k], leg_pfele_res_eta[k], etaMin, etaMax, 0.0, 50, "#eta", "(resolution in E)/E (%)", false, false);
2182 leg_pfele_res_eta[k]->Draw();
2183 pave->Draw();
2184
2185 TString s_ptrange = Form("pt_%.0f_",ptVals.at(k));
2186
2187 c_pfele_res_eta[k]->Print(pdfOutput,"pdf");
2188 c_pfele_res_eta[k]->Print(figPath+"img_pfele_res_"+s_ptrange+"eta.pdf","pdf");
2189 c_pfele_res_eta[k]->Print(figPath+"img_pfele_res_"+s_ptrange+"eta.png","png");
2190
2191 }
2192
2193 /////////////////
2194 // PF - Pions //
2195 /////////////////
2196
2197 TMultiGraph *mg_pfpi_res_e[n_etabins];
2198 TMultiGraph *mg_pfpi_res_eta[n_ptbins];
2199
2200 TLegend *leg_pfpi_res_e[n_etabins];
2201 TLegend *leg_pfpi_res_eta[n_ptbins];
2202
2203 TGraphErrors *gr_pfpi_res_e = new TGraphErrors[n_etabins];
2204 TGraphErrors *gr_pfpi_res_eta = new TGraphErrors[n_ptbins];
2205
2206 TGraphErrors *gr_trkpi_res_e = new TGraphErrors[n_etabins];
2207 TGraphErrors *gr_trkpi_res_eeta = new TGraphErrors[n_ptbins];
2208
2209 std::vector<resolPlot> *plots_pfpi_res_e = new std::vector<resolPlot>[n_etabins];
2210 std::vector<resolPlot> *plots_pfpi_res_eta = new std::vector<resolPlot>[n_ptbins];
2211 std::vector<resolPlot> *plots_trkpi_res_e = new std::vector<resolPlot>[n_etabins];
2212 std::vector<resolPlot> *plots_trkpi_res_eeta = new std::vector<resolPlot>[n_ptbins];
2213
2214 TCanvas *c_pfpi_res_e[n_etabins];
2215 TCanvas *c_pfpi_res_eta[n_ptbins];
2216
2217
2218 // loop over eta bins
2219 for (k = 0; k < etaVals.size()-1; k++)
2220 {
2221 mg_pfpi_res_e[k] = new TMultiGraph("","");
2222 leg_pfpi_res_e[k] = new TLegend(0.40,0.60,0.75,0.90);
2223
2224 HistogramsCollection(&plots_pfpi_res_e[k], TMath::Log10(ptMin), TMath::Log10(ptMax), "pfpi");
2225 GetEres<Track>(&plots_pfpi_res_e[k], branchPion, branchParticlePion, 211, etaVals.at(k), etaVals.at(k+1), treeReaderPion);
2226 gr_pfpi_res_e[k] = EresGraph(&plots_pfpi_res_e[k]);
2227
2228 HistogramsCollection(&plots_trkpi_res_e[k], TMath::Log10(ptMin), TMath::Log10(ptMax), "trkpi");
2229 GetEres<Track>(&plots_trkpi_res_e[k], branchTrackPion, branchParticlePion, 211, etaVals.at(k), etaVals.at(k+1), treeReaderPion);
2230 gr_trkpi_res_e[k] = EresGraph(&plots_trkpi_res_e[k]);
2231
2232
2233 s_etaMin = Form("%.1f",etaVals.at(k));
2234 s_etaMax = Form("%.1f",etaVals.at(k+1));
2235 s_eta = "#pi^{ #pm}, "+ s_etaMin + " < | #eta | < " + s_etaMax;
2236
2237 leg_pfpi_res_e[k]->SetTextFont(132);
2238 leg_pfpi_res_e[k]->SetHeader(s_eta);
2239
2240 addResoGraph(mg_pfpi_res_e[k], &gr_hcal_res_e[k], leg_pfpi_res_e[k], markerStyles.at(0), colors.at(0), "HCAL");
2241 addResoGraph(mg_pfpi_res_e[k], &gr_trkpi_res_e[k], leg_pfpi_res_e[k], markerStyles.at(1), colors.at(1), "Track");
2242 addResoGraph(mg_pfpi_res_e[k], &gr_pfpi_res_e[k], leg_pfpi_res_e[k], markerStyles.at(2), colors.at(2), "Particle-flow");
2243
2244 c_pfpi_res_e[k] = new TCanvas("","", 800, 600);
2245
2246 mg_pfpi_res_e[k]->Draw("APE");
2247 //DrawAxis(mg_pfpi_res_e[k], leg_pfpi_res_e[k], ptMin, ptMax, 0.1, 100, "E [GeV]", "(resolution in E)/E (%)", true, true);
2248 DrawAxis(mg_pfpi_res_e[k], leg_pfpi_res_e[k], ptMin, ptMax, 0.1, 50, "E [GeV]", "(resolution in E)/E (%)", true, false);
2249 leg_pfpi_res_e[k]->Draw();
2250 pave->Draw();
2251
2252 TString s_etarange = "eta_"+s_etaMin+"_"+s_etaMax+"_";
2253
2254 c_pfpi_res_e[k]->Print(pdfOutput,"pdf");
2255 c_pfpi_res_e[k]->Print(figPath+"img_pfpi_res_"+s_etarange+"e.pdf","pdf");
2256 c_pfpi_res_e[k]->Print(figPath+"img_pfpi_res_"+s_etarange+"e.png","png");
2257
2258 }
2259
2260
2261 // loop over eta bins
2262 for (k = 0; k < ptVals.size(); k++)
2263 {
2264
2265 mg_pfpi_res_eta[k] = new TMultiGraph("","");
2266 leg_pfpi_res_eta[k] = new TLegend(0.40,0.60,0.75,0.90);
2267
2268 HistogramsCollectionVsEta(&plots_pfpi_res_eta[k], etaMin, etaMax, 0.5*ptVals.at(k), 2.0*ptVals.at(k), "pfpi", 0.0, 2.0);
2269 GetEresVsEta<Track>(&plots_pfpi_res_eta[k], branchPion, branchParticlePion, 211, 0.5*ptVals.at(k), 2.0*ptVals.at(k), treeReaderPion);
2270 gr_pfpi_res_eta[k] = EresGraphVsEta(&plots_pfpi_res_eta[k]);
2271
2272 HistogramsCollectionVsEta(&plots_trkpi_res_eeta[k], etaMin, etaMax, 0.5*ptVals.at(k), 2.0*ptVals.at(k), "trkpi", 0.0, 2.0);
2273 GetEresVsEta<Track>(&plots_trkpi_res_eeta[k], branchPion, branchParticlePion, 211, 0.5*ptVals.at(k), 2.0*ptVals.at(k), treeReaderPion);
2274 gr_trkpi_res_eeta[k] = EresGraphVsEta(&plots_trkpi_res_eeta[k]);
2275
2276
2277 s_e = Form("#pi^{ #pm}, E = %.0f GeV",ptVals.at(k));
2278 if(ptVals.at(k) >= 1000.) s_e = Form("#pi^{ #pm}, E = %.0f TeV",ptVals.at(k)/1000.);
2279
2280 leg_pfpi_res_eta[k]->SetTextFont(132);
2281 leg_pfpi_res_eta[k]->SetHeader(s_e);
2282
2283 addResoGraph(mg_pfpi_res_eta[k], &gr_hcal_res_eta[k], leg_pfpi_res_eta[k], markerStyles.at(0), colors.at(0), "HCAL");
2284 addResoGraph(mg_pfpi_res_eta[k], &gr_trkpi_res_eeta[k], leg_pfpi_res_eta[k], markerStyles.at(1), colors.at(1), "Track");
2285 addResoGraph(mg_pfpi_res_eta[k], &gr_pfpi_res_eta[k], leg_pfpi_res_eta[k], markerStyles.at(2), colors.at(2), "Particle-flow");
2286
2287 c_pfpi_res_eta[k] = new TCanvas("","", 800, 600);
2288
2289 mg_pfpi_res_eta[k]->Draw("APE");
2290 //DrawAxis(mg_pfpi_res_eta[k], leg_pfpi_res_eta[k], etaMin, etaMax, 0.1, 1000, "#eta", "(resolution in E)/E (%)", false, true);
2291 DrawAxis(mg_pfpi_res_eta[k], leg_pfpi_res_eta[k], etaMin, etaMax, 0.0, 50, "#eta", "(resolution in E)/E (%)", false, false);
2292 leg_pfpi_res_eta[k]->Draw();
2293 pave->Draw();
2294
2295 TString s_ptrange = Form("pt_%.0f_",ptVals.at(k));
2296
2297 c_pfpi_res_eta[k]->Print(pdfOutput,"pdf");
2298 c_pfpi_res_eta[k]->Print(figPath+"img_pfpi_res_"+s_ptrange+"eta.pdf","pdf");
2299 c_pfpi_res_eta[k]->Print(figPath+"img_pfpi_res_"+s_ptrange+"eta.png","png");
2300
2301 }
2302
2303
2304 /////////////////
2305 // PF - jets //
2306 /////////////////
2307
2308 TMultiGraph *mg_pfjet_res_e[n_etabins];
2309 TMultiGraph *mg_pfjet_res_eta[n_ptbins];
2310
2311 TLegend *leg_pfjet_res_e[n_etabins];
2312 TLegend *leg_pfjet_res_eta[n_ptbins];
2313
2314 TGraphErrors *gr_pfjet_res_e = new TGraphErrors[n_etabins];
2315 TGraphErrors *gr_pfjet_res_eta = new TGraphErrors[n_ptbins];
2316
2317 TGraphErrors *gr_cajet_res_e = new TGraphErrors[n_etabins];
2318 TGraphErrors *gr_cajet_res_eta = new TGraphErrors[n_ptbins];
2319
2320 std::vector<resolPlot> *plots_pfjet_res_e = new std::vector<resolPlot>[n_etabins];
2321 std::vector<resolPlot> *plots_pfjet_res_eta = new std::vector<resolPlot>[n_ptbins];
2322 std::vector<resolPlot> *plots_cajet_res_e = new std::vector<resolPlot>[n_etabins];
2323 std::vector<resolPlot> *plots_cajet_res_eta = new std::vector<resolPlot>[n_ptbins];
2324
2325 TCanvas *c_pfjet_res_e[n_etabins];
2326 TCanvas *c_pfjet_res_eta[n_ptbins];
2327
2328
2329 // loop over eta bins
2330 for (k = 0; k < etaVals.size()-1; k++)
2331 {
2332
2333 mg_pfjet_res_e[k] = new TMultiGraph("","");
2334 leg_pfjet_res_e[k] = new TLegend(0.40,0.70,0.90,0.90);
2335
2336 HistogramsCollection(&plots_pfjet_res_e[k], TMath::Log10(ptMin), TMath::Log10(ptMax), "pfjet");
2337 HistogramsCollection(&plots_cajet_res_e[k], TMath::Log10(ptMin), TMath::Log10(ptMax), "cajet");
2338
2339 GetJetsEres(&plots_pfjet_res_e[k], branchPFJet, branchGenJet, treeReaderJet, etaVals.at(k), etaVals.at(k+1));
2340 GetJetsEres(&plots_cajet_res_e[k], branchCaloJet, branchGenJet, treeReaderJet, etaVals.at(k), etaVals.at(k+1));
2341
2342 gr_pfjet_res_e[k] = EresGraph(&plots_pfjet_res_e[k]);
2343 gr_cajet_res_e[k] = EresGraph(&plots_cajet_res_e[k]);
2344
2345 s_etaMin = Form("%.1f",etaVals.at(k));
2346 s_etaMax = Form("%.1f",etaVals.at(k+1));
2347 s_eta = "anti-k_{T}, R = 0.4, "+ s_etaMin + " < | #eta | < " + s_etaMax;
2348
2349 leg_pfjet_res_e[k]->SetTextFont(132);
2350 leg_pfjet_res_e[k]->SetHeader(s_eta);
2351
2352 addResoGraph(mg_pfjet_res_e[k], &gr_cajet_res_e[k], leg_pfjet_res_e[k], markerStyles.at(0), colors.at(0), "Calorimeter Jets");
2353 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");
2354
2355 c_pfjet_res_e[k] = new TCanvas("","", 800, 600);
2356
2357 mg_pfjet_res_e[k]->Draw("APE");
2358 //DrawAxis(mg_pfjet_res_e[k], leg_pfjet_res_e[k], 10, ptMax, 0.5, 100, "E [GeV]", "(resolution in E)/E (%)", true, true);
2359 DrawAxis(mg_pfjet_res_e[k], leg_pfjet_res_e[k], 10, ptMax, 0.0, 30, "E [GeV]", "(resolution in E)/E (%)", true, false);
2360 leg_pfjet_res_e[k]->Draw();
2361 pave->Draw();
2362
2363 TString s_etarange = "eta_"+s_etaMin+"_"+s_etaMax+"_";
2364
2365 c_pfjet_res_e[k]->Print(pdfOutput,"pdf");
2366 c_pfjet_res_e[k]->Print(figPath+"img_pfjet_res_"+s_etarange+"e.pdf","pdf");
2367 c_pfjet_res_e[k]->Print(figPath+"img_pfjet_res_"+s_etarange+"e.png","png");
2368
2369 }
2370
2371
2372 // loop over eta bins
2373 for (k = 0; k < ptVals.size(); k++)
2374 {
2375
2376 mg_pfjet_res_eta[k] = new TMultiGraph("","");
2377 leg_pfjet_res_eta[k] = new TLegend(0.30,0.70,0.85,0.90);
2378
2379 HistogramsCollectionVsEta(&plots_pfjet_res_eta[k], etaMin, etaMax, 0.5*ptVals.at(k), 2.0*ptVals.at(k), "pfjet", 0.0, 2.0);
2380 HistogramsCollectionVsEta(&plots_cajet_res_eta[k], etaMin, etaMax, 0.5*ptVals.at(k), 2.0*ptVals.at(k), "cajet", 0.0, 2.0);
2381
2382 GetJetsEresVsEta(&plots_pfjet_res_eta[k], branchPFJet, branchGenJet, treeReaderJet, 0.5*ptVals.at(k), 2.0*ptVals.at(k));
2383 GetJetsEresVsEta(&plots_cajet_res_eta[k], branchCaloJet, branchGenJet, treeReaderJet, 0.5*ptVals.at(k), 2.0*ptVals.at(k));
2384
2385 gr_pfjet_res_eta[k] = EresGraphVsEta(&plots_pfjet_res_eta[k]);
2386 gr_cajet_res_eta[k] = EresGraphVsEta(&plots_cajet_res_eta[k]);
2387
2388 s_e = Form("anti-k_{T}, R = 0.4, jets, E = %.0f GeV",ptVals.at(k));
2389 if(ptVals.at(k) >= 1000.) s_e = Form("anti-k_{T}, R = 0.4, E = %.0f TeV",ptVals.at(k)/1000.);
2390
2391 leg_pfjet_res_eta[k]->SetTextFont(132);
2392 leg_pfjet_res_eta[k]->SetHeader(s_e);
2393
2394 addResoGraph(mg_pfjet_res_eta[k], &gr_cajet_res_eta[k], leg_pfjet_res_eta[k], markerStyles.at(0), colors.at(0), "Calorimeter Jets");
2395 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");
2396
2397 c_pfjet_res_eta[k] = new TCanvas("","", 800, 600);
2398
2399 mg_pfjet_res_eta[k]->Draw("APE");
2400 //DrawAxis(mg_pfjet_res_eta[k], leg_pfjet_res_eta[k], etaMin, etaMax, 0.1, 1000, "#eta", "(resolution in E)/E (%)", false, true);
2401 DrawAxis(mg_pfjet_res_eta[k], leg_pfjet_res_eta[k], etaMin, etaMax, 0.0, 50, "#eta", "(resolution in E)/E (%)", false, false);
2402 leg_pfjet_res_eta[k]->Draw();
2403 pave->Draw();
2404
2405 TString s_ptrange = Form("pt_%.0f_",ptVals.at(k));
2406
2407 c_pfjet_res_eta[k]->Print(pdfOutput,"pdf");
2408 c_pfjet_res_eta[k]->Print(figPath+"img_pfjet_res_"+s_ptrange+"eta.pdf","pdf");
2409 c_pfjet_res_eta[k]->Print(figPath+"img_pfjet_res_"+s_ptrange+"eta.png","png");
2410
2411 }
2412
2413
2414 /////////////////////
2415 // PF Missing ET ///
2416 /////////////////////
2417
2418 TMultiGraph *mg_met_res_ht = new TMultiGraph("","");
2419 TLegend *leg_met_res_ht = new TLegend(0.60,0.22,0.90,0.42);
2420
2421 std::vector<resolPlot> plots_pfmet, plots_camet;
2422
2423 HistogramsCollection(&plots_pfmet, TMath::Log10(ptMin), TMath::Log10(ptMax), "pfMET", -500, 500);
2424 HistogramsCollection(&plots_camet, TMath::Log10(ptMin), TMath::Log10(ptMax), "caMET", -500, 500);
2425
2426 GetMetres(&plots_pfmet, branchGenScalarHT, branchMet, branchPFJet, treeReaderJet);
2427 GetMetres(&plots_camet, branchGenScalarHT, branchCaloMet, branchCaloJet, treeReaderJet);
2428
2429 TGraphErrors gr_pfmet_res_ht = MetResGraph(&plots_pfmet, true);
2430 TGraphErrors gr_camet_res_ht = MetResGraph(&plots_camet, true);
2431
2432 addResoGraph(mg_met_res_ht, &gr_camet_res_ht, leg_met_res_ht, markerStyles.at(0), colors.at(0), "Calorimeter E_{T}^{miss}");
2433 addResoGraph(mg_met_res_ht, &gr_pfmet_res_ht, leg_met_res_ht, markerStyles.at(1), colors.at(1), "Particle-flow E_{T}^{miss}");
2434
2435 TCanvas *c_met_res_ht = new TCanvas("","", 800, 600);
2436
2437 mg_met_res_ht->Draw("APE");
2438 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);
2439
2440 leg_met_res_ht->Draw();
2441 pave->Draw();
2442 c_met_res_ht->Print(pdfOutput,"pdf");
2443 c_met_res_ht->Print(figPath+"img_met_res_ht.pdf","pdf");
2444 c_met_res_ht->Print(figPath+"img_met_res_ht.png","png");
2445
2446
2447 /////////////////////////////////////////
2448 // Electron Reconstruction Efficiency ///
2449 /////////////////////////////////////////
2450
2451 TMultiGraph *mg_recele_eff_pt = new TMultiGraph("","");
2452 TMultiGraph *mg_recele_eff_eta = new TMultiGraph("","");
2453
2454 TLegend *leg_recele_eff_pt = new TLegend(0.55,0.22,0.90,0.48);
2455 TLegend *leg_recele_eff_eta = new TLegend(0.55,0.22,0.90,0.48);
2456
2457 TGraphErrors *gr_recele_eff_pt = new TGraphErrors[n_etabins];
2458 TGraphErrors *gr_recele_eff_eta = new TGraphErrors[n_ptbins];
2459 TH1D* h_recele_eff_pt, *h_recele_eff_eta;
2460
2461 // loop over eta bins
2462 for (k = 0; k < etaVals.size()-1; k++)
2463 {
2464
2465 h_recele_eff_pt = GetEffPt<Electron>(branchElectron, branchParticleElectron, "Electron", 11, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderElectron);
2466 gr_recele_eff_pt[k] = TGraphErrors(h_recele_eff_pt);
2467
2468 s_etaMin = Form("%.1f",etaVals.at(k));
2469 s_etaMax = Form("%.1f",etaVals.at(k+1));
2470
2471 s_eta = "e^{ #pm} , " + s_etaMin + " < | #eta | < "+s_etaMax;
2472
2473 gr_recele_eff_pt[k].SetName("recEff_"+s_etaMin+"_"+s_etaMax);
2474
2475 addResoGraph(mg_recele_eff_pt, &gr_recele_eff_pt[k], leg_recele_eff_pt, markerStyles.at(k), colors.at(k), s_eta);
2476 }
2477
2478 // loop over pt
2479 for (k = 0; k < ptVals.size(); k++)
2480 {
2481 h_recele_eff_eta = GetEffEta<Electron>(branchElectron, branchParticleElectron, "Electron", 11, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderElectron);
2482 gr_recele_eff_eta[k] = TGraphErrors(h_recele_eff_eta);
2483
2484 s_pt = Form("e^{ #pm} , p_{T} = %.0f GeV",ptVals.at(k));
2485 if(ptVals.at(k) >= 1000.) s_pt = Form("e^{ #pm} , p_{T} = %.0f TeV",ptVals.at(k)/1000.);
2486
2487 addResoGraph(mg_recele_eff_eta, &gr_recele_eff_eta[k], leg_recele_eff_eta, markerStyles.at(k), colors.at(k), s_pt );
2488 }
2489
2490 TCanvas *c_recele_eff_pt = new TCanvas("","", 800, 600);
2491
2492 mg_recele_eff_pt->Draw("APE");
2493 DrawAxis(mg_recele_eff_pt, leg_recele_eff_pt, ptMin, ptMax, 0.0, 100, "p_{T} [GeV]", "reconstruction efficiency (%)", true, false);
2494 leg_recele_eff_pt->Draw();
2495 pave->Draw();
2496
2497 c_recele_eff_pt->Print(pdfOutput,"pdf");
2498 c_recele_eff_pt->Print(figPath+"img_recele_eff_pt.pdf","pdf");
2499 c_recele_eff_pt->Print(figPath+"img_recele_eff_pt.png","png");
2500
2501 TCanvas *c_recele_eff_eta = new TCanvas("","", 800, 600);
2502
2503 mg_recele_eff_eta->Draw("APE");
2504 DrawAxis(mg_recele_eff_eta, leg_recele_eff_eta, etaMin, etaMax, 0.0, 100, " #eta ", "reconstruction efficiency (%)", false, false);
2505 leg_recele_eff_eta->Draw();
2506 pave->Draw();
2507
2508 c_recele_eff_eta->Print(pdfOutput,"pdf");
2509 c_recele_eff_eta->Print(figPath+"img_recele_eff_eta.pdf","pdf");
2510 c_recele_eff_eta->Print(figPath+"img_recele_eff_eta.png","png");
2511
2512
2513 /////////////////////////////////////////
2514 // Muon Reconstruction Efficiency ///
2515 /////////////////////////////////////////
2516
2517 TMultiGraph *mg_recmu_eff_pt = new TMultiGraph("","");
2518 TMultiGraph *mg_recmu_eff_eta = new TMultiGraph("","");
2519
2520 TLegend *leg_recmu_eff_pt = new TLegend(0.55,0.22,0.90,0.48);
2521 TLegend *leg_recmu_eff_eta = new TLegend(0.55,0.22,0.90,0.48);
2522
2523 TGraphErrors *gr_recmu_eff_pt = new TGraphErrors[n_etabins];
2524 TGraphErrors *gr_recmu_eff_eta = new TGraphErrors[n_ptbins];
2525 TH1D* h_recmu_eff_pt, *h_recmu_eff_eta;
2526
2527 // loop over eta bins
2528 for (k = 0; k < etaVals.size()-1; k++)
2529 {
2530
2531 h_recmu_eff_pt = GetEffPt<Muon>(branchMuon, branchParticleMuon, "muon", 13, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderMuon);
2532 gr_recmu_eff_pt[k] = TGraphErrors(h_recmu_eff_pt);
2533
2534 s_etaMin = Form("%.1f",etaVals.at(k));
2535 s_etaMax = Form("%.1f",etaVals.at(k+1));
2536
2537 s_eta = "#mu^{ #pm} , " + s_etaMin + " < | #eta | < "+s_etaMax;
2538
2539 gr_recmu_eff_pt[k].SetName("recEff_"+s_etaMin+"_"+s_etaMax);
2540
2541 addResoGraph(mg_recmu_eff_pt, &gr_recmu_eff_pt[k], leg_recmu_eff_pt, markerStyles.at(k), colors.at(k), s_eta);
2542 }
2543
2544 // loop over pt
2545 for (k = 0; k < ptVals.size(); k++)
2546 {
2547 h_recmu_eff_eta = GetEffEta<Muon>(branchMuon, branchParticleMuon, "muon", 13, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderMuon);
2548 gr_recmu_eff_eta[k] = TGraphErrors(h_recmu_eff_eta);
2549
2550 s_pt = Form("#mu^{ #pm} , p_{T} = %.0f GeV",ptVals.at(k));
2551 if(ptVals.at(k) >= 1000.) s_pt = Form("#mu^{ #pm} , p_{T} = %.0f TeV",ptVals.at(k)/1000.);
2552
2553 addResoGraph(mg_recmu_eff_eta, &gr_recmu_eff_eta[k], leg_recmu_eff_eta, markerStyles.at(k), colors.at(k), s_pt );
2554 }
2555
2556 TCanvas *c_recmu_eff_pt = new TCanvas("","", 800, 600);
2557
2558 mg_recmu_eff_pt->Draw("APE");
2559 DrawAxis(mg_recmu_eff_pt, leg_recmu_eff_pt, ptMin, ptMax, 0.0, 100, "p_{T} [GeV]", "reconstruction efficiency (%)", true, false);
2560 leg_recmu_eff_pt->Draw();
2561 pave->Draw();
2562
2563 c_recmu_eff_pt->Print(pdfOutput,"pdf");
2564 c_recmu_eff_pt->Print(figPath+"img_recmu_eff_pt.pdf","pdf");
2565 c_recmu_eff_pt->Print(figPath+"img_recmu_eff_pt.png","png");
2566
2567 TCanvas *c_recmu_eff_eta = new TCanvas("","", 800, 600);
2568
2569 mg_recmu_eff_eta->Draw("APE");
2570 DrawAxis(mg_recmu_eff_eta, leg_recmu_eff_eta, etaMin, etaMax, 0.0, 100, " #eta ", "reconstruction efficiency (%)", false, false);
2571 leg_recmu_eff_eta->Draw();
2572 pave->Draw();
2573
2574 c_recmu_eff_eta->Print(pdfOutput,"pdf");
2575 c_recmu_eff_eta->Print(figPath+"img_recmu_eff_eta.pdf","pdf");
2576 c_recmu_eff_eta->Print(figPath+"img_recmu_eff_eta.png","png");
2577
2578
2579 /////////////////////////////////////////
2580 // Photon Reconstruction Efficiency ///
2581 /////////////////////////////////////////
2582
2583 TMultiGraph *mg_recpho_eff_pt = new TMultiGraph("","");
2584 TMultiGraph *mg_recpho_eff_eta = new TMultiGraph("","");
2585
2586 TLegend *leg_recpho_eff_pt = new TLegend(0.55,0.22,0.90,0.48);
2587 TLegend *leg_recpho_eff_eta = new TLegend(0.55,0.22,0.90,0.48);
2588
2589 TGraphErrors *gr_recpho_eff_pt = new TGraphErrors[n_etabins];
2590 TGraphErrors *gr_recpho_eff_eta = new TGraphErrors[n_ptbins];
2591 TH1D* h_recpho_eff_pt, *h_recpho_eff_eta;
2592
2593 // loop over eta bins
2594 for (k = 0; k < etaVals.size()-1; k++)
2595 {
2596
2597 h_recpho_eff_pt = GetEffPt<Photon>(branchPhoton, branchParticlePhoton, "Photon", 22, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderPhoton);
2598 gr_recpho_eff_pt[k] = TGraphErrors(h_recpho_eff_pt);
2599
2600 s_etaMin = Form("%.1f",etaVals.at(k));
2601 s_etaMax = Form("%.1f",etaVals.at(k+1));
2602
2603 s_eta = "#gamma , " + s_etaMin + " < | #eta | < "+s_etaMax;
2604
2605 gr_recpho_eff_pt[k].SetName("recEff_"+s_etaMin+"_"+s_etaMax);
2606
2607 addResoGraph(mg_recpho_eff_pt, &gr_recpho_eff_pt[k], leg_recpho_eff_pt, markerStyles.at(k), colors.at(k), s_eta);
2608 }
2609
2610 // loop over pt
2611 for (k = 0; k < ptVals.size(); k++)
2612 {
2613 h_recpho_eff_eta = GetEffEta<Photon>(branchPhoton, branchParticlePhoton, "Photon", 22, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderPhoton);
2614 gr_recpho_eff_eta[k] = TGraphErrors(h_recpho_eff_eta);
2615
2616 s_pt = Form("#gamma , p_{T} = %.0f GeV",ptVals.at(k));
2617 if(ptVals.at(k) >= 1000.) s_pt = Form("#gamma , p_{T} = %.0f TeV",ptVals.at(k)/1000.);
2618
2619 addResoGraph(mg_recpho_eff_eta, &gr_recpho_eff_eta[k], leg_recpho_eff_eta, markerStyles.at(k), colors.at(k), s_pt );
2620 }
2621
2622 TCanvas *c_recpho_eff_pt = new TCanvas("","", 800, 600);
2623
2624 mg_recpho_eff_pt->Draw("APE");
2625 DrawAxis(mg_recpho_eff_pt, leg_recpho_eff_pt, ptMin, ptMax, 0.0, 100, "p_{T} [GeV]", "reconstruction efficiency (%)", true, false);
2626 leg_recpho_eff_pt->Draw();
2627 pave->Draw();
2628
2629 c_recpho_eff_pt->Print(pdfOutput,"pdf");
2630 c_recpho_eff_pt->Print(figPath+"img_recpho_eff_pt.pdf","pdf");
2631 c_recpho_eff_pt->Print(figPath+"img_recpho_eff_pt.png","png");
2632
2633 TCanvas *c_recpho_eff_eta = new TCanvas("","", 800, 600);
2634
2635 mg_recpho_eff_eta->Draw("APE");
2636 DrawAxis(mg_recpho_eff_eta, leg_recpho_eff_eta, etaMin, etaMax, 0.0, 100, " #eta ", "reconstruction efficiency (%)", false, false);
2637 leg_recpho_eff_eta->Draw();
2638 pave->Draw();
2639
2640 c_recpho_eff_eta->Print(pdfOutput,"pdf");
2641 c_recpho_eff_eta->Print(figPath+"img_recpho_eff_eta.pdf","pdf");
2642 c_recpho_eff_eta->Print(figPath+"img_recpho_eff_eta.png","png");
2643
2644 /////////////////////////////////////////
2645 // B-jets Efficiency/ mistag rates ///
2646 /////////////////////////////////////////
2647
2648 TMultiGraph *mg_recbjet_eff_pt = new TMultiGraph("","");
2649 TMultiGraph *mg_recbjet_eff_eta = new TMultiGraph("","");
2650
2651 TLegend *leg_recbjet_eff_pt = new TLegend(0.50,0.22,0.90,0.48);
2652 TLegend *leg_recbjet_eff_eta = new TLegend(0.50,0.22,0.90,0.48);
2653
2654 TGraphErrors *gr_recbjet_eff_pt = new TGraphErrors[n_etabins];
2655 TGraphErrors *gr_recbjet_eff_eta = new TGraphErrors[n_ptbins];
2656 TH1D* h_recbjet_eff_pt, *h_recbjet_eff_eta;
2657
2658 // loop over eta bins
2659 for (k = 0; k < etaVals.size()-1; k++)
2660 {
2661
2662 h_recbjet_eff_pt = GetJetEffPt<Jet>(branchPFBJet, "BJet", 5, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderBJet);
2663 //h_recbjet_eff_pt = GetEffPt<Jet>(branchPFBJet, branchParticleBJet, "BJet", 5, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderBJet);
2664 gr_recbjet_eff_pt[k] = TGraphErrors(h_recbjet_eff_pt);
2665
2666 s_etaMin = Form("%.1f",etaVals.at(k));
2667 s_etaMax = Form("%.1f",etaVals.at(k+1));
2668
2669 s_eta = "b-jet , " + s_etaMin + " < | #eta | < "+s_etaMax;
2670
2671 gr_recbjet_eff_pt[k].SetName("recEff_"+s_etaMin+"_"+s_etaMax);
2672
2673 addResoGraph(mg_recbjet_eff_pt, &gr_recbjet_eff_pt[k], leg_recbjet_eff_pt, markerStyles.at(k), colors.at(k), s_eta);
2674 }
2675
2676 // loop over pt
2677 for (k = 0; k < ptVals.size(); k++)
2678 {
2679 h_recbjet_eff_eta = GetJetEffEta<Jet>(branchPFBJet, "BJet", 5, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderBJet);
2680 //h_recbjet_eff_eta = GetEffEta<Jet>(branchPFBJet, branchParticleBJet, "BJet", 5, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderBJet);
2681 gr_recbjet_eff_eta[k] = TGraphErrors(h_recbjet_eff_eta);
2682
2683 s_pt = Form("b-jet , p_{T} = %.0f GeV",ptVals.at(k));
2684 if(ptVals.at(k) >= 1000.) s_pt = Form("b-jet , p_{T} = %.0f TeV",ptVals.at(k)/1000.);
2685
2686 addResoGraph(mg_recbjet_eff_eta, &gr_recbjet_eff_eta[k], leg_recbjet_eff_eta, markerStyles.at(k), colors.at(k), s_pt );
2687 }
2688
2689 TCanvas *c_recbjet_eff_pt = new TCanvas("","", 800, 600);
2690
2691 mg_recbjet_eff_pt->Draw("APE");
2692 DrawAxis(mg_recbjet_eff_pt, leg_recbjet_eff_pt, ptMin, ptMax, 0.0, 100, "p_{T} [GeV]", "b - tag efficiency (%)", true, false);
2693 leg_recbjet_eff_pt->Draw();
2694 pave->Draw();
2695
2696 c_recbjet_eff_pt->Print(pdfOutput,"pdf");
2697 c_recbjet_eff_pt->Print(figPath+"img_recbjet_eff_pt.pdf","pdf");
2698 c_recbjet_eff_pt->Print(figPath+"img_recbjet_eff_pt.png","png");
2699
2700 TCanvas *c_recbjet_eff_eta = new TCanvas("","", 800, 600);
2701
2702 mg_recbjet_eff_eta->Draw("APE");
2703 DrawAxis(mg_recbjet_eff_eta, leg_recbjet_eff_eta, etaMin, etaMax, 0.0, 100, " #eta ", "b - tag efficiency (%)", false, false);
2704 leg_recbjet_eff_eta->Draw();
2705 pave->Draw();
2706
2707 c_recbjet_eff_eta->Print(pdfOutput,"pdf");
2708 c_recbjet_eff_eta->Print(figPath+"img_recbjet_eff_eta.pdf","pdf");
2709 c_recbjet_eff_eta->Print(figPath+"img_recbjet_eff_eta.png","png");
2710
2711 // ------ c - mistag ------
2712
2713 TMultiGraph *mg_recbjet_cmis_pt = new TMultiGraph("","");
2714 TMultiGraph *mg_recbjet_cmis_eta = new TMultiGraph("","");
2715
2716 TLegend *leg_recbjet_cmis_pt = new TLegend(0.50,0.64,0.90,0.90);
2717 TLegend *leg_recbjet_cmis_eta = new TLegend(0.50,0.64,0.90,0.90);
2718
2719 TGraphErrors *gr_recbjet_cmis_pt = new TGraphErrors[n_etabins];
2720 TGraphErrors *gr_recbjet_cmis_eta = new TGraphErrors[n_ptbins];
2721 TH1D* h_recbjet_cmis_pt, *h_recbjet_cmis_eta;
2722
2723 // loop over eta bins
2724 for (k = 0; k < etaVals.size()-1; k++)
2725 {
2726
2727 h_recbjet_cmis_pt = GetJetEffPt<Jet>(branchPFCJet, "CJet", 4, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderCJet);
2728 //h_recbjet_cmis_pt = GetEffPt<Jet>(branchPFCJet, branchParticleCJet, "CJet", 4, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderCJet);
2729 gr_recbjet_cmis_pt[k] = TGraphErrors(h_recbjet_cmis_pt);
2730
2731 s_etaMin = Form("%.1f",etaVals.at(k));
2732 s_etaMax = Form("%.1f",etaVals.at(k+1));
2733
2734 s_eta = "c-jet , " + s_etaMin + " < | #eta | < "+s_etaMax;
2735
2736 gr_recbjet_cmis_pt[k].SetName("recEff_"+s_etaMin+"_"+s_etaMax);
2737
2738 addResoGraph(mg_recbjet_cmis_pt, &gr_recbjet_cmis_pt[k], leg_recbjet_cmis_pt, markerStyles.at(k), colors.at(k), s_eta);
2739 }
2740
2741 // loop over pt
2742 for (k = 0; k < ptVals.size(); k++)
2743 {
2744 h_recbjet_cmis_eta = GetJetEffEta<Jet>(branchPFCJet, "CJet", 4, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderCJet);
2745 //h_recbjet_cmis_eta = GetEffEta<Jet>(branchPFCJet, branchParticleCJet, "CJet", 4, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderCJet);
2746 gr_recbjet_cmis_eta[k] = TGraphErrors(h_recbjet_cmis_eta);
2747
2748 s_pt = Form("c-jet , p_{T} = %.0f GeV",ptVals.at(k));
2749 if(ptVals.at(k) >= 1000.) s_pt = Form("c-jet , p_{T} = %.0f TeV",ptVals.at(k)/1000.);
2750
2751 addResoGraph(mg_recbjet_cmis_eta, &gr_recbjet_cmis_eta[k], leg_recbjet_cmis_eta, markerStyles.at(k), colors.at(k), s_pt );
2752 }
2753
2754 TCanvas *c_recbjet_cmis_pt = new TCanvas("","", 800, 600);
2755
2756 mg_recbjet_cmis_pt->Draw("APE");
2757 DrawAxis(mg_recbjet_cmis_pt, leg_recbjet_cmis_pt, ptMin, ptMax, 0.0, 20, "p_{T} [GeV]", "c - mistag rate (%)", true, false);
2758 leg_recbjet_cmis_pt->Draw();
2759 pave->Draw();
2760
2761 c_recbjet_cmis_pt->Print(pdfOutput,"pdf");
2762 c_recbjet_cmis_pt->Print(figPath+"img_recbjet_cmis_pt.pdf","pdf");
2763 c_recbjet_cmis_pt->Print(figPath+"img_recbjet_cmis_pt.png","png");
2764
2765 TCanvas *c_recbjet_cmis_eta = new TCanvas("","", 800, 600);
2766
2767 mg_recbjet_cmis_eta->Draw("APE");
2768 DrawAxis(mg_recbjet_cmis_eta, leg_recbjet_cmis_eta, etaMin, etaMax, 0.0, 20, " #eta ", "c - mistag rate (%)", false, false);
2769 leg_recbjet_cmis_eta->Draw();
2770 pave->Draw();
2771
2772 c_recbjet_cmis_eta->Print(pdfOutput,"pdf");
2773 c_recbjet_cmis_eta->Print(figPath+"img_recbjet_cmis_eta.pdf","pdf");
2774 c_recbjet_cmis_eta->Print(figPath+"img_recbjet_cmis_eta.png","png");
2775
2776 // ------ light - mistag ------
2777
2778 TMultiGraph *mg_recbjet_lmis_pt = new TMultiGraph("","");
2779 TMultiGraph *mg_recbjet_lmis_eta = new TMultiGraph("","");
2780
2781 TLegend *leg_recbjet_lmis_pt = new TLegend(0.50,0.64,0.90,0.90);
2782 TLegend *leg_recbjet_lmis_eta = new TLegend(0.50,0.64,0.90,0.90);
2783
2784 TGraphErrors *gr_recbjet_lmis_pt = new TGraphErrors[n_etabins];
2785 TGraphErrors *gr_recbjet_lmis_eta = new TGraphErrors[n_ptbins];
2786 TH1D* h_recbjet_lmis_pt, *h_recbjet_lmis_eta;
2787
2788 // loop over eta bins
2789 for (k = 0; k < etaVals.size()-1; k++)
2790 {
2791
2792 h_recbjet_lmis_pt = GetJetEffPt<Jet>(branchJet, "Jet", 1, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderJet);
2793 //h_recbjet_lmis_pt = GetEffPt<Jet>(branchJet, branchParticleJet, "Jet", 1, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderJet);
2794 gr_recbjet_lmis_pt[k] = TGraphErrors(h_recbjet_lmis_pt);
2795
2796 s_etaMin = Form("%.1f",etaVals.at(k));
2797 s_etaMax = Form("%.1f",etaVals.at(k+1));
2798
2799 s_eta = "uds-jet , " + s_etaMin + " < | #eta | < "+s_etaMax;
2800
2801 gr_recbjet_lmis_pt[k].SetName("recEff_"+s_etaMin+"_"+s_etaMax);
2802
2803 addResoGraph(mg_recbjet_lmis_pt, &gr_recbjet_lmis_pt[k], leg_recbjet_lmis_pt, markerStyles.at(k), colors.at(k), s_eta);
2804 }
2805
2806 // loop over pt
2807 for (k = 0; k < ptVals.size(); k++)
2808 {
2809 h_recbjet_lmis_eta = GetJetEffEta<Jet>(branchJet, "Jet", 1, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderJet);
2810 //h_recbjet_lmis_eta = GetEffEta<Jet>(branchJet, branchParticleJet, "Jet", 1, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderJet);
2811 gr_recbjet_lmis_eta[k] = TGraphErrors(h_recbjet_lmis_eta);
2812
2813 s_pt = Form("uds-jet , p_{T} = %.0f GeV",ptVals.at(k));
2814 if(ptVals.at(k) >= 1000.) s_pt = Form("uds-jet , p_{T} = %.0f TeV",ptVals.at(k)/1000.);
2815
2816 addResoGraph(mg_recbjet_lmis_eta, &gr_recbjet_lmis_eta[k], leg_recbjet_lmis_eta, markerStyles.at(k), colors.at(k), s_pt );
2817 }
2818
2819 TCanvas *c_recbjet_lmis_pt = new TCanvas("","", 800, 600);
2820
2821 mg_recbjet_lmis_pt->Draw("APE");
2822
2823 DrawAxis(mg_recbjet_lmis_pt, leg_recbjet_lmis_pt, ptMin, ptMax, 0.0, 1.0, "p_{T} [GeV]", "light - mistag rate (%)", true, false);
2824
2825 leg_recbjet_lmis_pt->Draw();
2826 pave->Draw();
2827
2828 c_recbjet_lmis_pt->Print(pdfOutput,"pdf");
2829 c_recbjet_lmis_pt->Print(figPath+"img_recbjet_lmis_pt.pdf","pdf");
2830 c_recbjet_lmis_pt->Print(figPath+"img_recbjet_lmis_pt.png","png");
2831
2832 TCanvas *c_recbjet_lmis_eta = new TCanvas("","", 800, 600);
2833
2834 mg_recbjet_lmis_eta->Draw("APE");
2835 DrawAxis(mg_recbjet_lmis_eta, leg_recbjet_lmis_eta, etaMin, etaMax, 0.0, 1.0, " #eta ", "light - mistag rate (%)", false, false);
2836 leg_recbjet_lmis_eta->Draw();
2837 pave->Draw();
2838
2839 c_recbjet_lmis_eta->Print(pdfOutput,"pdf");
2840 c_recbjet_lmis_eta->Print(figPath+"img_recbjet_lmis_eta.pdf","pdf");
2841 c_recbjet_lmis_eta->Print(figPath+"img_recbjet_lmis_eta.png","png");
2842
2843
2844 ///////////////////////////////////////////
2845 // tau-jets Efficiency/ mistag rates ///
2846 ///////////////////////////////////////////
2847
2848 TMultiGraph *mg_rectaujet_eff_pt = new TMultiGraph("","");
2849 TMultiGraph *mg_rectaujet_eff_eta = new TMultiGraph("","");
2850
2851 TLegend *leg_rectaujet_eff_pt = new TLegend(0.50,0.22,0.90,0.48);
2852 TLegend *leg_rectaujet_eff_eta = new TLegend(0.50,0.22,0.90,0.48);
2853
2854 TGraphErrors *gr_rectaujet_eff_pt = new TGraphErrors[n_etabins];
2855 TGraphErrors *gr_rectaujet_eff_eta = new TGraphErrors[n_ptbins];
2856 TH1D* h_rectaujet_eff_pt, *h_rectaujet_eff_eta;
2857
2858 // loop over eta bins
2859 for (k = 0; k < etaVals.size()-1; k++)
2860 {
2861
2862 h_rectaujet_eff_pt = GetTauEffPt<Jet>(branchPFTauJet, branchParticleTauJet, "TauJet", 15, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderTauJet);
2863 gr_rectaujet_eff_pt[k] = TGraphErrors(h_rectaujet_eff_pt);
2864
2865 s_etaMin = Form("%.1f",etaVals.at(k));
2866 s_etaMax = Form("%.1f",etaVals.at(k+1));
2867
2868 s_eta = "#tau-jet , " + s_etaMin + " < | #eta | < "+s_etaMax;
2869
2870 gr_rectaujet_eff_pt[k].SetName("recEff_"+s_etaMin+"_"+s_etaMax);
2871
2872 addResoGraph(mg_rectaujet_eff_pt, &gr_rectaujet_eff_pt[k], leg_rectaujet_eff_pt, markerStyles.at(k), colors.at(k), s_eta);
2873 }
2874
2875 // loop over pt
2876 for (k = 0; k < ptVals.size(); k++)
2877 {
2878 h_rectaujet_eff_eta = GetTauEffEta<Jet>(branchPFTauJet, branchParticleTauJet, "TauJet", 15, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderTauJet);
2879 gr_rectaujet_eff_eta[k] = TGraphErrors(h_rectaujet_eff_eta);
2880
2881 s_pt = Form("#tau-jet , p_{T} = %.0f GeV",ptVals.at(k));
2882 if(ptVals.at(k) >= 1000.) s_pt = Form("#tau-jet , p_{T} = %.0f TeV",ptVals.at(k)/1000.);
2883
2884 addResoGraph(mg_rectaujet_eff_eta, &gr_rectaujet_eff_eta[k], leg_rectaujet_eff_eta, markerStyles.at(k), colors.at(k), s_pt );
2885 }
2886
2887
2888 TCanvas *c_rectaujet_eff_pt = new TCanvas("","", 800, 600);
2889
2890 mg_rectaujet_eff_pt->Draw("APE");
2891 DrawAxis(mg_rectaujet_eff_pt, leg_rectaujet_eff_pt, ptMin, ptMax, 0.0, 100, "p_{T} [GeV]", "#tau - tag efficiency (%)", true, false);
2892 leg_rectaujet_eff_pt->Draw();
2893 pave->Draw();
2894
2895 c_rectaujet_eff_pt->Print(pdfOutput,"pdf");
2896 c_rectaujet_eff_pt->Print(figPath+"img_rectaujet_eff_pt.pdf","pdf");
2897 c_rectaujet_eff_pt->Print(figPath+"img_rectaujet_eff_pt.png","png");
2898
2899 TCanvas *c_rectaujet_eff_eta = new TCanvas("","", 800, 600);
2900
2901 mg_rectaujet_eff_eta->Draw("APE");
2902 DrawAxis(mg_rectaujet_eff_eta, leg_rectaujet_eff_eta, etaMin, etaMax, 0.0, 100., " #eta ", "#tau - tag efficiency (%)", false, false);
2903 leg_rectaujet_eff_eta->Draw();
2904 pave->Draw();
2905
2906 c_rectaujet_eff_eta->Print(pdfOutput,"pdf");
2907 c_rectaujet_eff_eta->Print(figPath+"img_rectaujet_eff_eta.pdf","pdf");
2908 c_rectaujet_eff_eta->Print(figPath+"img_rectaujet_eff_eta.png","png");
2909
2910
2911 //--------------- tau mistag rate ----------
2912
2913 TMultiGraph *mg_rectaujet_mis_pt = new TMultiGraph("","");
2914 TMultiGraph *mg_rectaujet_mis_eta = new TMultiGraph("","");
2915
2916 TLegend *leg_rectaujet_mis_pt = new TLegend(0.50,0.64,0.90,0.90);
2917 TLegend *leg_rectaujet_mis_eta = new TLegend(0.50,0.64,0.90,0.90);
2918
2919 TGraphErrors *gr_rectaujet_mis_pt = new TGraphErrors[n_etabins];
2920 TGraphErrors *gr_rectaujet_mis_eta = new TGraphErrors[n_ptbins];
2921 TH1D* h_rectaujet_mis_pt, *h_rectaujet_mis_eta;
2922
2923 // loop over eta bins
2924 for (k = 0; k < etaVals.size()-1; k++)
2925 {
2926
2927 h_rectaujet_mis_pt = GetTauEffPt<Jet>(branchJet, branchParticleJet, "TauJet", 1, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderJet);
2928 gr_rectaujet_mis_pt[k] = TGraphErrors(h_rectaujet_mis_pt);
2929
2930 s_etaMin = Form("%.1f",etaVals.at(k));
2931 s_etaMax = Form("%.1f",etaVals.at(k+1));
2932
2933 s_eta = "uds-jet , " + s_etaMin + " < | #eta | < "+s_etaMax;
2934
2935 gr_rectaujet_mis_pt[k].SetName("recEff_"+s_etaMin+"_"+s_etaMax);
2936
2937 addResoGraph(mg_rectaujet_mis_pt, &gr_rectaujet_mis_pt[k], leg_rectaujet_mis_pt, markerStyles.at(k), colors.at(k), s_eta);
2938 }
2939
2940 // loop over pt
2941 for (k = 0; k < ptVals.size(); k++)
2942 {
2943 h_rectaujet_mis_eta = GetTauEffEta<Jet>(branchJet, branchParticleJet, "TauJet", 1, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderJet);
2944 gr_rectaujet_mis_eta[k] = TGraphErrors(h_rectaujet_mis_eta);
2945
2946 s_pt = Form("uds-jet , p_{T} = %.0f GeV",ptVals.at(k));
2947 if(ptVals.at(k) >= 1000.) s_pt = Form("uds-jet , p_{T} = %.0f TeV",ptVals.at(k)/1000.);
2948
2949 addResoGraph(mg_rectaujet_mis_eta, &gr_rectaujet_mis_eta[k], leg_rectaujet_mis_eta, markerStyles.at(k), colors.at(k), s_pt );
2950 }
2951
2952 TCanvas *c_rectaujet_mis_pt = new TCanvas("","", 800, 600);
2953
2954 mg_rectaujet_mis_pt->Draw("APE");
2955 DrawAxis(mg_rectaujet_mis_pt, leg_rectaujet_mis_pt, ptMin, ptMax, 0.0, 5., "p_{T} [GeV]", "#tau - mistag(%)", true, false);
2956 leg_rectaujet_mis_pt->Draw();
2957 pave->Draw();
2958
2959 c_rectaujet_mis_pt->Print(pdfOutput,"pdf");
2960 c_rectaujet_mis_pt->Print(figPath+"img_rectaujet_mis_pt.pdf","pdf");
2961 c_rectaujet_mis_pt->Print(figPath+"img_rectaujet_mis_pt.png","png");
2962
2963 TCanvas *c_rectaujet_mis_eta = new TCanvas("","", 800, 600);
2964
2965 mg_rectaujet_mis_eta->Draw("APE");
2966 DrawAxis(mg_rectaujet_mis_eta, leg_rectaujet_mis_eta, etaMin, etaMax, 0.0, 5., " #eta ", "#tau - mistag (%)", false, false);
2967 leg_rectaujet_mis_eta->Draw();
2968 pave->Draw();
2969
2970 c_rectaujet_mis_eta->Print(pdfOutput+")","pdf");
2971 c_rectaujet_mis_eta->Print(figPath+"img_rectaujet_mis_eta.pdf","pdf");
2972 c_rectaujet_mis_eta->Print(figPath+"img_rectaujet_mis_eta.png","png");
2973
2974
2975// ---- store resolution histograms in the output (for leave efficiencies out) ---
2976
2977
2978 TFile *fout = new TFile(outputFile,"recreate");
2979
2980 for (int bin = 0; bin < Nbins; bin++)
2981 {
2982
2983 for (k = 0; k < etaVals.size()-1; k++)
2984 {
2985 plots_trkpi_res_pt[k].at(bin).resolHist->Write();
2986 plots_trkele_res_pt[k].at(bin).resolHist->Write();
2987 plots_trkmu_res_pt[k].at(bin).resolHist->Write();
2988 plots_ecal_res_e[k].at(bin).resolHist->Write();
2989 plots_hcal_res_e[k].at(bin).resolHist->Write();
2990 plots_pfele_res_e[k].at(bin).resolHist->Write();
2991 plots_pfpi_res_e[k].at(bin).resolHist->Write();
2992 plots_pfjet_res_e[k].at(bin).resolHist->Write();
2993 plots_cajet_res_e[k].at(bin).resolHist->Write();
2994
2995 }
2996
2997 for (k = 0; k < ptVals.size(); k++)
2998 {
2999 plots_trkpi_res_eta[k].at(bin).resolHist->Write();
3000 plots_trkele_res_eta[k].at(bin).resolHist->Write();
3001 plots_trkmu_res_eta[k].at(bin).resolHist->Write();
3002 plots_ecal_res_eta[k].at(bin).resolHist->Write();
3003 plots_hcal_res_eta[k].at(bin).resolHist->Write();
3004 plots_pfele_res_eta[k].at(bin).resolHist->Write();
3005 plots_pfpi_res_eta[k].at(bin).resolHist->Write();
3006 plots_pfjet_res_eta[k].at(bin).resolHist->Write();
3007 plots_cajet_res_eta[k].at(bin).resolHist->Write();
3008
3009 }
3010
3011 plots_pfmet.at(bin).resolHist->Write();
3012 plots_camet.at(bin).resolHist->Write();
3013
3014 }
3015
3016 fout->Write();
3017
3018
3019 cout << "** Exiting..." << endl;
3020
3021 delete treeReaderElectron;
3022 delete treeReaderMuon;
3023 delete treeReaderPhoton;
3024 delete treeReaderJet;
3025 delete treeReaderBJet;
3026 delete treeReaderTauJet;
3027 delete chainElectron;
3028 delete chainMuon;
3029 delete chainPhoton;
3030 delete chainJet;
3031 delete chainBJet;
3032 delete chainTauJet;
3033}
3034
3035//------------------------------------------------------------------------------
3036
3037int main(int argc, char *argv[])
3038{
3039 char *appName = "DelphesValidation";
3040
3041 if(argc != 12)
3042 {
3043 cout << " Usage: " << appName << " input_file_electron input_file_muon input_file_photon input_file_jet input_file_bjet input_file_taujet output_file version" << endl;
3044 cout << " input_file_pion - input file in ROOT format ('Delphes' tree)," << endl;
3045 cout << " input_file_electron - input file in ROOT format ('Delphes' tree)," << endl;
3046 cout << " input_file_muon - input file in ROOT format ('Delphes' tree)," << endl;
3047 cout << " input_file_photon - input file in ROOT format ('Delphes' tree)," << endl;
3048 cout << " input_file_neutralhadron - input file in ROOT format ('Delphes' tree)," << endl;
3049 cout << " input_file_jet - input file in ROOT format ('Delphes' tree)," << endl;
3050 cout << " input_file_bjet - input file in ROOT format ('Delphes' tree)," << endl;
3051 cout << " input_file_cjet - input file in ROOT format ('Delphes' tree)," << endl;
3052 cout << " input_file_taujet - input file in ROOT format ('Delphes' tree)," << endl;
3053 cout << " output_file - output file in ROOT format," << endl;
3054 cout << " version - Delphes version" << endl;
3055
3056 return 1;
3057 }
3058
3059 gROOT->SetBatch();
3060
3061 int appargc = 1;
3062 char *appargv[] = {appName};
3063 TApplication app(appName, &appargc, appargv);
3064
3065 DelphesValidation(argv[1], argv[2], argv[3], argv[4], argv[5], argv[6], argv[7], argv[8], argv[9], argv[10], argv[11]);
3066}
3067
3068
3069
Note: See TracBrowser for help on using the repository browser.