Fork me on GitHub

source: git/validation/DelphesValidation.cpp@ e1e72fd

Last change on this file since e1e72fd was e09b9bf, checked in by Pavel Demin <pavel.demin@…>, 4 years ago

add include for TF1.h to DelphesValidation.cpp

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