Fork me on GitHub

source: git/examples/Validation.cpp@ b4ec6ac

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

fixed style in validation plots
  • Property mode set to 100644
File size: 98.4 KB
Line 
1/*
2 * Delphes: a framework for fast simulation of a generic collider experiment
3 * Copyright (C) 2012-2014 Universite catholique de Louvain (UCL), Belgium
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19
20#include <iostream>
21#include <utility>
22#include <vector>
23#include <typeinfo>
24
25#include "TROOT.h"
26#include "TSystem.h"
27#include "TApplication.h"
28
29#include "TString.h"
30
31#include "TH1.h"
32#include "TH2.h"
33#include "TMath.h"
34#include "TStyle.h"
35#include "TGraph.h"
36#include "TCanvas.h"
37#include "THStack.h"
38#include "TLegend.h"
39#include "TPaveText.h"
40#include "TClonesArray.h"
41#include "TLorentzVector.h"
42#include "TGraphErrors.h"
43#include "TMultiGraph.h"
44
45#include "classes/DelphesClasses.h"
46
47#include "ExRootAnalysis/ExRootTreeReader.h"
48#include "ExRootAnalysis/ExRootTreeWriter.h"
49#include "ExRootAnalysis/ExRootTreeBranch.h"
50#include "ExRootAnalysis/ExRootResult.h"
51#include "ExRootAnalysis/ExRootUtilities.h"
52
53using namespace std;
54
55//------------------------------------------------------------------------------
56
57static const int Nbins = 50;
58
59int objStyle = 1;
60int trackStyle = 7;
61int towerStyle = 3;
62
63Color_t objColor = kBlack;
64Color_t trackColor = kBlack;
65Color_t towerColor = kBlack;
66
67double effLegXmin = 0.22;
68double effLegXmax = 0.7;
69double effLegYmin = 0.22;
70double effLegYmax = 0.5;
71
72double resLegXmin = 0.62;
73double resLegXmax = 0.9;
74double resLegYmin = 0.52;
75double resLegYmax = 0.85;
76
77double topLeftLegXmin = 0.22;
78double topLeftLegXmax = 0.7;
79double topLeftLegYmin = 0.52;
80double topLeftLegYmax = 0.85;
81
82unsigned int k;
83
84
85struct resolPlot
86{
87 TH1* resolHist;
88 double ptmin;
89 double ptmax;
90 double etamin;
91 double etamax;
92 double xmin;
93 double xmax;
94 TString obj;
95
96 resolPlot();
97 resolPlot(double ptdown, double ptup, TString object);
98 resolPlot(double etadown, double etaup, double ptdown, double ptup, TString object);
99 void set(double ptdown, double ptup, TString object, double xmin = 0, double xmax = 2);
100 void set(double etadown, double etaup, double ptdown, double ptup, TString object, double xmin = 0, double xmax = 2);
101 void print(){std::cout << ptmin << std::endl;}
102};
103
104
105resolPlot::resolPlot()
106{
107}
108
109resolPlot::resolPlot(double ptdown, double ptup, TString object)
110{
111 this->set(ptdown,ptup,object);
112}
113
114resolPlot::resolPlot(double etadown, double etaup, double ptdown, double ptup, TString object)
115{
116 this->set(etadown, etaup, ptdown, ptup, object);
117}
118
119void resolPlot::set(double ptdown, double ptup, TString object, double xmin, double xmax)
120{
121 ptmin = ptdown;
122 ptmax = ptup;
123 obj = object;
124
125 resolHist = new TH1D(obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax), obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax), 1000, xmin, xmax);
126}
127
128void resolPlot::set(double etadown, double etaup, double ptdown, double ptup, TString object, double xmin, double xmax)
129{
130 etamin = etadown;
131 etamax = etaup;
132 ptmin = ptdown;
133 ptmax = ptup;
134 obj = object;
135
136 resolHist = new TH1D(obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_"+Form("%4.2f",etamin)+"_"+Form("%4.2f",etamax), obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_"+Form("%4.2f",etamin)+"_"+Form("%4.2f",etamax), 1000, xmin, xmax);
137}
138
139
140
141void HistogramsCollection(std::vector<resolPlot> *histos, double ptmin, double ptmax, TString obj, double xmin = 0, double xmax = 2)
142{
143 double width;
144 double ptdown;
145 double ptup;
146 resolPlot ptemp;
147
148 for (int i = 0; i < Nbins; i++)
149 {
150 width = (ptmax - ptmin) / Nbins;
151 ptdown = TMath::Power(10,ptmin + i * width );
152 ptup = TMath::Power(10,ptmin + (i+1) * width );
153 ptemp.set(ptdown, ptup, obj, xmin, xmax);
154 histos->push_back(ptemp);
155 }
156}
157
158
159void HistogramsCollectionVsEta(std::vector<resolPlot> *histos, double etamin, double etamax, double ptmin, double ptmax, TString obj, double xmin = 0, double xmax = 2)
160{
161 resolPlot ptemp;
162 double width;
163 double etadown;
164 double etaup;
165
166 for (int i = 0; i < Nbins; i++)
167 {
168 width = (etamax - etamin) / Nbins;
169 etadown = etamin + i * width;
170 etaup = etamin + (i+1) * width;
171
172 ptemp.set(etadown, etaup, ptmin, ptmax, obj, xmin, xmax);
173 histos->push_back(ptemp);
174 }
175}
176
177
178//------------------------------------------------------------------------------
179
180class ExRootResult;
181class ExRootTreeReader;
182
183//------------------------------------------------------------------------------
184
185void BinLogX(TH1*h)
186{
187 TAxis *axis = h->GetXaxis();
188 int bins = axis->GetNbins();
189
190 Axis_t from = axis->GetXmin();
191 Axis_t to = axis->GetXmax();
192 Axis_t width = (to - from) / bins;
193 Axis_t *new_bins = new Axis_t[bins + 1];
194
195 for (int i = 0; i <= bins; i++)
196 {
197 new_bins[i] = TMath::Power(10, from + i * width);
198 }
199 axis->Set(bins, new_bins);
200 delete new_bins;
201}
202
203
204//------------------------------------------------------------------------------
205
206template<typename T>
207TH1D* GetEffPt(TClonesArray *branchReco, TClonesArray *branchParticle, TString name, int pdgID, double ptmin, double ptmax, double etamin, double etamax, ExRootTreeReader *treeReader)
208{
209
210 cout << "** Computing Efficiency of reconstructing "<< branchReco->GetName() << " induced by " << branchParticle->GetName() << " with PID " << pdgID << endl;
211
212 Long64_t allEntries = treeReader->GetEntries();
213
214 GenParticle *particle;
215 T *recoObj;
216
217 TLorentzVector recoMomentum, genMomentum, bestRecoMomentum;
218
219 Float_t deltaR;
220 Float_t pt, eta;
221 Long64_t entry;
222
223 Int_t i, j;
224
225 TH1D *histGenPt = new TH1D(name+" gen spectra Pt",name+" gen spectra cen", Nbins, TMath::Log10(ptmin), TMath::Log10(ptmax));
226 TH1D *histRecoPt = new TH1D(name+" reco spectra Pt",name+" reco spectra cen", Nbins, TMath::Log10(ptmin), TMath::Log10(ptmax));
227
228 histGenPt->SetDirectory(0);
229 histRecoPt->SetDirectory(0);
230
231 BinLogX(histGenPt);
232 BinLogX(histRecoPt);
233
234 // Loop over all events
235 for(entry = 0; entry < allEntries; ++entry)
236 {
237 // Load selected branches with data from specified event
238 treeReader->ReadEntry(entry);
239
240 // Loop over all generated particle in event
241 for(i = 0; i < branchParticle->GetEntriesFast(); ++i)
242 {
243
244 particle = (GenParticle*) branchParticle->At(i);
245 genMomentum = particle->P4();
246
247 deltaR = 999;
248
249 pt = genMomentum.Pt();
250 eta = TMath::Abs(genMomentum.Eta());
251
252 if(eta > etamax || eta < etamin ) continue;
253
254 if (particle->PID == pdgID && genMomentum.Pt() > ptmin && genMomentum.Pt() < ptmax )
255 {
256 // Loop over all reco object in event
257 for(j = 0; j < branchReco->GetEntriesFast(); ++j)
258 {
259 recoObj = (T*)branchReco->At(j);
260 recoMomentum = recoObj->P4();
261 // this is simply to avoid warnings from initial state particle
262 // having infite rapidity ...
263 //if(Momentum.Px() == 0 && genMomentum.Py() == 0) continue;
264
265 // 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
1268 //gPad->SetGridx();
1269 //gPad->SetGridy();
1270 gPad->SetBottomMargin(0.2);
1271 gPad->SetLeftMargin(0.2);
1272 gPad->Modified();
1273 gPad->Update();
1274
1275}
1276
1277
1278void Validation(const char *inputFilePion,
1279 const char *inputFileElectron,
1280 const char *inputFileMuon,
1281 const char *inputFilePhoton,
1282 const char *inputFileNeutralHadron,
1283 const char *inputFileJet,
1284 const char *inputFileBJet,
1285 const char *inputFileCJet,
1286 const char *inputFileTauJet,
1287 const char *outputFile,
1288 const char *version)
1289{
1290
1291 TChain *chainPion = new TChain("Delphes");
1292 chainPion->Add(inputFilePion);
1293 ExRootTreeReader *treeReaderPion = new ExRootTreeReader(chainPion);
1294
1295 TChain *chainElectron = new TChain("Delphes");
1296 chainElectron->Add(inputFileElectron);
1297 ExRootTreeReader *treeReaderElectron = new ExRootTreeReader(chainElectron);
1298
1299 TChain *chainMuon = new TChain("Delphes");
1300 chainMuon->Add(inputFileMuon);
1301 ExRootTreeReader *treeReaderMuon = new ExRootTreeReader(chainMuon);
1302
1303 TChain *chainPhoton = new TChain("Delphes");
1304 chainPhoton->Add(inputFilePhoton);
1305 ExRootTreeReader *treeReaderPhoton = new ExRootTreeReader(chainPhoton);
1306
1307 TChain *chainNeutralHadron = new TChain("Delphes");
1308 chainNeutralHadron->Add(inputFileNeutralHadron);
1309 ExRootTreeReader *treeReaderNeutralHadron = new ExRootTreeReader(chainNeutralHadron);
1310
1311 TChain *chainJet = new TChain("Delphes");
1312 chainJet->Add(inputFileJet);
1313 ExRootTreeReader *treeReaderJet = new ExRootTreeReader(chainJet);
1314
1315 TChain *chainBJet = new TChain("Delphes");
1316 chainBJet->Add(inputFileBJet);
1317 ExRootTreeReader *treeReaderBJet = new ExRootTreeReader(chainBJet);
1318
1319 TChain *chainCJet = new TChain("Delphes");
1320 chainCJet->Add(inputFileCJet);
1321 ExRootTreeReader *treeReaderCJet = new ExRootTreeReader(chainCJet);
1322
1323 TChain *chainTauJet = new TChain("Delphes");
1324 chainTauJet->Add(inputFileTauJet);
1325 ExRootTreeReader *treeReaderTauJet = new ExRootTreeReader(chainTauJet);
1326
1327 TClonesArray *branchParticleElectron = treeReaderElectron->UseBranch("Particle");
1328 TClonesArray *branchTrackElectron = treeReaderElectron->UseBranch("Track");
1329 TClonesArray *branchElectron = treeReaderElectron->UseBranch("Electron");
1330 TClonesArray *branchElectronPF = treeReaderElectron->UseBranch("ElectronPF");
1331
1332 TClonesArray *branchParticleMuon = treeReaderMuon->UseBranch("Particle");
1333 TClonesArray *branchTrackMuon = treeReaderMuon->UseBranch("Track");
1334 TClonesArray *branchMuon = treeReaderMuon->UseBranch("Muon");
1335
1336 TClonesArray *branchParticlePion = treeReaderPion->UseBranch("Particle");
1337 TClonesArray *branchTrackPion = treeReaderPion->UseBranch("Track");
1338 TClonesArray *branchPion = treeReaderPion->UseBranch("Pion");
1339
1340 TClonesArray *branchParticlePhoton = treeReaderPhoton->UseBranch("Particle");
1341 TClonesArray *branchTowerPhoton = treeReaderPhoton->UseBranch("Tower");
1342 TClonesArray *branchPhoton = treeReaderPhoton->UseBranch("Photon");
1343
1344 TClonesArray *branchParticleNeutralHadron = treeReaderNeutralHadron->UseBranch("Particle");
1345 TClonesArray *branchTowerNeutralHadron = treeReaderNeutralHadron->UseBranch("Tower");
1346
1347 TClonesArray *branchGenJet = treeReaderJet->UseBranch("GenJet");
1348 TClonesArray *branchParticleJet = treeReaderJet->UseBranch("Particle");
1349 TClonesArray *branchPFJet = treeReaderJet->UseBranch("Jet");
1350 TClonesArray *branchCaloJet = treeReaderJet->UseBranch("CaloJet");
1351
1352 TClonesArray *branchParticleBJet = treeReaderBJet->UseBranch("Particle");
1353 TClonesArray *branchPFBJet = treeReaderBJet->UseBranch("Jet");
1354
1355 TClonesArray *branchParticleCJet = treeReaderCJet->UseBranch("Particle");
1356 TClonesArray *branchPFCJet = treeReaderCJet->UseBranch("Jet");
1357
1358 TClonesArray *branchParticleTauJet = treeReaderTauJet->UseBranch("Particle");
1359 TClonesArray *branchPFTauJet = treeReaderTauJet->UseBranch("Jet");
1360
1361 TClonesArray *branchGenScalarHT = treeReaderJet->UseBranch("GenScalarHT");
1362 TClonesArray *branchMet = treeReaderJet->UseBranch("MissingET");
1363 TClonesArray *branchCaloMet = treeReaderJet->UseBranch("CaloMissingET");
1364
1365 std::vector<Color_t> colors;
1366
1367 colors.push_back(kBlack);
1368 colors.push_back(kBlue);
1369 colors.push_back(kRed);
1370 colors.push_back(kGreen+1);
1371 colors.push_back(kMagenta+1);
1372 colors.push_back(kOrange);
1373
1374 std::vector<Int_t> markerStyles;
1375
1376 markerStyles.push_back(20);
1377 markerStyles.push_back(21);
1378 markerStyles.push_back(22);
1379 markerStyles.push_back(23);
1380 markerStyles.push_back(33);
1381 markerStyles.push_back(34);
1382
1383 TString pdfOutput(outputFile);
1384 pdfOutput.ReplaceAll(".root", ".pdf");
1385
1386 TString figPath = inputFilePion;
1387 figPath.ReplaceAll(".root", "");
1388 figPath.ReplaceAll("root", "www/fig");
1389 Int_t lastSlash = figPath.Last('/');
1390 Int_t sizePath = figPath.Length();
1391 figPath.Remove(lastSlash+1,sizePath);
1392
1393 TString header = pdfOutput;
1394 header.ReplaceAll(".pdf", "");
1395 header.ReplaceAll("validation_", "");
1396 lastSlash = header.Last('/');
1397 sizePath = header.Length();
1398 header.Remove(0,lastSlash+1);
1399
1400 TString vrs(version);
1401
1402 TPaveText *pave = new TPaveText(0.0, 0.89, 0.94, 0.94,"NDC");
1403 pave->SetTextAlign(kHAlignRight);
1404 pave->SetTextFont(132);
1405 pave->SetBorderSize(0);
1406 pave->SetShadowColor(0);
1407 pave->SetFillColor(0);
1408 pave->SetFillStyle(0);
1409 pave->AddText("Delphes "+vrs+" - "+header);
1410
1411 TString s_etaMin, s_etaMax, s_eta, s_pt, s_e;
1412
1413 Double_t ptMin = 1.;
1414 Double_t ptMax = 50000.;
1415
1416 Double_t etaMin = -6.;
1417 Double_t etaMax = 6.;
1418
1419 std::vector<Double_t> ptVals;
1420
1421 ptVals.push_back(1.);
1422 ptVals.push_back(10.);
1423 ptVals.push_back(100.);
1424 ptVals.push_back(1000.);
1425 ptVals.push_back(10000.);
1426
1427 std::vector<Double_t> etaVals;
1428
1429 etaVals.push_back(0.);
1430 etaVals.push_back(1.5);
1431 etaVals.push_back(2.5);
1432 etaVals.push_back(4.0);
1433 etaVals.push_back(6.0);
1434
1435 const int n_etabins = etaVals.size()-1;
1436 const int n_ptbins = ptVals.size();
1437
1438
1439 //////////////////////////
1440 // Tracking performance //
1441 //////////////////////////
1442
1443 // --------- Pion Tracks --------- //
1444
1445 TMultiGraph *mg_trkpi_res_pt = new TMultiGraph("","");
1446 TMultiGraph *mg_trkpi_eff_pt = new TMultiGraph("","");
1447 TMultiGraph *mg_trkpi_res_eta = new TMultiGraph("","");
1448 TMultiGraph *mg_trkpi_eff_eta = new TMultiGraph("","");
1449
1450 TLegend *leg_trkpi_res_pt = new TLegend(0.55,0.22,0.90,0.48);
1451 TLegend *leg_trkpi_eff_pt = (TLegend*)leg_trkpi_res_pt->Clone();
1452 TLegend *leg_trkpi_res_eta = (TLegend*)leg_trkpi_res_pt->Clone();
1453 TLegend *leg_trkpi_eff_eta = (TLegend*)leg_trkpi_res_eta->Clone();
1454
1455
1456 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];
1457 TH1D* h_trkpi_eff_pt, *h_trkpi_eff_eta;
1458
1459 std::vector<resolPlot> plots_trkpi_res_pt[n_etabins], plots_trkpi_res_eta[n_ptbins];
1460
1461 // loop over eta bins
1462 for (k = 0; k < etaVals.size()-1; k++)
1463 {
1464 HistogramsCollection(&plots_trkpi_res_pt[k], TMath::Log10(ptMin), TMath::Log10(ptMax), "trkpi");
1465 GetPtres<Track>(&plots_trkpi_res_pt[k], branchTrackPion, branchParticlePion, 211, etaVals.at(k), etaVals.at(k+1), treeReaderPion);
1466 gr_trkpi_res_pt[k] = EresGraph(&plots_trkpi_res_pt[k]);
1467
1468 h_trkpi_eff_pt = GetEffPt<Track>(branchTrackPion, branchParticlePion, "Pion", 211, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderPion);
1469 gr_trkpi_eff_pt[k] = TGraphErrors(h_trkpi_eff_pt);
1470
1471 s_etaMin = Form("%.1f",etaVals.at(k));
1472 s_etaMax = Form("%.1f",etaVals.at(k+1));
1473
1474 s_eta = "#pi^{ #pm} , " + s_etaMin + " < | #eta | < "+s_etaMax;
1475
1476 gr_trkpi_res_pt[k].SetName("trkRes_"+s_etaMin+"_"+s_etaMax);
1477 gr_trkpi_eff_pt[k].SetName("trkEff_"+s_etaMin+"_"+s_etaMax);
1478
1479 addResoGraph(mg_trkpi_res_pt, &gr_trkpi_res_pt[k], leg_trkpi_res_pt, markerStyles.at(k), colors.at(k), s_eta);
1480 addResoGraph(mg_trkpi_eff_pt, &gr_trkpi_eff_pt[k], leg_trkpi_eff_pt, markerStyles.at(k), colors.at(k), s_eta);
1481 }
1482
1483 // loop over pt
1484 for (k = 0; k < ptVals.size(); k++)
1485 {
1486 HistogramsCollectionVsEta(&plots_trkpi_res_eta[k], etaMin, etaMax, 0.5*ptVals.at(k), 2.0*ptVals.at(k), "trkpi", 0.0, 2.0);
1487 GetPtresVsEta<Track>(&plots_trkpi_res_eta[k], branchTrackPion, branchParticlePion, 211, 0.5*ptVals.at(k), 2.0*ptVals.at(k), treeReaderPion);
1488 gr_trkpi_res_eta[k] = EresGraphVsEta(&plots_trkpi_res_eta[k]);
1489
1490 h_trkpi_eff_eta = GetEffEta<Track>(branchTrackPion, branchParticlePion, "Pion", 211, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderPion);
1491 gr_trkpi_eff_eta[k] = TGraphErrors(h_trkpi_eff_eta);
1492
1493 s_pt = Form("#pi^{ #pm} , p_{T} = %.0f GeV",ptVals.at(k));
1494 if(ptVals.at(k) >= 1000.) s_pt = Form("#pi^{ #pm} , p_{T} = %.0f TeV",ptVals.at(k)/1000.);
1495
1496 addResoGraph(mg_trkpi_res_eta, &gr_trkpi_res_eta[k], leg_trkpi_res_eta, markerStyles.at(k), colors.at(k), s_pt );
1497 addResoGraph(mg_trkpi_eff_eta, &gr_trkpi_eff_eta[k], leg_trkpi_eff_eta, markerStyles.at(k), colors.at(k), s_pt );
1498
1499 }
1500
1501 TCanvas *c_trkpi_res_pt = new TCanvas("","", 800, 600);
1502
1503 mg_trkpi_res_pt->Draw("APE");
1504 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);
1505 leg_trkpi_res_pt->Draw();
1506 pave->Draw();
1507
1508 c_trkpi_res_pt->Print(pdfOutput+"(","pdf");
1509 c_trkpi_res_pt->Print(figPath+"img_trkpi_res_pt.pdf","pdf");
1510 c_trkpi_res_pt->Print(figPath+"img_trkpi_res_pt.png","png");
1511
1512 TCanvas *c_trkpi_res_eta = new TCanvas("","", 800, 600);
1513
1514 mg_trkpi_res_eta->Draw("APE");
1515 DrawAxis(mg_trkpi_res_eta, leg_trkpi_res_eta, etaMin, etaMax, 0.01, 100, " #eta ", "(track resolution in p_{T})/p_{T} (%)", false, true);
1516 leg_trkpi_res_eta->Draw();
1517 pave->Draw();
1518
1519 c_trkpi_res_eta->Print(pdfOutput,"pdf");
1520 c_trkpi_res_eta->Print(figPath+"img_trkpi_res_eta.pdf","pdf");
1521 c_trkpi_res_eta->Print(figPath+"img_trkpi_res_eta.png","png");
1522
1523 TCanvas *c_trkpi_eff_pt = new TCanvas("","", 800, 600);
1524
1525 mg_trkpi_eff_pt->Draw("APE");
1526 DrawAxis(mg_trkpi_eff_pt, leg_trkpi_eff_pt, ptMin, ptMax, 0.0, 100, "p_{T} [GeV]", "tracking efficiency (%)", true, false);
1527 leg_trkpi_eff_pt->Draw();
1528 pave->Draw();
1529
1530 c_trkpi_eff_pt->Print(pdfOutput,"pdf");
1531 c_trkpi_eff_pt->Print(figPath+"img_trkpi_eff_pt.pdf","pdf");
1532 c_trkpi_eff_pt->Print(figPath+"img_trkpi_eff_pt.png","png");
1533
1534 TCanvas *c_trkpi_eff_eta = new TCanvas("","", 800, 600);
1535
1536 mg_trkpi_eff_eta->Draw("APE");
1537 DrawAxis(mg_trkpi_eff_eta, leg_trkpi_eff_eta, etaMin, etaMax, 0.0, 100, " #eta ", "tracking efficiency (%)", false, false);
1538 leg_trkpi_eff_eta->Draw();
1539 pave->Draw();
1540
1541
1542 c_trkpi_eff_eta->Print(pdfOutput,"pdf");
1543 c_trkpi_eff_eta->Print(figPath+"img_trkpi_eff_eta.pdf","pdf");
1544 c_trkpi_eff_eta->Print(figPath+"img_trkpi_eff_eta.png","png");
1545
1546 // --------- Electron Tracks --------- //
1547
1548 TMultiGraph *mg_trkele_res_pt = new TMultiGraph("","");
1549 TMultiGraph *mg_trkele_eff_pt = new TMultiGraph("","");
1550 TMultiGraph *mg_trkele_res_eta = new TMultiGraph("","");
1551 TMultiGraph *mg_trkele_eff_eta = new TMultiGraph("","");
1552
1553 TLegend *leg_trkele_res_pt = new TLegend(0.55,0.22,0.90,0.48);
1554 TLegend *leg_trkele_eff_pt = (TLegend*)leg_trkele_res_pt->Clone();
1555 TLegend *leg_trkele_res_eta = (TLegend*)leg_trkele_res_pt->Clone();
1556 TLegend *leg_trkele_eff_eta = (TLegend*)leg_trkele_res_eta->Clone();
1557
1558 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];
1559 TH1D* h_trkele_eff_pt, *h_trkele_eff_eta;
1560
1561 std::vector<resolPlot> plots_trkele_res_pt[n_etabins], plots_trkele_res_eta[n_ptbins];
1562
1563 // loop over eta bins
1564 for (k = 0; k < etaVals.size()-1; k++)
1565 {
1566 HistogramsCollection(&plots_trkele_res_pt[k], TMath::Log10(ptMin), TMath::Log10(ptMax), "trkele");
1567 GetPtres<Track>(&plots_trkele_res_pt[k], branchTrackElectron, branchParticleElectron, 11, etaVals.at(k), etaVals.at(k+1), treeReaderElectron);
1568 gr_trkele_res_pt[k] = EresGraph(&plots_trkele_res_pt[k]);
1569
1570 h_trkele_eff_pt = GetEffPt<Track>(branchTrackElectron, branchParticleElectron, "Electron", 11, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderElectron);
1571 gr_trkele_eff_pt[k] = TGraphErrors(h_trkele_eff_pt);
1572
1573 s_etaMin = Form("%.1f",etaVals.at(k));
1574 s_etaMax = Form("%.1f",etaVals.at(k+1));
1575
1576 s_eta = "e^{ #pm} , " + s_etaMin + " < | #eta | < "+s_etaMax;
1577
1578 gr_trkele_res_pt[k].SetName("trkRes_"+s_etaMin+"_"+s_etaMax);
1579 gr_trkele_eff_pt[k].SetName("trkEff_"+s_etaMin+"_"+s_etaMax);
1580
1581 addResoGraph(mg_trkele_res_pt, &gr_trkele_res_pt[k], leg_trkele_res_pt, markerStyles.at(k), colors.at(k), s_eta);
1582 addResoGraph(mg_trkele_eff_pt, &gr_trkele_eff_pt[k], leg_trkele_eff_pt, markerStyles.at(k), colors.at(k), s_eta);
1583 }
1584
1585 // loop over pt
1586 for (k = 0; k < ptVals.size(); k++)
1587 {
1588 HistogramsCollectionVsEta(&plots_trkele_res_eta[k], etaMin, etaMax, 0.5*ptVals.at(k), 2.0*ptVals.at(k), "trkele", 0.0, 2.0);
1589 GetPtresVsEta<Track>(&plots_trkele_res_eta[k], branchTrackElectron, branchParticleElectron, 11, 0.5*ptVals.at(k), 2.0*ptVals.at(k), treeReaderElectron);
1590 gr_trkele_res_eta[k] = EresGraphVsEta(&plots_trkele_res_eta[k]);
1591
1592 h_trkele_eff_eta = GetEffEta<Track>(branchTrackElectron, branchParticleElectron, "Electron", 11, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderElectron);
1593 gr_trkele_eff_eta[k] = TGraphErrors(h_trkele_eff_eta);
1594
1595 s_pt = Form("e^{ #pm} , p_{T} = %.0f GeV",ptVals.at(k));
1596 if(ptVals.at(k) >= 1000.) s_pt = Form("e^{ #pm} , p_{T} = %.0f TeV",ptVals.at(k)/1000.);
1597
1598 addResoGraph(mg_trkele_res_eta, &gr_trkele_res_eta[k], leg_trkele_res_eta, markerStyles.at(k), colors.at(k), s_pt );
1599 addResoGraph(mg_trkele_eff_eta, &gr_trkele_eff_eta[k], leg_trkele_eff_eta, markerStyles.at(k), colors.at(k), s_pt );
1600
1601 }
1602
1603 TCanvas *c_trkele_res_pt = new TCanvas("","", 800, 600);
1604
1605 mg_trkele_res_pt->Draw("APE");
1606 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);
1607 leg_trkele_res_pt->Draw();
1608 pave->Draw();
1609
1610 c_trkele_res_pt->Print(pdfOutput,"pdf");
1611 c_trkele_res_pt->Print(figPath+"img_trkele_res_pt.pdf","pdf");
1612 c_trkele_res_pt->Print(figPath+"img_trkele_res_pt.png","png");
1613
1614 TCanvas *c_trkele_res_eta = new TCanvas("","", 800, 600);
1615
1616 mg_trkele_res_eta->Draw("APE");
1617 DrawAxis(mg_trkele_res_eta, leg_trkele_res_eta, etaMin, etaMax, 0.01, 100, " #eta ", "(track resolution in p_{T})/p_{T} (%)", false, true);
1618 leg_trkele_res_eta->Draw();
1619 pave->Draw();
1620
1621 c_trkele_res_eta->Print(pdfOutput,"pdf");
1622 c_trkele_res_eta->Print(figPath+"img_trkele_res_eta.pdf","pdf");
1623 c_trkele_res_eta->Print(figPath+"img_trkele_res_eta.png","png");
1624
1625 TCanvas *c_trkele_eff_pt = new TCanvas("","", 800, 600);
1626
1627 mg_trkele_eff_pt->Draw("APE");
1628 DrawAxis(mg_trkele_eff_pt, leg_trkele_eff_pt, ptMin, ptMax, 0.0, 100, "p_{T} [GeV]", "tracking efficiency (%)", true, false);
1629 leg_trkele_eff_pt->Draw();
1630 pave->Draw();
1631
1632 c_trkele_eff_pt->Print(pdfOutput,"pdf");
1633 c_trkele_eff_pt->Print(figPath+"img_trkele_eff_pt.pdf","pdf");
1634 c_trkele_eff_pt->Print(figPath+"img_trkele_eff_pt.png","png");
1635
1636 TCanvas *c_trkele_eff_eta = new TCanvas("","", 800, 600);
1637
1638 mg_trkele_eff_eta->Draw("APE");
1639 DrawAxis(mg_trkele_eff_eta, leg_trkele_eff_eta, etaMin, etaMax, 0.0, 100, " #eta ", "tracking efficiency (%)", false, false);
1640 leg_trkele_eff_eta->Draw();
1641 pave->Draw();
1642
1643 c_trkele_eff_eta->Print(pdfOutput,"pdf");
1644 c_trkele_eff_eta->Print(figPath+"img_trkele_eff_eta.pdf","pdf");
1645 c_trkele_eff_eta->Print(figPath+"img_trkele_eff_eta.png","png");
1646
1647
1648 // --------- Muon Tracks --------- //
1649
1650 TMultiGraph *mg_trkmu_res_pt = new TMultiGraph("","");
1651 TMultiGraph *mg_trkmu_eff_pt = new TMultiGraph("","");
1652 TMultiGraph *mg_trkmu_res_eta = new TMultiGraph("","");
1653 TMultiGraph *mg_trkmu_eff_eta = new TMultiGraph("","");
1654
1655 TLegend *leg_trkmu_res_pt = new TLegend(0.55,0.22,0.90,0.48);
1656 TLegend *leg_trkmu_eff_pt = (TLegend*)leg_trkmu_res_pt->Clone();
1657
1658 TLegend *leg_trkmu_res_eta = (TLegend*)leg_trkmu_res_pt->Clone();
1659 TLegend *leg_trkmu_eff_eta = (TLegend*)leg_trkmu_res_eta->Clone();
1660
1661
1662 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];
1663 TH1D* h_trkmu_eff_pt, *h_trkmu_eff_eta;
1664
1665 std::vector<resolPlot> plots_trkmu_res_pt[n_etabins], plots_trkmu_res_eta[n_ptbins];
1666
1667 // loop over eta bins
1668 for (k = 0; k < etaVals.size()-1; k++)
1669 {
1670 HistogramsCollection(&plots_trkmu_res_pt[k], TMath::Log10(ptMin), TMath::Log10(ptMax), "trkmu");
1671 GetPtres<Track>(&plots_trkmu_res_pt[k], branchTrackMuon, branchParticleMuon, 13, etaVals.at(k), etaVals.at(k+1), treeReaderMuon);
1672 gr_trkmu_res_pt[k] = EresGraph(&plots_trkmu_res_pt[k]);
1673
1674 h_trkmu_eff_pt = GetEffPt<Track>(branchTrackMuon, branchParticleMuon, "Muon", 13, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderMuon);
1675 gr_trkmu_eff_pt[k] = TGraphErrors(h_trkmu_eff_pt);
1676
1677 s_etaMin = Form("%.1f",etaVals.at(k));
1678 s_etaMax = Form("%.1f",etaVals.at(k+1));
1679
1680 s_eta = "#mu^{ #pm} , " + s_etaMin + " < | #eta | < "+s_etaMax;
1681
1682 gr_trkmu_res_pt[k].SetName("trkRes_"+s_etaMin+"_"+s_etaMax);
1683 gr_trkmu_eff_pt[k].SetName("trkEff_"+s_etaMin+"_"+s_etaMax);
1684
1685 addResoGraph(mg_trkmu_res_pt, &gr_trkmu_res_pt[k], leg_trkmu_res_pt, markerStyles.at(k), colors.at(k), s_eta);
1686 addResoGraph(mg_trkmu_eff_pt, &gr_trkmu_eff_pt[k], leg_trkmu_eff_pt, markerStyles.at(k), colors.at(k), s_eta);
1687 }
1688
1689 // loop over pt
1690 for (k = 0; k < ptVals.size(); k++)
1691 {
1692 HistogramsCollectionVsEta(&plots_trkmu_res_eta[k], etaMin, etaMax, 0.5*ptVals.at(k), 2.0*ptVals.at(k), "trkmu", 0.0, 2.0);
1693 GetPtresVsEta<Track>(&plots_trkmu_res_eta[k], branchTrackMuon, branchParticleMuon, 13, 0.5*ptVals.at(k), 2.0*ptVals.at(k), treeReaderMuon);
1694 gr_trkmu_res_eta[k] = EresGraphVsEta(&plots_trkmu_res_eta[k]);
1695
1696 h_trkmu_eff_eta = GetEffEta<Track>(branchTrackMuon, branchParticleMuon, "Muon", 13, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderMuon);
1697 gr_trkmu_eff_eta[k] = TGraphErrors(h_trkmu_eff_eta);
1698
1699 s_pt = Form("#mu^{ #pm} , p_{T} = %.0f GeV",ptVals.at(k));
1700 if(ptVals.at(k) >= 1000.) s_pt = Form("#mu^{ #pm} , p_{T} = %.0f TeV",ptVals.at(k)/1000.);
1701
1702 addResoGraph(mg_trkmu_res_eta, &gr_trkmu_res_eta[k], leg_trkmu_res_eta, markerStyles.at(k), colors.at(k), s_pt );
1703 addResoGraph(mg_trkmu_eff_eta, &gr_trkmu_eff_eta[k], leg_trkmu_eff_eta, markerStyles.at(k), colors.at(k), s_pt );
1704
1705 }
1706
1707 TCanvas *c_trkmu_res_pt = new TCanvas("","", 800, 600);
1708
1709 mg_trkmu_res_pt->Draw("APE");
1710 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);
1711 leg_trkmu_res_pt->Draw();
1712 pave->Draw();
1713
1714 c_trkmu_res_pt->Print(pdfOutput,"pdf");
1715 c_trkmu_res_pt->Print(figPath+"img_trkmu_res_pt.pdf","pdf");
1716 c_trkmu_res_pt->Print(figPath+"img_trkmu_res_pt.png","png");
1717
1718 TCanvas *c_trkmu_res_eta = new TCanvas("","", 800, 600);
1719
1720 mg_trkmu_res_eta->Draw("APE");
1721 DrawAxis(mg_trkmu_res_eta, leg_trkmu_res_eta, etaMin, etaMax, 0.01, 100, " #eta ", "(track resolution in p_{T})/p_{T} (%)", false, true);
1722 leg_trkmu_res_eta->Draw();
1723 pave->Draw();
1724
1725 c_trkmu_res_eta->Print(pdfOutput,"pdf");
1726 c_trkmu_res_eta->Print(figPath+"img_trkmu_res_eta.pdf","pdf");
1727 c_trkmu_res_eta->Print(figPath+"img_trkmu_res_eta.png","png");
1728
1729 TCanvas *c_trkmu_eff_pt = new TCanvas("","", 800, 600);
1730
1731 mg_trkmu_eff_pt->Draw("APE");
1732 DrawAxis(mg_trkmu_eff_pt, leg_trkmu_eff_pt, ptMin, ptMax, 0.0, 100, "p_{T} [GeV]", "tracking efficiency (%)", true, false);
1733 leg_trkmu_eff_pt->Draw();
1734 pave->Draw();
1735
1736 c_trkmu_eff_pt->Print(pdfOutput,"pdf");
1737 c_trkmu_eff_pt->Print(figPath+"img_trkmu_eff_pt.pdf","pdf");
1738 c_trkmu_eff_pt->Print(figPath+"img_trkmu_eff_pt.png","png");
1739
1740 TCanvas *c_trkmu_eff_eta = new TCanvas("","", 800, 600);
1741
1742 mg_trkmu_eff_eta->Draw("APE");
1743 DrawAxis(mg_trkmu_eff_eta, leg_trkmu_eff_eta, etaMin, etaMax, 0.0, 100, " #eta ", "tracking efficiency (%)", false, false);
1744 leg_trkmu_eff_eta->Draw();
1745 pave->Draw();
1746
1747 c_trkmu_eff_eta->Print(pdfOutput,"pdf");
1748 c_trkmu_eff_eta->Print(figPath+"img_trkmu_eff_eta.pdf","pdf");
1749 c_trkmu_eff_eta->Print(figPath+"img_trkmu_eff_eta.png","png");
1750
1751
1752 //////////////////////
1753 // Ecal performance //
1754 //////////////////////
1755
1756
1757 TMultiGraph *mg_ecal_res_e = new TMultiGraph("","");
1758 TMultiGraph *mg_ecal_res_eta = new TMultiGraph("","");
1759
1760 TLegend *leg_ecal_res_e = new TLegend(0.55,0.64,0.90,0.90);
1761 TLegend *leg_ecal_res_eta = new TLegend(0.60,0.59,0.95,0.90);
1762
1763 TGraphErrors gr_ecal_res_e[n_etabins], gr_ecal_res_eta[n_ptbins];
1764
1765 std::vector<resolPlot> plots_ecal_res_e[n_etabins], plots_ecal_res_eta[n_ptbins];
1766
1767 // loop over eta bins
1768 for (k = 0; k < etaVals.size()-1; k++)
1769 {
1770 HistogramsCollection(&plots_ecal_res_e[k], TMath::Log10(ptMin), TMath::Log10(ptMax), "ecal");
1771 GetEres<Tower>(&plots_ecal_res_e[k], branchTowerPhoton, branchParticlePhoton, 22, etaVals.at(k), etaVals.at(k+1), treeReaderPhoton);
1772 gr_ecal_res_e[k] = EresGraph(&plots_ecal_res_e[k]);
1773
1774 s_etaMin = Form("%.1f",etaVals.at(k));
1775 s_etaMax = Form("%.1f",etaVals.at(k+1));
1776
1777 s_eta = "#gamma , " + s_etaMin + " < | #eta | < "+s_etaMax;
1778
1779 gr_ecal_res_e[k].SetName("trkRes_"+s_etaMin+"_"+s_etaMax);
1780
1781 addResoGraph(mg_ecal_res_e, &gr_ecal_res_e[k], leg_ecal_res_e, markerStyles.at(k), colors.at(k), s_eta);
1782 }
1783
1784 // loop over pt
1785 for (k = 0; k < ptVals.size(); k++)
1786 {
1787 HistogramsCollectionVsEta(&plots_ecal_res_eta[k], etaMin, etaMax, 0.5*ptVals.at(k), 2.0*ptVals.at(k), "ecal", 0.0, 2.0);
1788 GetEresVsEta<Tower>(&plots_ecal_res_eta[k], branchTowerPhoton, branchParticlePhoton, 22, 0.5*ptVals.at(k), 2.0*ptVals.at(k), treeReaderPhoton);
1789 gr_ecal_res_eta[k] = EresGraphVsEta(&plots_ecal_res_eta[k]);
1790
1791 s_e = Form("#gamma , E = %.0f GeV",ptVals.at(k));
1792 if(ptVals.at(k) >= 1000.) s_e = Form("#gamma , E = %.0f TeV",ptVals.at(k)/1000.);
1793
1794 addResoGraph(mg_ecal_res_eta, &gr_ecal_res_eta[k], leg_ecal_res_eta, markerStyles.at(k), colors.at(k), s_e );
1795 }
1796
1797 TCanvas *c_ecal_res_e = new TCanvas("","", 800, 600);
1798
1799 mg_ecal_res_e->Draw("APE");
1800 // DrawAxis(mg_ecal_res_e, leg_ecal_res_e, ptMin, ptMax, 0.5, 100, "E [GeV]", "(ECAL resolution in E)/E (%)", true, true);
1801 DrawAxis(mg_ecal_res_e, leg_ecal_res_e, ptMin, ptMax, 0.0, 20, "E [GeV]", "(ECAL resolution in E)/E (%)", true, false);
1802 leg_ecal_res_e->Draw();
1803 pave->Draw();
1804
1805 c_ecal_res_e->Print(pdfOutput,"pdf");
1806 c_ecal_res_e->Print(figPath+"img_ecal_res_e.pdf","pdf");
1807 c_ecal_res_e->Print(figPath+"img_ecal_res_e.png","png");
1808
1809 TCanvas *c_ecal_res_eta = new TCanvas("","", 800, 600);
1810
1811 mg_ecal_res_eta->Draw("APE");
1812 //DrawAxis(mg_ecal_res_eta, leg_ecal_res_eta, etaMin, etaMax, 0.5, 100, " #eta ", "(ECAL resolution in E)/E (%)", false, true);
1813 DrawAxis(mg_ecal_res_eta, leg_ecal_res_eta, etaMin, etaMax, 0.0, 20, " #eta ", "(ECAL resolution in E)/E (%)", false, false);
1814 leg_ecal_res_eta->Draw();
1815 pave->Draw();
1816
1817 c_ecal_res_eta->Print(pdfOutput,"pdf");
1818 c_ecal_res_eta->Print(figPath+"img_ecal_res_eta.pdf","pdf");
1819 c_ecal_res_eta->Print(figPath+"img_ecal_res_eta.png","png");
1820
1821 //////////////////////
1822 // Hcal performance //
1823 //////////////////////
1824
1825
1826 TMultiGraph *mg_hcal_res_e = new TMultiGraph("","");
1827 TMultiGraph *mg_hcal_res_eta = new TMultiGraph("","");
1828
1829 TLegend *leg_hcal_res_e = new TLegend(0.55,0.64,0.90,0.90);
1830 TLegend *leg_hcal_res_eta = new TLegend(0.60,0.59,0.95,0.90);
1831
1832 TGraphErrors gr_hcal_res_e[n_etabins], gr_hcal_res_eta[n_ptbins];
1833
1834 std::vector<resolPlot> plots_hcal_res_e[n_etabins], plots_hcal_res_eta[n_ptbins];
1835
1836 // loop over eta bins
1837 for (k = 0; k < etaVals.size()-1; k++)
1838 {
1839 HistogramsCollection(&plots_hcal_res_e[k], TMath::Log10(ptMin), TMath::Log10(ptMax), "hcal");
1840 GetEres<Tower>(&plots_hcal_res_e[k], branchTowerNeutralHadron, branchParticleNeutralHadron, 2112, etaVals.at(k), etaVals.at(k+1), treeReaderNeutralHadron);
1841
1842 gr_hcal_res_e[k] = EresGraph(&plots_hcal_res_e[k]);
1843
1844 s_etaMin = Form("%.1f",etaVals.at(k));
1845 s_etaMax = Form("%.1f",etaVals.at(k+1));
1846
1847 s_eta = "n , " + s_etaMin + " < | #eta | < "+s_etaMax;
1848
1849 gr_hcal_res_e[k].SetName("trkRes_"+s_etaMin+"_"+s_etaMax);
1850
1851 addResoGraph(mg_hcal_res_e, &gr_hcal_res_e[k], leg_hcal_res_e, markerStyles.at(k), colors.at(k), s_eta);
1852 }
1853
1854 // loop over pt
1855 for (k = 0; k < ptVals.size(); k++)
1856 {
1857 HistogramsCollectionVsEta(&plots_hcal_res_eta[k], etaMin, etaMax, 0.5*ptVals.at(k), 2.0*ptVals.at(k), "hcal", 0.0, 2.0);
1858 GetEresVsEta<Tower>(&plots_hcal_res_eta[k], branchTowerNeutralHadron, branchParticleNeutralHadron, 2112, 0.5*ptVals.at(k), 2.0*ptVals.at(k), treeReaderNeutralHadron);
1859 gr_hcal_res_eta[k] = EresGraphVsEta(&plots_hcal_res_eta[k]);
1860
1861 s_e = Form("n , E = %.0f GeV",ptVals.at(k));
1862 if(ptVals.at(k) >= 1000.) s_e = Form("n , E = %.0f TeV",ptVals.at(k)/1000.);
1863
1864 addResoGraph(mg_hcal_res_eta, &gr_hcal_res_eta[k], leg_hcal_res_eta, markerStyles.at(k), colors.at(k), s_e );
1865 }
1866
1867
1868 TCanvas *c_hcal_res_e = new TCanvas("","", 800, 600);
1869
1870 mg_hcal_res_e->Draw("APE");
1871 //DrawAxis(mg_hcal_res_e, leg_hcal_res_e, ptMin, ptMax, 1, 100, "E [GeV]", "(HCAL resolution in E)/E (%)", true, true);
1872 DrawAxis(mg_hcal_res_e, leg_hcal_res_e, ptMin, ptMax, 0.0, 50, "E [GeV]", "(HCAL resolution in E)/E (%)", true, false);
1873 leg_hcal_res_e->Draw();
1874 pave->Draw();
1875
1876 c_hcal_res_e->Print(pdfOutput,"pdf");
1877 c_hcal_res_e->Print(figPath+"img_hcal_res_e.pdf","pdf");
1878 c_hcal_res_e->Print(figPath+"img_hcal_res_e.png","png");
1879
1880 TCanvas *c_hcal_res_eta = new TCanvas("","", 800, 600);
1881
1882 mg_hcal_res_eta->Draw("APE");
1883 //DrawAxis(mg_hcal_res_eta, leg_hcal_res_eta, etaMin, etaMax, 1, 100, " #eta ", "(HCAL resolution in E)/E (%)", false, true);
1884 DrawAxis(mg_hcal_res_eta, leg_hcal_res_eta, etaMin, etaMax, 0.0, 50, " #eta ", "(HCAL resolution in E)/E (%)", false, false);
1885 leg_hcal_res_eta->Draw();
1886 pave->Draw();
1887
1888 c_hcal_res_eta->Print(pdfOutput,"pdf");
1889 c_hcal_res_eta->Print(figPath+"img_hcal_res_eta.pdf","pdf");
1890 c_hcal_res_eta->Print(figPath+"img_hcal_res_eta.png","png");
1891
1892 ////////////////////
1893 // PF - electrons //
1894 ////////////////////
1895
1896 TMultiGraph *mg_pfele_res_e[n_etabins];
1897 TMultiGraph *mg_pfele_res_eta[n_ptbins];
1898
1899 TLegend *leg_pfele_res_e[n_etabins];
1900 TLegend *leg_pfele_res_eta[n_ptbins];
1901
1902 TGraphErrors gr_pfele_res_e[n_etabins];
1903 TGraphErrors gr_pfele_res_eta[n_ptbins];
1904
1905 std::vector<resolPlot> plots_pfele_res_e[n_etabins], plots_pfele_res_eta[n_ptbins];
1906
1907 TCanvas *c_pfele_res_e[n_etabins];
1908 TCanvas *c_pfele_res_eta[n_ptbins];
1909
1910
1911 // loop over eta bins
1912 for (k = 0; k < etaVals.size()-1; k++)
1913 {
1914 mg_pfele_res_e[k] = new TMultiGraph("","");
1915 leg_pfele_res_e[k] = new TLegend(0.40,0.60,0.75,0.90);
1916
1917 HistogramsCollection(&plots_pfele_res_e[k], TMath::Log10(ptMin), TMath::Log10(ptMax), "pfele");
1918 GetEres<Electron>(&plots_pfele_res_e[k], branchElectronPF, branchParticleElectron, 11, etaVals.at(k), etaVals.at(k+1), treeReaderElectron);
1919 gr_pfele_res_e[k] = EresGraph(&plots_pfele_res_e[k]);
1920
1921 s_etaMin = Form("%.1f",etaVals.at(k));
1922 s_etaMax = Form("%.1f",etaVals.at(k+1));
1923 s_eta = "e^{ #pm}, "+ s_etaMin + " < | #eta | < " + s_etaMax;
1924
1925 leg_pfele_res_e[k]->SetTextFont(132);
1926 leg_pfele_res_e[k]->SetHeader(s_eta);
1927
1928 addResoGraph(mg_pfele_res_e[k], &gr_ecal_res_e[k], leg_pfele_res_e[k], markerStyles.at(0), colors.at(0), "ECAL");
1929 addResoGraph(mg_pfele_res_e[k], &gr_trkele_res_pt[k], leg_pfele_res_e[k], markerStyles.at(1), colors.at(1), "Track");
1930 addResoGraph(mg_pfele_res_e[k], &gr_pfele_res_e[k], leg_pfele_res_e[k], markerStyles.at(2), colors.at(2), "Particle-flow");
1931
1932 c_pfele_res_e[k] = new TCanvas("","", 800, 600);
1933
1934 mg_pfele_res_e[k]->Draw("APE");
1935 //DrawAxis(mg_pfele_res_e[k], leg_pfele_res_e[k], ptMin, ptMax, 0.1, 100, "E [GeV]", "(resolution in E)/E (%)", true, true);
1936 DrawAxis(mg_pfele_res_e[k], leg_pfele_res_e[k], ptMin, ptMax, 0.0, 20, "E [GeV]", "(resolution in E)/E (%)", true, false);
1937 leg_pfele_res_e[k]->Draw();
1938 pave->Draw();
1939
1940 TString s_etarange = "eta_"+s_etaMin+"_"+s_etaMax+"_";
1941
1942 c_pfele_res_e[k]->Print(pdfOutput,"pdf");
1943 c_pfele_res_e[k]->Print(figPath+"img_pfele_res_"+s_etarange+"e.pdf","pdf");
1944 c_pfele_res_e[k]->Print(figPath+"img_pfele_res_"+s_etarange+"e.png","png");
1945
1946 }
1947
1948
1949 // loop over eta bins
1950 for (k = 0; k < ptVals.size(); k++)
1951 {
1952
1953 mg_pfele_res_eta[k] = new TMultiGraph("","");
1954 leg_pfele_res_eta[k] = new TLegend(0.40,0.60,0.75,0.90);
1955
1956 HistogramsCollectionVsEta(&plots_pfele_res_eta[k], etaMin, etaMax, 0.5*ptVals.at(k), 2.0*ptVals.at(k), "pfele", 0.0, 2.0);
1957 GetEresVsEta<Electron>(&plots_pfele_res_eta[k], branchElectronPF, branchParticleElectron, 11, 0.5*ptVals.at(k), 2.0*ptVals.at(k), treeReaderElectron);
1958 gr_pfele_res_eta[k] = EresGraphVsEta(&plots_pfele_res_eta[k]);
1959
1960 s_e = Form("e^{ #pm}, E = %.0f GeV",ptVals.at(k));
1961 if(ptVals.at(k) >= 1000.) s_e = Form("e^{ #pm}, E = %.0f TeV",ptVals.at(k)/1000.);
1962
1963
1964 leg_pfele_res_eta[k]->SetTextFont(132);
1965 leg_pfele_res_eta[k]->SetHeader(s_e);
1966
1967 addResoGraph(mg_pfele_res_eta[k], &gr_ecal_res_eta[k], leg_pfele_res_eta[k], markerStyles.at(0), colors.at(0), "ECAL");
1968 addResoGraph(mg_pfele_res_eta[k], &gr_trkele_res_eta[k], leg_pfele_res_eta[k], markerStyles.at(1), colors.at(1), "Track");
1969 addResoGraph(mg_pfele_res_eta[k], &gr_pfele_res_eta[k], leg_pfele_res_eta[k], markerStyles.at(2), colors.at(2), "Particle-flow");
1970
1971 c_pfele_res_eta[k] = new TCanvas("","", 800, 600);
1972
1973 mg_pfele_res_eta[k]->Draw("APE");
1974 //DrawAxis(mg_pfele_res_eta[k], leg_pfele_res_eta[k], etaMin, etaMax, 0.1, 1000, "#eta", "(resolution in E)/E (%)", false, true);
1975 DrawAxis(mg_pfele_res_eta[k], leg_pfele_res_eta[k], etaMin, etaMax, 0.0, 50, "#eta", "(resolution in E)/E (%)", false, false);
1976 leg_pfele_res_eta[k]->Draw();
1977 pave->Draw();
1978
1979 TString s_ptrange = Form("pt_%.0f_",ptVals.at(k));
1980
1981 c_pfele_res_eta[k]->Print(pdfOutput,"pdf");
1982 c_pfele_res_eta[k]->Print(figPath+"img_pfele_res_"+s_ptrange+"eta.pdf","pdf");
1983 c_pfele_res_eta[k]->Print(figPath+"img_pfele_res_"+s_ptrange+"eta.png","png");
1984
1985 }
1986
1987 /////////////////
1988 // PF - Pions //
1989 /////////////////
1990
1991 TMultiGraph *mg_pfpi_res_e[n_etabins];
1992 TMultiGraph *mg_pfpi_res_eta[n_ptbins];
1993
1994 TLegend *leg_pfpi_res_e[n_etabins];
1995 TLegend *leg_pfpi_res_eta[n_ptbins];
1996
1997 TGraphErrors gr_pfpi_res_e[n_etabins];
1998 TGraphErrors gr_pfpi_res_eta[n_ptbins];
1999
2000 std::vector<resolPlot> plots_pfpi_res_e[n_etabins], plots_pfpi_res_eta[n_ptbins];
2001
2002 TCanvas *c_pfpi_res_e[n_etabins];
2003 TCanvas *c_pfpi_res_eta[n_ptbins];
2004
2005
2006 // loop over eta bins
2007 for (k = 0; k < etaVals.size()-1; k++)
2008 {
2009 mg_pfpi_res_e[k] = new TMultiGraph("","");
2010 leg_pfpi_res_e[k] = new TLegend(0.40,0.60,0.75,0.90);
2011
2012 HistogramsCollection(&plots_pfpi_res_e[k], TMath::Log10(ptMin), TMath::Log10(ptMax), "pfpi");
2013 GetEres<Track>(&plots_pfpi_res_e[k], branchPion, branchParticlePion, 211, etaVals.at(k), etaVals.at(k+1), treeReaderPion);
2014 gr_pfpi_res_e[k] = EresGraph(&plots_pfpi_res_e[k]);
2015
2016 s_etaMin = Form("%.1f",etaVals.at(k));
2017 s_etaMax = Form("%.1f",etaVals.at(k+1));
2018 s_eta = "#pi^{ #pm}, "+ s_etaMin + " < | #eta | < " + s_etaMax;
2019
2020 leg_pfpi_res_e[k]->SetTextFont(132);
2021 leg_pfpi_res_e[k]->SetHeader(s_eta);
2022
2023 addResoGraph(mg_pfpi_res_e[k], &gr_hcal_res_e[k], leg_pfpi_res_e[k], markerStyles.at(0), colors.at(0), "HCAL");
2024 addResoGraph(mg_pfpi_res_e[k], &gr_trkpi_res_pt[k], leg_pfpi_res_e[k], markerStyles.at(1), colors.at(1), "Track");
2025 addResoGraph(mg_pfpi_res_e[k], &gr_pfpi_res_e[k], leg_pfpi_res_e[k], markerStyles.at(2), colors.at(2), "Particle-flow");
2026
2027 c_pfpi_res_e[k] = new TCanvas("","", 800, 600);
2028
2029 mg_pfpi_res_e[k]->Draw("APE");
2030 //DrawAxis(mg_pfpi_res_e[k], leg_pfpi_res_e[k], ptMin, ptMax, 0.1, 100, "E [GeV]", "(resolution in E)/E (%)", true, true);
2031 DrawAxis(mg_pfpi_res_e[k], leg_pfpi_res_e[k], ptMin, ptMax, 0.1, 50, "E [GeV]", "(resolution in E)/E (%)", true, false);
2032 leg_pfpi_res_e[k]->Draw();
2033 pave->Draw();
2034
2035 TString s_etarange = "eta_"+s_etaMin+"_"+s_etaMax+"_";
2036
2037 c_pfpi_res_e[k]->Print(pdfOutput,"pdf");
2038 c_pfpi_res_e[k]->Print(figPath+"img_pfpi_res_"+s_etarange+"e.pdf","pdf");
2039 c_pfpi_res_e[k]->Print(figPath+"img_pfpi_res_"+s_etarange+"e.png","png");
2040
2041 }
2042
2043
2044 // loop over eta bins
2045 for (k = 0; k < ptVals.size(); k++)
2046 {
2047
2048 mg_pfpi_res_eta[k] = new TMultiGraph("","");
2049 leg_pfpi_res_eta[k] = new TLegend(0.40,0.60,0.75,0.90);
2050
2051 HistogramsCollectionVsEta(&plots_pfpi_res_eta[k], etaMin, etaMax, 0.5*ptVals.at(k), 2.0*ptVals.at(k), "pfpi", 0.0, 2.0);
2052 GetEresVsEta<Track>(&plots_pfpi_res_eta[k], branchPion, branchParticlePion, 211, 0.5*ptVals.at(k), 2.0*ptVals.at(k), treeReaderPion);
2053 gr_pfpi_res_eta[k] = EresGraphVsEta(&plots_pfpi_res_eta[k]);
2054
2055 s_e = Form("#pi^{ #pm}, E = %.0f GeV",ptVals.at(k));
2056 if(ptVals.at(k) >= 1000.) s_e = Form("#pi^{ #pm}, E = %.0f TeV",ptVals.at(k)/1000.);
2057
2058 leg_pfpi_res_eta[k]->SetTextFont(132);
2059 leg_pfpi_res_eta[k]->SetHeader(s_e);
2060
2061 addResoGraph(mg_pfpi_res_eta[k], &gr_hcal_res_eta[k], leg_pfpi_res_eta[k], markerStyles.at(0), colors.at(0), "HCAL");
2062 addResoGraph(mg_pfpi_res_eta[k], &gr_trkpi_res_eta[k], leg_pfpi_res_eta[k], markerStyles.at(1), colors.at(1), "Track");
2063 addResoGraph(mg_pfpi_res_eta[k], &gr_pfpi_res_eta[k], leg_pfpi_res_eta[k], markerStyles.at(2), colors.at(2), "Particle-flow");
2064
2065 c_pfpi_res_eta[k] = new TCanvas("","", 800, 600);
2066
2067 mg_pfpi_res_eta[k]->Draw("APE");
2068 //DrawAxis(mg_pfpi_res_eta[k], leg_pfpi_res_eta[k], etaMin, etaMax, 0.1, 1000, "#eta", "(resolution in E)/E (%)", false, true);
2069 DrawAxis(mg_pfpi_res_eta[k], leg_pfpi_res_eta[k], etaMin, etaMax, 0.0, 50, "#eta", "(resolution in E)/E (%)", false, false);
2070 leg_pfpi_res_eta[k]->Draw();
2071 pave->Draw();
2072
2073 TString s_ptrange = Form("pt_%.0f_",ptVals.at(k));
2074
2075 c_pfpi_res_eta[k]->Print(pdfOutput,"pdf");
2076 c_pfpi_res_eta[k]->Print(figPath+"img_pfpi_res_"+s_ptrange+"eta.pdf","pdf");
2077 c_pfpi_res_eta[k]->Print(figPath+"img_pfpi_res_"+s_ptrange+"eta.png","png");
2078
2079 }
2080
2081
2082 /////////////////
2083 // PF - jets //
2084 /////////////////
2085
2086 TMultiGraph *mg_pfjet_res_e[n_etabins];
2087 TMultiGraph *mg_pfjet_res_eta[n_ptbins];
2088
2089 TLegend *leg_pfjet_res_e[n_etabins];
2090 TLegend *leg_pfjet_res_eta[n_ptbins];
2091
2092 TGraphErrors gr_pfjet_res_e[n_etabins];
2093 TGraphErrors gr_pfjet_res_eta[n_ptbins];
2094
2095 TGraphErrors gr_cajet_res_e[n_etabins];
2096 TGraphErrors gr_cajet_res_eta[n_ptbins];
2097
2098 std::vector<resolPlot> plots_pfjet_res_e[n_etabins], plots_pfjet_res_eta[n_ptbins];
2099 std::vector<resolPlot> plots_cajet_res_e[n_etabins], plots_cajet_res_eta[n_ptbins];
2100
2101 TCanvas *c_pfjet_res_e[n_etabins];
2102 TCanvas *c_pfjet_res_eta[n_ptbins];
2103
2104
2105 // loop over eta bins
2106 for (k = 0; k < etaVals.size()-1; k++)
2107 {
2108
2109 mg_pfjet_res_e[k] = new TMultiGraph("","");
2110 leg_pfjet_res_e[k] = new TLegend(0.40,0.70,0.90,0.90);
2111
2112 HistogramsCollection(&plots_pfjet_res_e[k], TMath::Log10(ptMin), TMath::Log10(ptMax), "pfjet");
2113 HistogramsCollection(&plots_cajet_res_e[k], TMath::Log10(ptMin), TMath::Log10(ptMax), "cajet");
2114
2115 GetJetsEres(&plots_pfjet_res_e[k], branchPFJet, branchGenJet, treeReaderJet, etaVals.at(k), etaVals.at(k+1));
2116 GetJetsEres(&plots_cajet_res_e[k], branchCaloJet, branchGenJet, treeReaderJet, etaVals.at(k), etaVals.at(k+1));
2117
2118 gr_pfjet_res_e[k] = EresGraph(&plots_pfjet_res_e[k]);
2119 gr_cajet_res_e[k] = EresGraph(&plots_cajet_res_e[k]);
2120
2121 s_etaMin = Form("%.1f",etaVals.at(k));
2122 s_etaMax = Form("%.1f",etaVals.at(k+1));
2123 s_eta = "anti-k_{T}, R = 0.4, "+ s_etaMin + " < | #eta | < " + s_etaMax;
2124
2125 leg_pfjet_res_e[k]->SetTextFont(132);
2126 leg_pfjet_res_e[k]->SetHeader(s_eta);
2127
2128 addResoGraph(mg_pfjet_res_e[k], &gr_cajet_res_e[k], leg_pfjet_res_e[k], markerStyles.at(0), colors.at(0), "Calorimeter Jets");
2129 addResoGraph(mg_pfjet_res_e[k], &gr_pfjet_res_e[k], leg_pfjet_res_e[k], markerStyles.at(1), colors.at(1), "Particle-flow Jets");
2130
2131 c_pfjet_res_e[k] = new TCanvas("","", 800, 600);
2132
2133 mg_pfjet_res_e[k]->Draw("APE");
2134 //DrawAxis(mg_pfjet_res_e[k], leg_pfjet_res_e[k], 10, ptMax, 0.5, 100, "E [GeV]", "(resolution in E)/E (%)", true, true);
2135 DrawAxis(mg_pfjet_res_e[k], leg_pfjet_res_e[k], 10, ptMax, 0.0, 30, "E [GeV]", "(resolution in E)/E (%)", true, false);
2136 leg_pfjet_res_e[k]->Draw();
2137 pave->Draw();
2138
2139 TString s_etarange = "eta_"+s_etaMin+"_"+s_etaMax+"_";
2140
2141 c_pfjet_res_e[k]->Print(pdfOutput,"pdf");
2142 c_pfjet_res_e[k]->Print(figPath+"img_pfjet_res_"+s_etarange+"e.pdf","pdf");
2143 c_pfjet_res_e[k]->Print(figPath+"img_pfjet_res_"+s_etarange+"e.png","png");
2144
2145 }
2146
2147
2148 // loop over eta bins
2149 for (k = 0; k < ptVals.size(); k++)
2150 {
2151
2152 mg_pfjet_res_eta[k] = new TMultiGraph("","");
2153 leg_pfjet_res_eta[k] = new TLegend(0.30,0.70,0.85,0.90);
2154
2155 HistogramsCollectionVsEta(&plots_pfjet_res_eta[k], etaMin, etaMax, 0.5*ptVals.at(k), 2.0*ptVals.at(k), "pfjet", 0.0, 2.0);
2156 HistogramsCollectionVsEta(&plots_cajet_res_eta[k], etaMin, etaMax, 0.5*ptVals.at(k), 2.0*ptVals.at(k), "cajet", 0.0, 2.0);
2157
2158 GetJetsEresVsEta(&plots_pfjet_res_eta[k], branchPFJet, branchGenJet, treeReaderJet, 0.5*ptVals.at(k), 2.0*ptVals.at(k));
2159 GetJetsEresVsEta(&plots_cajet_res_eta[k], branchCaloJet, branchGenJet, treeReaderJet, 0.5*ptVals.at(k), 2.0*ptVals.at(k));
2160
2161 gr_pfjet_res_eta[k] = EresGraphVsEta(&plots_pfjet_res_eta[k]);
2162 gr_cajet_res_eta[k] = EresGraphVsEta(&plots_cajet_res_eta[k]);
2163
2164 s_e = Form("anti-k_{T}, R = 0.4, jets, E = %.0f GeV",ptVals.at(k));
2165 if(ptVals.at(k) >= 1000.) s_e = Form("anti-k_{T}, R = 0.4, E = %.0f TeV",ptVals.at(k)/1000.);
2166
2167 leg_pfjet_res_eta[k]->SetTextFont(132);
2168 leg_pfjet_res_eta[k]->SetHeader(s_e);
2169
2170 addResoGraph(mg_pfjet_res_eta[k], &gr_cajet_res_eta[k], leg_pfjet_res_eta[k], markerStyles.at(0), colors.at(0), "Calorimeter Jets");
2171 addResoGraph(mg_pfjet_res_eta[k], &gr_pfjet_res_eta[k], leg_pfjet_res_eta[k], markerStyles.at(1), colors.at(1), "Particle-flow Jets");
2172
2173 c_pfjet_res_eta[k] = new TCanvas("","", 800, 600);
2174
2175 mg_pfjet_res_eta[k]->Draw("APE");
2176 //DrawAxis(mg_pfjet_res_eta[k], leg_pfjet_res_eta[k], etaMin, etaMax, 0.1, 1000, "#eta", "(resolution in E)/E (%)", false, true);
2177 DrawAxis(mg_pfjet_res_eta[k], leg_pfjet_res_eta[k], etaMin, etaMax, 0.0, 50, "#eta", "(resolution in E)/E (%)", false, false);
2178 leg_pfjet_res_eta[k]->Draw();
2179 pave->Draw();
2180
2181 TString s_ptrange = Form("pt_%.0f_",ptVals.at(k));
2182
2183 c_pfjet_res_eta[k]->Print(pdfOutput,"pdf");
2184 c_pfjet_res_eta[k]->Print(figPath+"img_pfjet_res_"+s_ptrange+"eta.pdf","pdf");
2185 c_pfjet_res_eta[k]->Print(figPath+"img_pfjet_res_"+s_ptrange+"eta.png","png");
2186
2187 }
2188
2189
2190 /////////////////////
2191 // PF Missing ET ///
2192 /////////////////////
2193
2194 TMultiGraph *mg_met_res_ht = new TMultiGraph("","");
2195 TLegend *leg_met_res_ht = new TLegend(0.60,0.22,0.90,0.42);
2196
2197 std::vector<resolPlot> plots_pfmet, plots_camet;
2198
2199 HistogramsCollection(&plots_pfmet, TMath::Log10(ptMin), TMath::Log10(ptMax), "pfMET", -500, 500);
2200 HistogramsCollection(&plots_camet, TMath::Log10(ptMin), TMath::Log10(ptMax), "caMET", -500, 500);
2201
2202 GetMetres(&plots_pfmet, branchGenScalarHT, branchMet, branchPFJet, treeReaderJet);
2203 GetMetres(&plots_camet, branchGenScalarHT, branchCaloMet, branchCaloJet, treeReaderJet);
2204
2205 TGraphErrors gr_pfmet_res_ht = MetResGraph(&plots_pfmet, true);
2206 TGraphErrors gr_camet_res_ht = MetResGraph(&plots_camet, true);
2207
2208 addResoGraph(mg_met_res_ht, &gr_camet_res_ht, leg_met_res_ht, markerStyles.at(0), colors.at(0), "Calorimeter E_{T}^{miss}");
2209 addResoGraph(mg_met_res_ht, &gr_pfmet_res_ht, leg_met_res_ht, markerStyles.at(1), colors.at(1), "Particle-flow E_{T}^{miss}");
2210
2211 TCanvas *c_met_res_ht = new TCanvas("","", 800, 600);
2212
2213 mg_met_res_ht->Draw("APE");
2214 DrawAxis(mg_met_res_ht, leg_met_res_ht, 10, 10000, 0.1, 1000, " #sum p_{T} [GeV]", "resolution in E_{x,y}^{miss} [GeV]", true, true);
2215
2216 leg_met_res_ht->Draw();
2217 pave->Draw();
2218 c_met_res_ht->Print(pdfOutput,"pdf");
2219 c_met_res_ht->Print(figPath+"img_met_res_ht.pdf","pdf");
2220 c_met_res_ht->Print(figPath+"img_met_res_ht.png","png");
2221
2222
2223 /////////////////////////////////////////
2224 // Electron Reconstruction Efficiency ///
2225 /////////////////////////////////////////
2226
2227 TMultiGraph *mg_recele_eff_pt = new TMultiGraph("","");
2228 TMultiGraph *mg_recele_eff_eta = new TMultiGraph("","");
2229
2230 TLegend *leg_recele_eff_pt = new TLegend(0.55,0.22,0.90,0.48);
2231 TLegend *leg_recele_eff_eta = new TLegend(0.55,0.22,0.90,0.48);
2232
2233 TGraphErrors gr_recele_eff_pt[n_etabins], gr_recele_eff_eta[n_ptbins];
2234 TH1D* h_recele_eff_pt, *h_recele_eff_eta;
2235
2236 // loop over eta bins
2237 for (k = 0; k < etaVals.size()-1; k++)
2238 {
2239
2240 h_recele_eff_pt = GetEffPt<Electron>(branchElectron, branchParticleElectron, "Electron", 11, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderElectron);
2241 gr_recele_eff_pt[k] = TGraphErrors(h_recele_eff_pt);
2242
2243 s_etaMin = Form("%.1f",etaVals.at(k));
2244 s_etaMax = Form("%.1f",etaVals.at(k+1));
2245
2246 s_eta = "e^{ #pm} , " + s_etaMin + " < | #eta | < "+s_etaMax;
2247
2248 gr_recele_eff_pt[k].SetName("recEff_"+s_etaMin+"_"+s_etaMax);
2249
2250 addResoGraph(mg_recele_eff_pt, &gr_recele_eff_pt[k], leg_recele_eff_pt, markerStyles.at(k), colors.at(k), s_eta);
2251 }
2252
2253 // loop over pt
2254 for (k = 0; k < ptVals.size(); k++)
2255 {
2256 h_recele_eff_eta = GetEffEta<Electron>(branchElectron, branchParticleElectron, "Electron", 11, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderElectron);
2257 gr_recele_eff_eta[k] = TGraphErrors(h_recele_eff_eta);
2258
2259 s_pt = Form("e^{ #pm} , p_{T} = %.0f GeV",ptVals.at(k));
2260 if(ptVals.at(k) >= 1000.) s_pt = Form("e^{ #pm} , p_{T} = %.0f TeV",ptVals.at(k)/1000.);
2261
2262 addResoGraph(mg_recele_eff_eta, &gr_recele_eff_eta[k], leg_recele_eff_eta, markerStyles.at(k), colors.at(k), s_pt );
2263 }
2264
2265 TCanvas *c_recele_eff_pt = new TCanvas("","", 800, 600);
2266
2267 mg_recele_eff_pt->Draw("APE");
2268 DrawAxis(mg_recele_eff_pt, leg_recele_eff_pt, ptMin, ptMax, 0.0, 100, "p_{T} [GeV]", "reconstruction efficiency (%)", true, false);
2269 leg_recele_eff_pt->Draw();
2270 pave->Draw();
2271
2272 c_recele_eff_pt->Print(pdfOutput,"pdf");
2273 c_recele_eff_pt->Print(figPath+"img_recele_eff_pt.pdf","pdf");
2274 c_recele_eff_pt->Print(figPath+"img_recele_eff_pt.png","png");
2275
2276 TCanvas *c_recele_eff_eta = new TCanvas("","", 800, 600);
2277
2278 mg_recele_eff_eta->Draw("APE");
2279 DrawAxis(mg_recele_eff_eta, leg_recele_eff_eta, etaMin, etaMax, 0.0, 100, " #eta ", "reconstruction efficiency (%)", false, false);
2280 leg_recele_eff_eta->Draw();
2281 pave->Draw();
2282
2283 c_recele_eff_eta->Print(pdfOutput,"pdf");
2284 c_recele_eff_eta->Print(figPath+"img_recele_eff_eta.pdf","pdf");
2285 c_recele_eff_eta->Print(figPath+"img_recele_eff_eta.png","png");
2286
2287
2288 /////////////////////////////////////////
2289 // Muon Reconstruction Efficiency ///
2290 /////////////////////////////////////////
2291
2292 TMultiGraph *mg_recmu_eff_pt = new TMultiGraph("","");
2293 TMultiGraph *mg_recmu_eff_eta = new TMultiGraph("","");
2294
2295 TLegend *leg_recmu_eff_pt = new TLegend(0.55,0.22,0.90,0.48);
2296 TLegend *leg_recmu_eff_eta = new TLegend(0.55,0.22,0.90,0.48);
2297
2298 TGraphErrors gr_recmu_eff_pt[n_etabins], gr_recmu_eff_eta[n_ptbins];
2299 TH1D* h_recmu_eff_pt, *h_recmu_eff_eta;
2300
2301 // loop over eta bins
2302 for (k = 0; k < etaVals.size()-1; k++)
2303 {
2304
2305 h_recmu_eff_pt = GetEffPt<Muon>(branchMuon, branchParticleMuon, "muon", 13, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderMuon);
2306 gr_recmu_eff_pt[k] = TGraphErrors(h_recmu_eff_pt);
2307
2308 s_etaMin = Form("%.1f",etaVals.at(k));
2309 s_etaMax = Form("%.1f",etaVals.at(k+1));
2310
2311 s_eta = "#mu^{ #pm} , " + s_etaMin + " < | #eta | < "+s_etaMax;
2312
2313 gr_recmu_eff_pt[k].SetName("recEff_"+s_etaMin+"_"+s_etaMax);
2314
2315 addResoGraph(mg_recmu_eff_pt, &gr_recmu_eff_pt[k], leg_recmu_eff_pt, markerStyles.at(k), colors.at(k), s_eta);
2316 }
2317
2318 // loop over pt
2319 for (k = 0; k < ptVals.size(); k++)
2320 {
2321 h_recmu_eff_eta = GetEffEta<Muon>(branchMuon, branchParticleMuon, "muon", 13, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderMuon);
2322 gr_recmu_eff_eta[k] = TGraphErrors(h_recmu_eff_eta);
2323
2324 s_pt = Form("#mu^{ #pm} , p_{T} = %.0f GeV",ptVals.at(k));
2325 if(ptVals.at(k) >= 1000.) s_pt = Form("#mu^{ #pm} , p_{T} = %.0f TeV",ptVals.at(k)/1000.);
2326
2327 addResoGraph(mg_recmu_eff_eta, &gr_recmu_eff_eta[k], leg_recmu_eff_eta, markerStyles.at(k), colors.at(k), s_pt );
2328 }
2329
2330 TCanvas *c_recmu_eff_pt = new TCanvas("","", 800, 600);
2331
2332 mg_recmu_eff_pt->Draw("APE");
2333 DrawAxis(mg_recmu_eff_pt, leg_recmu_eff_pt, ptMin, ptMax, 0.0, 100, "p_{T} [GeV]", "reconstruction efficiency (%)", true, false);
2334 leg_recmu_eff_pt->Draw();
2335 pave->Draw();
2336
2337 c_recmu_eff_pt->Print(pdfOutput,"pdf");
2338 c_recmu_eff_pt->Print(figPath+"img_recmu_eff_pt.pdf","pdf");
2339 c_recmu_eff_pt->Print(figPath+"img_recmu_eff_pt.png","png");
2340
2341 TCanvas *c_recmu_eff_eta = new TCanvas("","", 800, 600);
2342
2343 mg_recmu_eff_eta->Draw("APE");
2344 DrawAxis(mg_recmu_eff_eta, leg_recmu_eff_eta, etaMin, etaMax, 0.0, 100, " #eta ", "reconstruction efficiency (%)", false, false);
2345 leg_recmu_eff_eta->Draw();
2346 pave->Draw();
2347
2348 c_recmu_eff_eta->Print(pdfOutput,"pdf");
2349 c_recmu_eff_eta->Print(figPath+"img_recmu_eff_eta.pdf","pdf");
2350 c_recmu_eff_eta->Print(figPath+"img_recmu_eff_eta.png","png");
2351
2352
2353 /////////////////////////////////////////
2354 // Photon Reconstruction Efficiency ///
2355 /////////////////////////////////////////
2356
2357 TMultiGraph *mg_recpho_eff_pt = new TMultiGraph("","");
2358 TMultiGraph *mg_recpho_eff_eta = new TMultiGraph("","");
2359
2360 TLegend *leg_recpho_eff_pt = new TLegend(0.55,0.22,0.90,0.48);
2361 TLegend *leg_recpho_eff_eta = new TLegend(0.55,0.22,0.90,0.48);
2362
2363 TGraphErrors gr_recpho_eff_pt[n_etabins], gr_recpho_eff_eta[n_ptbins];
2364 TH1D* h_recpho_eff_pt, *h_recpho_eff_eta;
2365
2366 // loop over eta bins
2367 for (k = 0; k < etaVals.size()-1; k++)
2368 {
2369
2370 h_recpho_eff_pt = GetEffPt<Photon>(branchPhoton, branchParticlePhoton, "Photon", 22, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderPhoton);
2371 gr_recpho_eff_pt[k] = TGraphErrors(h_recpho_eff_pt);
2372
2373 s_etaMin = Form("%.1f",etaVals.at(k));
2374 s_etaMax = Form("%.1f",etaVals.at(k+1));
2375
2376 s_eta = "#gamma , " + s_etaMin + " < | #eta | < "+s_etaMax;
2377
2378 gr_recpho_eff_pt[k].SetName("recEff_"+s_etaMin+"_"+s_etaMax);
2379
2380 addResoGraph(mg_recpho_eff_pt, &gr_recpho_eff_pt[k], leg_recpho_eff_pt, markerStyles.at(k), colors.at(k), s_eta);
2381 }
2382
2383 // loop over pt
2384 for (k = 0; k < ptVals.size(); k++)
2385 {
2386 h_recpho_eff_eta = GetEffEta<Photon>(branchPhoton, branchParticlePhoton, "Photon", 22, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderPhoton);
2387 gr_recpho_eff_eta[k] = TGraphErrors(h_recpho_eff_eta);
2388
2389 s_pt = Form("#gamma , p_{T} = %.0f GeV",ptVals.at(k));
2390 if(ptVals.at(k) >= 1000.) s_pt = Form("#gamma , p_{T} = %.0f TeV",ptVals.at(k)/1000.);
2391
2392 addResoGraph(mg_recpho_eff_eta, &gr_recpho_eff_eta[k], leg_recpho_eff_eta, markerStyles.at(k), colors.at(k), s_pt );
2393 }
2394
2395 TCanvas *c_recpho_eff_pt = new TCanvas("","", 800, 600);
2396
2397 mg_recpho_eff_pt->Draw("APE");
2398 DrawAxis(mg_recpho_eff_pt, leg_recpho_eff_pt, ptMin, ptMax, 0.0, 100, "p_{T} [GeV]", "reconstruction efficiency (%)", true, false);
2399 leg_recpho_eff_pt->Draw();
2400 pave->Draw();
2401
2402 c_recpho_eff_pt->Print(pdfOutput,"pdf");
2403 c_recpho_eff_pt->Print(figPath+"img_recpho_eff_pt.pdf","pdf");
2404 c_recpho_eff_pt->Print(figPath+"img_recpho_eff_pt.png","png");
2405
2406 TCanvas *c_recpho_eff_eta = new TCanvas("","", 800, 600);
2407
2408 mg_recpho_eff_eta->Draw("APE");
2409 DrawAxis(mg_recpho_eff_eta, leg_recpho_eff_eta, etaMin, etaMax, 0.0, 100, " #eta ", "reconstruction efficiency (%)", false, false);
2410 leg_recpho_eff_eta->Draw();
2411 pave->Draw();
2412
2413 c_recpho_eff_eta->Print(pdfOutput,"pdf");
2414 c_recpho_eff_eta->Print(figPath+"img_recpho_eff_eta.pdf","pdf");
2415 c_recpho_eff_eta->Print(figPath+"img_recpho_eff_eta.png","png");
2416
2417 /////////////////////////////////////////
2418 // B-jets Efficiency/ mistag rates ///
2419 /////////////////////////////////////////
2420
2421 TMultiGraph *mg_recbjet_eff_pt = new TMultiGraph("","");
2422 TMultiGraph *mg_recbjet_eff_eta = new TMultiGraph("","");
2423
2424 TLegend *leg_recbjet_eff_pt = new TLegend(0.50,0.22,0.90,0.48);
2425 TLegend *leg_recbjet_eff_eta = new TLegend(0.50,0.22,0.90,0.48);
2426
2427 TGraphErrors gr_recbjet_eff_pt[n_etabins], gr_recbjet_eff_eta[n_ptbins];
2428 TH1D* h_recbjet_eff_pt, *h_recbjet_eff_eta;
2429
2430 // loop over eta bins
2431 for (k = 0; k < etaVals.size()-1; k++)
2432 {
2433
2434 h_recbjet_eff_pt = GetEffPt<Jet>(branchPFBJet, branchParticleBJet, "BJet", 5, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderBJet);
2435 gr_recbjet_eff_pt[k] = TGraphErrors(h_recbjet_eff_pt);
2436
2437 s_etaMin = Form("%.1f",etaVals.at(k));
2438 s_etaMax = Form("%.1f",etaVals.at(k+1));
2439
2440 s_eta = "b-jet , " + s_etaMin + " < | #eta | < "+s_etaMax;
2441
2442 gr_recbjet_eff_pt[k].SetName("recEff_"+s_etaMin+"_"+s_etaMax);
2443
2444 addResoGraph(mg_recbjet_eff_pt, &gr_recbjet_eff_pt[k], leg_recbjet_eff_pt, markerStyles.at(k), colors.at(k), s_eta);
2445 }
2446
2447 // loop over pt
2448 for (k = 0; k < ptVals.size(); k++)
2449 {
2450 h_recbjet_eff_eta = GetEffEta<Jet>(branchPFBJet, branchParticleBJet, "BJet", 5, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderBJet);
2451 gr_recbjet_eff_eta[k] = TGraphErrors(h_recbjet_eff_eta);
2452
2453 s_pt = Form("b-jet , p_{T} = %.0f GeV",ptVals.at(k));
2454 if(ptVals.at(k) >= 1000.) s_pt = Form("b-jet , p_{T} = %.0f TeV",ptVals.at(k)/1000.);
2455
2456 addResoGraph(mg_recbjet_eff_eta, &gr_recbjet_eff_eta[k], leg_recbjet_eff_eta, markerStyles.at(k), colors.at(k), s_pt );
2457 }
2458
2459 TCanvas *c_recbjet_eff_pt = new TCanvas("","", 800, 600);
2460
2461 mg_recbjet_eff_pt->Draw("APE");
2462 DrawAxis(mg_recbjet_eff_pt, leg_recbjet_eff_pt, ptMin, ptMax, 0.0, 100, "p_{T} [GeV]", "b - tag efficiency (%)", true, false);
2463 leg_recbjet_eff_pt->Draw();
2464 pave->Draw();
2465
2466 c_recbjet_eff_pt->Print(pdfOutput,"pdf");
2467 c_recbjet_eff_pt->Print(figPath+"img_recbjet_eff_pt.pdf","pdf");
2468 c_recbjet_eff_pt->Print(figPath+"img_recbjet_eff_pt.png","png");
2469
2470 TCanvas *c_recbjet_eff_eta = new TCanvas("","", 800, 600);
2471
2472 mg_recbjet_eff_eta->Draw("APE");
2473 DrawAxis(mg_recbjet_eff_eta, leg_recbjet_eff_eta, etaMin, etaMax, 0.0, 100, " #eta ", "b - tag efficiency (%)", false, false);
2474 leg_recbjet_eff_eta->Draw();
2475 pave->Draw();
2476
2477 c_recbjet_eff_eta->Print(pdfOutput,"pdf");
2478 c_recbjet_eff_eta->Print(figPath+"img_recbjet_eff_eta.pdf","pdf");
2479 c_recbjet_eff_eta->Print(figPath+"img_recbjet_eff_eta.png","png");
2480
2481
2482
2483 // ------ c - mistag ------
2484
2485 TMultiGraph *mg_recbjet_cmis_pt = new TMultiGraph("","");
2486 TMultiGraph *mg_recbjet_cmis_eta = new TMultiGraph("","");
2487
2488 TLegend *leg_recbjet_cmis_pt = new TLegend(0.50,0.64,0.90,0.90);
2489 TLegend *leg_recbjet_cmis_eta = new TLegend(0.50,0.64,0.90,0.90);
2490
2491 TGraphErrors gr_recbjet_cmis_pt[n_etabins], gr_recbjet_cmis_eta[n_ptbins];
2492 TH1D* h_recbjet_cmis_pt, *h_recbjet_cmis_eta;
2493
2494 // loop over eta bins
2495 for (k = 0; k < etaVals.size()-1; k++)
2496 {
2497
2498 h_recbjet_cmis_pt = GetEffPt<Jet>(branchPFCJet, branchParticleCJet, "CJet", 4, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderCJet);
2499 gr_recbjet_cmis_pt[k] = TGraphErrors(h_recbjet_cmis_pt);
2500
2501 s_etaMin = Form("%.1f",etaVals.at(k));
2502 s_etaMax = Form("%.1f",etaVals.at(k+1));
2503
2504 s_eta = "c-jet , " + s_etaMin + " < | #eta | < "+s_etaMax;
2505
2506 gr_recbjet_cmis_pt[k].SetName("recEff_"+s_etaMin+"_"+s_etaMax);
2507
2508 addResoGraph(mg_recbjet_cmis_pt, &gr_recbjet_cmis_pt[k], leg_recbjet_cmis_pt, markerStyles.at(k), colors.at(k), s_eta);
2509 }
2510
2511 // loop over pt
2512 for (k = 0; k < ptVals.size(); k++)
2513 {
2514 h_recbjet_cmis_eta = GetEffEta<Jet>(branchPFCJet, branchParticleCJet, "CJet", 4, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderCJet);
2515 gr_recbjet_cmis_eta[k] = TGraphErrors(h_recbjet_cmis_eta);
2516
2517 s_pt = Form("c-jet , p_{T} = %.0f GeV",ptVals.at(k));
2518 if(ptVals.at(k) >= 1000.) s_pt = Form("c-jet , p_{T} = %.0f TeV",ptVals.at(k)/1000.);
2519
2520 addResoGraph(mg_recbjet_cmis_eta, &gr_recbjet_cmis_eta[k], leg_recbjet_cmis_eta, markerStyles.at(k), colors.at(k), s_pt );
2521 }
2522
2523 TCanvas *c_recbjet_cmis_pt = new TCanvas("","", 800, 600);
2524
2525 mg_recbjet_cmis_pt->Draw("APE");
2526 DrawAxis(mg_recbjet_cmis_pt, leg_recbjet_cmis_pt, ptMin, ptMax, 0.0, 20, "p_{T} [GeV]", "c - mistag rate (%)", true, false);
2527 leg_recbjet_cmis_pt->Draw();
2528 pave->Draw();
2529
2530 c_recbjet_cmis_pt->Print(pdfOutput,"pdf");
2531 c_recbjet_cmis_pt->Print(figPath+"img_recbjet_cmis_pt.pdf","pdf");
2532 c_recbjet_cmis_pt->Print(figPath+"img_recbjet_cmis_pt.png","png");
2533
2534 TCanvas *c_recbjet_cmis_eta = new TCanvas("","", 800, 600);
2535
2536 mg_recbjet_cmis_eta->Draw("APE");
2537 DrawAxis(mg_recbjet_cmis_eta, leg_recbjet_cmis_eta, etaMin, etaMax, 0.0, 20, " #eta ", "c - mistag rate (%)", false, false);
2538 leg_recbjet_cmis_eta->Draw();
2539 pave->Draw();
2540
2541 c_recbjet_cmis_eta->Print(pdfOutput,"pdf");
2542 c_recbjet_cmis_eta->Print(figPath+"img_recbjet_cmis_eta.pdf","pdf");
2543 c_recbjet_cmis_eta->Print(figPath+"img_recbjet_cmis_eta.png","png");
2544
2545 // ------ light - mistag ------
2546
2547 TMultiGraph *mg_recbjet_lmis_pt = new TMultiGraph("","");
2548 TMultiGraph *mg_recbjet_lmis_eta = new TMultiGraph("","");
2549
2550 TLegend *leg_recbjet_lmis_pt = new TLegend(0.50,0.64,0.90,0.90);
2551 TLegend *leg_recbjet_lmis_eta = new TLegend(0.50,0.64,0.90,0.90);
2552
2553 TGraphErrors gr_recbjet_lmis_pt[n_etabins], gr_recbjet_lmis_eta[n_ptbins];
2554 TH1D* h_recbjet_lmis_pt, *h_recbjet_lmis_eta;
2555
2556 // loop over eta bins
2557 for (k = 0; k < etaVals.size()-1; k++)
2558 {
2559
2560 h_recbjet_lmis_pt = GetEffPt<Jet>(branchPFJet, branchParticleJet, "Jet", 1, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderJet);
2561 gr_recbjet_lmis_pt[k] = TGraphErrors(h_recbjet_lmis_pt);
2562
2563 s_etaMin = Form("%.1f",etaVals.at(k));
2564 s_etaMax = Form("%.1f",etaVals.at(k+1));
2565
2566 s_eta = "uds-jet , " + s_etaMin + " < | #eta | < "+s_etaMax;
2567
2568 gr_recbjet_lmis_pt[k].SetName("recEff_"+s_etaMin+"_"+s_etaMax);
2569
2570 addResoGraph(mg_recbjet_lmis_pt, &gr_recbjet_lmis_pt[k], leg_recbjet_lmis_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_recbjet_lmis_eta = GetEffEta<Jet>(branchPFJet, branchParticleJet, "Jet", 1, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderJet);
2577 gr_recbjet_lmis_eta[k] = TGraphErrors(h_recbjet_lmis_eta);
2578
2579 s_pt = Form("uds-jet , p_{T} = %.0f GeV",ptVals.at(k));
2580 if(ptVals.at(k) >= 1000.) s_pt = Form("uds-jet , p_{T} = %.0f TeV",ptVals.at(k)/1000.);
2581
2582 addResoGraph(mg_recbjet_lmis_eta, &gr_recbjet_lmis_eta[k], leg_recbjet_lmis_eta, markerStyles.at(k), colors.at(k), s_pt );
2583 }
2584
2585 TCanvas *c_recbjet_lmis_pt = new TCanvas("","", 800, 600);
2586
2587 mg_recbjet_lmis_pt->Draw("APE");
2588 DrawAxis(mg_recbjet_lmis_pt, leg_recbjet_lmis_pt, ptMin, ptMax, 0.0, 0.1, "p_{T} [GeV]", "light - mistag rate (%)", true, false);
2589 leg_recbjet_lmis_pt->Draw();
2590 pave->Draw();
2591
2592 c_recbjet_lmis_pt->Print(pdfOutput,"pdf");
2593 c_recbjet_lmis_pt->Print(figPath+"img_recbjet_lmis_pt.pdf","pdf");
2594 c_recbjet_lmis_pt->Print(figPath+"img_recbjet_lmis_pt.png","png");
2595
2596 TCanvas *c_recbjet_lmis_eta = new TCanvas("","", 800, 600);
2597
2598 mg_recbjet_lmis_eta->Draw("APE");
2599 DrawAxis(mg_recbjet_lmis_eta, leg_recbjet_lmis_eta, etaMin, etaMax, 0.0, 0.1, " #eta ", "light - mistag rate (%)", false, false);
2600 leg_recbjet_lmis_eta->Draw();
2601 pave->Draw();
2602
2603 c_recbjet_lmis_eta->Print(pdfOutput,"pdf");
2604 c_recbjet_lmis_eta->Print(figPath+"img_recbjet_lmis_eta.pdf","pdf");
2605 c_recbjet_lmis_eta->Print(figPath+"img_recbjet_lmis_eta.png","png");
2606
2607
2608 ///////////////////////////////////////////
2609 // tau-jets Efficiency/ mistag rates ///
2610 ///////////////////////////////////////////
2611
2612 TMultiGraph *mg_rectaujet_eff_pt = new TMultiGraph("","");
2613 TMultiGraph *mg_rectaujet_eff_eta = new TMultiGraph("","");
2614
2615 TLegend *leg_rectaujet_eff_pt = new TLegend(0.50,0.22,0.90,0.48);
2616 TLegend *leg_rectaujet_eff_eta = new TLegend(0.50,0.22,0.90,0.48);
2617
2618 TGraphErrors gr_rectaujet_eff_pt[n_etabins], gr_rectaujet_eff_eta[n_ptbins];
2619 TH1D* h_rectaujet_eff_pt, *h_rectaujet_eff_eta;
2620
2621 // loop over eta bins
2622 for (k = 0; k < etaVals.size()-1; k++)
2623 {
2624
2625 h_rectaujet_eff_pt = GetTauEffPt<Jet>(branchPFTauJet, branchParticleTauJet, "TauJet", 15, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderTauJet);
2626 gr_rectaujet_eff_pt[k] = TGraphErrors(h_rectaujet_eff_pt);
2627
2628 s_etaMin = Form("%.1f",etaVals.at(k));
2629 s_etaMax = Form("%.1f",etaVals.at(k+1));
2630
2631 s_eta = "#tau-jet , " + s_etaMin + " < | #eta | < "+s_etaMax;
2632
2633 gr_rectaujet_eff_pt[k].SetName("recEff_"+s_etaMin+"_"+s_etaMax);
2634
2635 addResoGraph(mg_rectaujet_eff_pt, &gr_rectaujet_eff_pt[k], leg_rectaujet_eff_pt, markerStyles.at(k), colors.at(k), s_eta);
2636 }
2637
2638 // loop over pt
2639 for (k = 0; k < ptVals.size(); k++)
2640 {
2641 h_rectaujet_eff_eta = GetTauEffEta<Jet>(branchPFTauJet, branchParticleTauJet, "TauJet", 15, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderTauJet);
2642 gr_rectaujet_eff_eta[k] = TGraphErrors(h_rectaujet_eff_eta);
2643
2644 s_pt = Form("#tau-jet , p_{T} = %.0f GeV",ptVals.at(k));
2645 if(ptVals.at(k) >= 1000.) s_pt = Form("#tau-jet , p_{T} = %.0f TeV",ptVals.at(k)/1000.);
2646
2647 addResoGraph(mg_rectaujet_eff_eta, &gr_rectaujet_eff_eta[k], leg_rectaujet_eff_eta, markerStyles.at(k), colors.at(k), s_pt );
2648 }
2649
2650
2651 TCanvas *c_rectaujet_eff_pt = new TCanvas("","", 800, 600);
2652
2653 mg_rectaujet_eff_pt->Draw("APE");
2654 DrawAxis(mg_rectaujet_eff_pt, leg_rectaujet_eff_pt, ptMin, ptMax, 0.0, 100, "p_{T} [GeV]", "#tau - tag efficiency (%)", true, false);
2655 leg_rectaujet_eff_pt->Draw();
2656 pave->Draw();
2657
2658 c_rectaujet_eff_pt->Print(pdfOutput,"pdf");
2659 c_rectaujet_eff_pt->Print(figPath+"img_rectaujet_eff_pt.pdf","pdf");
2660 c_rectaujet_eff_pt->Print(figPath+"img_rectaujet_eff_pt.png","png");
2661
2662 TCanvas *c_rectaujet_eff_eta = new TCanvas("","", 800, 600);
2663
2664 mg_rectaujet_eff_eta->Draw("APE");
2665 DrawAxis(mg_rectaujet_eff_eta, leg_rectaujet_eff_eta, etaMin, etaMax, 0.0000001, 100, " #eta ", "#tau - tag efficiency (%)", false, true);
2666 leg_rectaujet_eff_eta->Draw();
2667 pave->Draw();
2668
2669 c_rectaujet_eff_eta->Print(pdfOutput,"pdf");
2670 c_rectaujet_eff_eta->Print(figPath+"img_rectaujet_eff_eta.pdf","pdf");
2671 c_rectaujet_eff_eta->Print(figPath+"img_rectaujet_eff_eta.png","png");
2672
2673
2674 //--------------- tau mistag rate ----------
2675
2676 TMultiGraph *mg_rectaujet_mis_pt = new TMultiGraph("","");
2677 TMultiGraph *mg_rectaujet_mis_eta = new TMultiGraph("","");
2678
2679 TLegend *leg_rectaujet_mis_pt = new TLegend(0.50,0.64,0.90,0.90);
2680 TLegend *leg_rectaujet_mis_eta = new TLegend(0.50,0.64,0.90,0.90);
2681
2682 TGraphErrors gr_rectaujet_mis_pt[n_etabins], gr_rectaujet_mis_eta[n_ptbins];
2683 TH1D* h_rectaujet_mis_pt, *h_rectaujet_mis_eta;
2684
2685 // loop over eta bins
2686 for (k = 0; k < etaVals.size()-1; k++)
2687 {
2688
2689 h_rectaujet_mis_pt = GetTauEffPt<Jet>(branchPFJet, branchParticleJet, "TauJet", 1, ptMin, ptMax, etaVals.at(k), etaVals.at(k+1), treeReaderJet);
2690 gr_rectaujet_mis_pt[k] = TGraphErrors(h_rectaujet_mis_pt);
2691
2692 s_etaMin = Form("%.1f",etaVals.at(k));
2693 s_etaMax = Form("%.1f",etaVals.at(k+1));
2694
2695 s_eta = "uds-jet , " + s_etaMin + " < | #eta | < "+s_etaMax;
2696
2697 gr_rectaujet_mis_pt[k].SetName("recEff_"+s_etaMin+"_"+s_etaMax);
2698
2699 addResoGraph(mg_rectaujet_mis_pt, &gr_rectaujet_mis_pt[k], leg_rectaujet_mis_pt, markerStyles.at(k), colors.at(k), s_eta);
2700 }
2701
2702 // loop over pt
2703 for (k = 0; k < ptVals.size(); k++)
2704 {
2705 h_rectaujet_mis_eta = GetTauEffEta<Jet>(branchPFJet, branchParticleJet, "TauJet", 1, 0.5*ptVals.at(k), 2.0*ptVals.at(k) ,etaMin, etaMax , treeReaderJet);
2706 gr_rectaujet_mis_eta[k] = TGraphErrors(h_rectaujet_mis_eta);
2707
2708 s_pt = Form("uds-jet , p_{T} = %.0f GeV",ptVals.at(k));
2709 if(ptVals.at(k) >= 1000.) s_pt = Form("uds-jet , p_{T} = %.0f TeV",ptVals.at(k)/1000.);
2710
2711 addResoGraph(mg_rectaujet_mis_eta, &gr_rectaujet_mis_eta[k], leg_rectaujet_mis_eta, markerStyles.at(k), colors.at(k), s_pt );
2712 }
2713
2714 TCanvas *c_rectaujet_mis_pt = new TCanvas("","", 800, 600);
2715
2716 mg_rectaujet_mis_pt->Draw("APE");
2717 DrawAxis(mg_rectaujet_mis_pt, leg_rectaujet_mis_pt, ptMin, ptMax, 0.0, 10., "p_{T} [GeV]", "#tau - mistag(%)", true, false);
2718 leg_rectaujet_mis_pt->Draw();
2719 pave->Draw();
2720
2721 c_rectaujet_mis_pt->Print(pdfOutput,"pdf");
2722 c_rectaujet_mis_pt->Print(figPath+"img_rectaujet_mis_pt.pdf","pdf");
2723 c_rectaujet_mis_pt->Print(figPath+"img_rectaujet_mis_pt.png","png");
2724
2725 TCanvas *c_rectaujet_mis_eta = new TCanvas("","", 800, 600);
2726
2727 mg_rectaujet_mis_eta->Draw("APE");
2728 DrawAxis(mg_rectaujet_mis_eta, leg_rectaujet_mis_eta, etaMin, etaMax, 0.0, 5., " #eta ", "#tau - mistag (%)", false, false);
2729 leg_rectaujet_mis_eta->Draw();
2730 pave->Draw();
2731
2732 c_rectaujet_mis_eta->Print(pdfOutput+")","pdf");
2733 c_rectaujet_mis_eta->Print(figPath+"img_rectaujet_mis_eta.pdf","pdf");
2734 c_rectaujet_mis_eta->Print(figPath+"img_rectaujet_mis_eta.png","png");
2735
2736
2737// ---- store resolution histograms in the output (for leave efficiencies out) ---
2738
2739
2740 TFile *fout = new TFile(outputFile,"recreate");
2741
2742 for (int bin = 0; bin < Nbins; bin++)
2743 {
2744
2745 for (k = 0; k < etaVals.size()-1; k++)
2746 {
2747 plots_trkpi_res_pt[k].at(bin).resolHist->Write();
2748 plots_trkele_res_pt[k].at(bin).resolHist->Write();
2749 plots_trkmu_res_pt[k].at(bin).resolHist->Write();
2750 plots_ecal_res_e[k].at(bin).resolHist->Write();
2751 plots_hcal_res_e[k].at(bin).resolHist->Write();
2752 plots_pfele_res_e[k].at(bin).resolHist->Write();
2753 plots_pfpi_res_e[k].at(bin).resolHist->Write();
2754 plots_pfjet_res_e[k].at(bin).resolHist->Write();
2755 plots_cajet_res_e[k].at(bin).resolHist->Write();
2756
2757 }
2758
2759 for (k = 0; k < ptVals.size(); k++)
2760 {
2761 plots_trkpi_res_eta[k].at(bin).resolHist->Write();
2762 plots_trkele_res_eta[k].at(bin).resolHist->Write();
2763 plots_trkmu_res_eta[k].at(bin).resolHist->Write();
2764 plots_ecal_res_eta[k].at(bin).resolHist->Write();
2765 plots_hcal_res_eta[k].at(bin).resolHist->Write();
2766 plots_pfele_res_eta[k].at(bin).resolHist->Write();
2767 plots_pfpi_res_eta[k].at(bin).resolHist->Write();
2768 plots_pfjet_res_eta[k].at(bin).resolHist->Write();
2769 plots_cajet_res_eta[k].at(bin).resolHist->Write();
2770
2771 }
2772
2773 plots_pfmet.at(bin).resolHist->Write();
2774 plots_camet.at(bin).resolHist->Write();
2775
2776 }
2777
2778 fout->Write();
2779
2780 cout << "** Exiting..." << endl;
2781
2782 delete treeReaderElectron;
2783 delete treeReaderMuon;
2784 delete treeReaderPhoton;
2785 delete treeReaderJet;
2786 delete treeReaderBJet;
2787 delete treeReaderTauJet;
2788 delete chainElectron;
2789 delete chainMuon;
2790 delete chainPhoton;
2791 delete chainJet;
2792 delete chainBJet;
2793 delete chainTauJet;
2794}
2795
2796//------------------------------------------------------------------------------
2797
2798int main(int argc, char *argv[])
2799{
2800 char *appName = "Validation";
2801
2802 if(argc != 12)
2803 {
2804 cout << " Usage: " << appName << " input_file_electron input_file_muon input_file_photon input_file_jet input_file_bjet input_file_taujet output_file" << endl;
2805 cout << " input_file_pion - input file in ROOT format ('Delphes' tree)," << endl;
2806 cout << " input_file_electron - input file in ROOT format ('Delphes' tree)," << endl;
2807 cout << " input_file_muon - input file in ROOT format ('Delphes' tree)," << endl;
2808 cout << " input_file_photon - input file in ROOT format ('Delphes' tree)," << endl;
2809 cout << " input_file_neutralhadron - input file in ROOT format ('Delphes' tree)," << endl;
2810 cout << " input_file_jet - input file in ROOT format ('Delphes' tree)," << endl;
2811 cout << " input_file_bjet - input file in ROOT format ('Delphes' tree)," << endl;
2812 cout << " input_file_cjet - input file in ROOT format ('Delphes' tree)," << endl;
2813 cout << " input_file_taujet - input file in ROOT format ('Delphes' tree)," << endl;
2814 cout << " output_file - output file in ROOT format" << endl;
2815 cout << " delphes version" << endl;
2816
2817 return 1;
2818 }
2819
2820 gROOT->SetBatch();
2821
2822 int appargc = 1;
2823 char *appargv[] = {appName};
2824 TApplication app(appName, &appargc, appargv);
2825
2826 Validation(argv[1], argv[2], argv[3], argv[4], argv[5], argv[6], argv[7], argv[8], argv[9], argv[10], argv[11]);
2827}
2828
2829
2830
Note: See TracBrowser for help on using the repository browser.