Fork me on GitHub

source: git/examples/Validation.cpp@ 6934875

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

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