Fork me on GitHub

source: git/examples/Validation.cpp@ 85ad2b9

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

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