Fork me on GitHub

source: git/examples/Validation.cpp@ fa068d3

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

added card and number_of_events parameters, and some instructions
  • Property mode set to 100644
File size: 41.5 KB
Line 
1/*
2 * Delphes: a framework for fast simulation of a generic collider experiment
3 * Copyright (C) 2012-2014 Universite catholique de Louvain (UCL), Belgium
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19
20#include <iostream>
21#include <utility>
22#include <vector>
23#include <typeinfo>
24
25#include "TROOT.h"
26#include "TSystem.h"
27#include "TApplication.h"
28
29#include "TString.h"
30
31#include "TH1.h"
32#include "TH2.h"
33#include "TMath.h"
34#include "TStyle.h"
35#include "TGraph.h"
36#include "TCanvas.h"
37#include "THStack.h"
38#include "TLegend.h"
39#include "TPaveText.h"
40#include "TClonesArray.h"
41#include "TLorentzVector.h"
42#include "TGraphErrors.h"
43#include "TMultiGraph.h"
44
45#include "classes/DelphesClasses.h"
46
47#include "ExRootAnalysis/ExRootTreeReader.h"
48#include "ExRootAnalysis/ExRootTreeWriter.h"
49#include "ExRootAnalysis/ExRootTreeBranch.h"
50#include "ExRootAnalysis/ExRootResult.h"
51#include "ExRootAnalysis/ExRootUtilities.h"
52
53using namespace std;
54
55//------------------------------------------------------------------------------
56
57double ptrangemin = 10;
58double ptrangemax = 10000;
59static const int Nbins = 20;
60
61int objStyle = 1;
62int trackStyle = 7;
63int towerStyle = 3;
64
65Color_t objColor = kBlack;
66Color_t trackColor = kBlack;
67Color_t towerColor = kBlack;
68
69double effLegXmin = 0.22;
70double effLegXmax = 0.7;
71double effLegYmin = 0.22;
72double effLegYmax = 0.5;
73
74double resLegXmin = 0.62;
75double resLegXmax = 0.9;
76double resLegYmin = 0.52;
77double resLegYmax = 0.85;
78
79double topLeftLegXmin = 0.22;
80double topLeftLegXmax = 0.7;
81double topLeftLegYmin = 0.52;
82double topLeftLegYmax = 0.85;
83
84
85struct resolPlot
86{
87 TH1 *cenResolHist;
88 TH1 *fwdResolHist;
89 double ptmin;
90 double ptmax;
91 double xmin;
92 double xmax;
93 TString obj;
94
95 resolPlot();
96 resolPlot(double ptdown, double ptup, TString object);
97 void set(double ptdown, double ptup, TString object, double xmin = 0, double xmax = 2);
98 void print(){std::cout << ptmin << std::endl;}
99};
100
101
102resolPlot::resolPlot()
103{
104}
105
106resolPlot::resolPlot(double ptdown, double ptup, TString object)
107{
108 this->set(ptdown,ptup,object);
109}
110
111void resolPlot::set(double ptdown, double ptup, TString object, double xmin, double xmax)
112{
113 ptmin = ptdown;
114 ptmax = ptup;
115 obj = object;
116
117 cenResolHist = new TH1D(obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_cen", obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_cen", 200, xmin, xmax);
118 fwdResolHist = new TH1D(obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_fwd", obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_fwd", 200, xmin, xmax);
119}
120
121void HistogramsCollection(std::vector<resolPlot> *histos, double ptmin, double ptmax, TString obj, double xmin = 0, double xmax = 2)
122{
123 double width;
124 double ptdown;
125 double ptup;
126 resolPlot ptemp;
127
128 for (int i = 0; i < Nbins; i++)
129 {
130 width = (ptmax - ptmin) / Nbins;
131 ptdown = TMath::Power(10,ptmin + i * width );
132 ptup = TMath::Power(10,ptmin + (i+1) * width );
133 ptemp.set(ptdown, ptup, obj, xmin, xmax);
134 histos->push_back(ptemp);
135 }
136}
137
138//------------------------------------------------------------------------------
139
140class ExRootResult;
141class ExRootTreeReader;
142
143//------------------------------------------------------------------------------
144
145void BinLogX(TH1*h)
146{
147 TAxis *axis = h->GetXaxis();
148 int bins = axis->GetNbins();
149
150 Axis_t from = axis->GetXmin();
151 Axis_t to = axis->GetXmax();
152 Axis_t width = (to - from) / bins;
153 Axis_t *new_bins = new Axis_t[bins + 1];
154
155 for (int i = 0; i <= bins; i++)
156 {
157 new_bins[i] = TMath::Power(10, from + i * width);
158 }
159 axis->Set(bins, new_bins);
160 delete new_bins;
161}
162
163
164//------------------------------------------------------------------------------
165
166template<typename T>
167std::pair<TH1D*, TH1D*> GetEff(TClonesArray *branchReco, TClonesArray *branchParticle, TString name, int pdgID, ExRootTreeReader *treeReader)
168{
169
170 cout << "** Computing Efficiency of reconstructing "<< branchReco->GetName() << " induced by " << branchParticle->GetName() << " with PID " << pdgID << endl;
171
172 Long64_t allEntries = treeReader->GetEntries();
173
174 GenParticle *particle;
175 T *recoObj;
176
177 TLorentzVector recoMomentum, genMomentum, bestRecoMomentum;
178
179 Float_t deltaR;
180 Float_t pt, eta;
181 Long64_t entry;
182
183 Int_t i, j;
184
185 TH1D *histGenPtcen = new TH1D(name+" gen spectra Pt",name+" gen spectra cen", Nbins, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax));
186 TH1D *histRecoPtcen = new TH1D(name+" reco spectra Pt",name+" reco spectra cen", Nbins, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax));
187 TH1D *histGenPtfwd = new TH1D(name+" gen spectra Eta",name+" gen spectra fwd", Nbins, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax));
188 TH1D *histRecoPtfwd = new TH1D(name+" reco spectra Eta",name+" reco spectra fwd", Nbins, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax));
189
190 histGenPtcen->SetDirectory(0);
191 histRecoPtcen->SetDirectory(0);
192 histGenPtfwd->SetDirectory(0);
193 histRecoPtfwd->SetDirectory(0);
194
195 BinLogX(histGenPtcen);
196 BinLogX(histRecoPtcen);
197 BinLogX(histGenPtfwd);
198 BinLogX(histRecoPtfwd);
199
200 // Loop over all events
201 for(entry = 0; entry < allEntries; ++entry)
202 {
203 // Load selected branches with data from specified event
204 treeReader->ReadEntry(entry);
205
206 // Loop over all generated particle in event
207 for(i = 0; i < branchParticle->GetEntriesFast(); ++i)
208 {
209
210 particle = (GenParticle*) branchParticle->At(i);
211 genMomentum = particle->P4();
212
213 deltaR = 999;
214
215 if (particle->PID == pdgID && genMomentum.Pt() > ptrangemin && genMomentum.Pt() < ptrangemax )
216 {
217
218 // Loop over all reco object in event
219 for(j = 0; j < branchReco->GetEntriesFast(); ++j)
220 {
221 recoObj = (T*)branchReco->At(j);
222 recoMomentum = recoObj->P4();
223 // this is simply to avoid warnings from initial state particle
224 // having infite rapidity ...
225 //if(Momentum.Px() == 0 && genMomentum.Py() == 0) continue;
226
227 // take the closest parton candidate
228 if(TMath::Abs(pdgID) == 5)
229 {
230 Jet *jet = (Jet *)recoObj;
231 if(jet->BTag != 1) continue;
232 }
233 if(TMath::Abs(pdgID) == 15)
234 {
235 Jet *jet = (Jet *)recoObj;
236 if(jet->TauTag != 1) continue;
237 }
238 if(genMomentum.DeltaR(recoMomentum) < deltaR)
239 {
240 deltaR = genMomentum.DeltaR(recoMomentum);
241 bestRecoMomentum = recoMomentum;
242 }
243 }
244
245 pt = genMomentum.Pt();
246 eta = genMomentum.Eta();
247
248 if (TMath::Abs(eta) < 1.5)
249 {
250 histGenPtcen->Fill(pt);
251 if(deltaR < 0.3) { histRecoPtcen->Fill(pt); }
252 }
253 else if (TMath::Abs(eta) < 2.5)
254 {
255 histGenPtfwd->Fill(pt);
256 if(deltaR < 0.3) { histRecoPtfwd->Fill(pt); }
257
258 }
259 }
260 }
261 }
262
263
264 std::pair<TH1D*,TH1D*> histos;
265
266 histRecoPtcen->Divide(histGenPtcen);
267 histRecoPtfwd->Divide(histGenPtfwd);
268
269 histos.first = histRecoPtcen;
270 histos.second = histRecoPtfwd;
271
272 return histos;
273}
274
275template<typename T>
276void GetEres(std::vector<resolPlot> *histos, TClonesArray *branchReco, TClonesArray *branchParticle, int pdgID, ExRootTreeReader *treeReader)
277{
278 Long64_t allEntries = treeReader->GetEntries();
279
280 cout << "** Computing resolution of " << branchReco->GetName() << " induced by " << branchParticle->GetName() << " with PID " << pdgID << endl;
281
282 GenParticle *particle;
283 T* recoObj;
284
285 TLorentzVector recoMomentum, genMomentum, bestGenMomentum;
286
287 Float_t deltaR;
288 Float_t pt, eta;
289 Long64_t entry;
290
291 Int_t i, j, bin;
292
293 // Loop over all events
294 for(entry = 0; entry < allEntries; ++entry)
295 {
296 // Load selected branches with data from specified event
297 treeReader->ReadEntry(entry);
298
299 // Loop over all reconstructed jets in event
300 for(i = 0; i < branchReco->GetEntriesFast(); ++i)
301 {
302 recoObj = (T*) branchReco->At(i);
303 recoMomentum = recoObj->P4();
304
305 deltaR = 999;
306
307 // Loop over all hard partons in event
308 for(j = 0; j < branchParticle->GetEntriesFast(); ++j)
309 {
310 particle = (GenParticle*) branchParticle->At(j);
311 if (particle->PID == pdgID && particle->Status == 1)
312 {
313 genMomentum = particle->P4();
314
315 // this is simply to avoid warnings from initial state particle
316 // having infite rapidity ...
317 if(genMomentum.Px() == 0 && genMomentum.Py() == 0) continue;
318
319 // take the closest parton candidate
320 if(genMomentum.DeltaR(recoMomentum) < deltaR)
321 {
322 deltaR = genMomentum.DeltaR(recoMomentum);
323 bestGenMomentum = genMomentum;
324 }
325 }
326 }
327
328 if(deltaR < 0.3)
329 {
330 pt = bestGenMomentum.E();
331 eta = TMath::Abs(bestGenMomentum.Eta());
332
333 for (bin = 0; bin < Nbins; bin++)
334 {
335 if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta < 2.5)
336 {
337 if (eta < 1.5) {histos->at(bin).cenResolHist->Fill(recoMomentum.E()/bestGenMomentum.E());}
338 else if (eta < 2.5) {histos->at(bin).fwdResolHist->Fill(recoMomentum.E()/bestGenMomentum.E());}
339 }
340 }
341 }
342 }
343 }
344}
345
346
347template<typename T>
348void GetPtres(std::vector<resolPlot> *histos, TClonesArray *branchReco, TClonesArray *branchParticle, int pdgID, ExRootTreeReader *treeReader)
349{
350 Long64_t allEntries = treeReader->GetEntries();
351
352 cout << "** Computing pt resolution of " << branchReco->GetName() << " induced by " << branchParticle->GetName() << " with PID " << pdgID << endl;
353
354 GenParticle *particle;
355 T* recoObj;
356
357 TLorentzVector recoMomentum, genMomentum, bestGenMomentum;
358
359 Float_t deltaR;
360 Float_t pt, eta;
361 Long64_t entry;
362
363 Int_t i, j, bin;
364
365 // Loop over all events
366 for(entry = 0; entry < allEntries; ++entry)
367 {
368 // Load selected branches with data from specified event
369 treeReader->ReadEntry(entry);
370
371 // Loop over all reconstructed jets in event
372 for(i = 0; i < branchReco->GetEntriesFast(); ++i)
373 {
374 recoObj = (T*) branchReco->At(i);
375 recoMomentum = recoObj->P4();
376
377 deltaR = 999;
378
379 // Loop over all hard partons in event
380 for(j = 0; j < branchParticle->GetEntriesFast(); ++j)
381 {
382 particle = (GenParticle*) branchParticle->At(j);
383 if (particle->PID == pdgID && particle->Status == 1)
384 {
385 genMomentum = particle->P4();
386
387 // this is simply to avoid warnings from initial state particle
388 // having infite rapidity ...
389 if(genMomentum.Px() == 0 && genMomentum.Py() == 0) continue;
390
391 // take the closest parton candidate
392 if(genMomentum.DeltaR(recoMomentum) < deltaR)
393 {
394 deltaR = genMomentum.DeltaR(recoMomentum);
395 bestGenMomentum = genMomentum;
396 }
397 }
398 }
399
400 if(deltaR < 0.3)
401 {
402 pt = bestGenMomentum.Pt();
403 eta = TMath::Abs(bestGenMomentum.Eta());
404
405 for (bin = 0; bin < Nbins; bin++)
406 {
407 if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta < 2.5)
408 {
409 if (eta < 1.5) {histos->at(bin).cenResolHist->Fill(recoMomentum.Pt()/bestGenMomentum.Pt());}
410 else if (eta < 2.5) {histos->at(bin).fwdResolHist->Fill(recoMomentum.Pt()/bestGenMomentum.Pt());}
411 }
412 }
413 }
414 }
415 }
416}
417
418
419void GetJetsEres(std::vector<resolPlot> *histos, TClonesArray *branchJet, TClonesArray *branchGenJet, ExRootTreeReader *treeReader)
420{
421
422 Long64_t allEntries = treeReader->GetEntries();
423
424 cout << "** Computing resolution of " << branchJet->GetName() << " induced by " << branchGenJet->GetName() << endl;
425
426 Jet *jet, *genjet;
427
428 TLorentzVector jetMomentum, genJetMomentum, bestGenJetMomentum;
429
430 Float_t deltaR;
431 Float_t pt, eta;
432 Long64_t entry;
433
434 Int_t i, j, bin;
435
436 // Loop over all events
437 for(entry = 0; entry < allEntries; ++entry)
438 {
439 // Load selected branches with data from specified event
440 treeReader->ReadEntry(entry);
441
442 if(entry%10000 == 0) cout << "Event number: "<< entry <<endl;
443
444 // Loop over all reconstructed jets in event
445 for(i = 0; i < TMath::Min(2,branchJet->GetEntriesFast()); ++i) //branchJet->GetEntriesFast(); ++i)
446 {
447
448 jet = (Jet*) branchJet->At(i);
449 jetMomentum = jet->P4();
450
451 deltaR = 999;
452
453 // Loop over all hard partons in event
454 for(j = 0; j < TMath::Min(2,branchGenJet->GetEntriesFast()); ++j)
455 {
456 genjet = (Jet*) branchGenJet->At(j);
457
458 genJetMomentum = genjet->P4();
459
460 // this is simply to avoid warnings from initial state particle
461 // having infite rapidity ...
462 if(genJetMomentum.Px() == 0 && genJetMomentum.Py() == 0) continue;
463
464 // take the closest parton candidate
465 if(genJetMomentum.DeltaR(jetMomentum) < deltaR)
466 {
467 deltaR = genJetMomentum.DeltaR(jetMomentum);
468 bestGenJetMomentum = genJetMomentum;
469 }
470 }
471
472 if(deltaR < 0.25)
473 {
474 pt = genJetMomentum.E();
475 eta = TMath::Abs(genJetMomentum.Eta());
476
477 for (bin = 0; bin < Nbins; bin++)
478 {
479 if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta < 1.5)
480 {
481 histos->at(bin).cenResolHist->Fill(jetMomentum.E()/bestGenJetMomentum.E());
482 }
483 else if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta < 2.5)
484 {
485 histos->at(bin).fwdResolHist->Fill(jetMomentum.E()/bestGenJetMomentum.E());
486 }
487 }
488 }
489 }
490 }
491}
492
493void GetMetres(std::vector<resolPlot> *histos, TClonesArray *branchScalarHT, TClonesArray *branchMet, TClonesArray *branchJet, ExRootTreeReader *treeReader)
494{
495
496 Long64_t allEntries = treeReader->GetEntries();
497
498 cout << "** Computing resolution of " << branchMet->GetName() << " vs " << branchScalarHT->GetName() << endl;
499
500 MissingET *met;
501 ScalarHT *scalarHT;
502
503 Long64_t entry;
504
505 Int_t bin;
506 Double_t ht;
507
508 Jet *jet;
509 TLorentzVector p1, p2;
510
511 // Loop over all events
512 for(entry = 0; entry < allEntries; ++entry)
513 {
514 // Load selected branches with data from specified event
515 treeReader->ReadEntry(entry);
516
517 if(entry%10000 == 0) cout << "Event number: "<< entry <<endl;
518
519 if (branchJet->GetEntriesFast() > 1)
520 {
521
522 jet = (Jet*) branchJet->At(0);
523 p1 = jet->P4();
524 jet = (Jet*) branchJet->At(1);
525 p2 = jet->P4();
526
527 met = (MissingET*) branchMet->At(0);
528 scalarHT = (ScalarHT*) branchScalarHT->At(0);
529 ht = scalarHT->HT;
530
531 if(p1.Pt() < 0.75*ht/2) continue;
532 if(p2.Pt() < 0.75*ht/2) continue;
533
534 for (bin = 0; bin < Nbins; bin++)
535 {
536 if(ht > histos->at(bin).ptmin && ht < histos->at(bin).ptmax )
537 {
538 histos->at(bin).cenResolHist->Fill(met->P4().Px());
539 histos->at(bin).fwdResolHist->Fill(met->P4().Py());
540 }
541 }
542 }
543 }
544}
545
546
547std::pair<Double_t, Double_t> GausFit(TH1* hist)
548{
549 TF1 *f1 = new TF1("f1", "gaus", hist->GetMean()-2*hist->GetRMS(), hist->GetMean()+2*hist->GetRMS());
550 hist->Fit("f1","RQ");
551
552 TF1 *f2 = new TF1("f2", "gaus", f1->GetParameter(1) - 2*f1->GetParameter(2), f1->GetParameter(1) + 2*f1->GetParameter(2));
553 hist->Fit("f2","RQ");
554
555 Double_t sig = f2->GetParameter(2);
556 Double_t sigErr = f2->GetParError(2);
557
558 delete f1;
559 delete f2;
560 return make_pair (sig, sigErr);
561}
562
563
564TGraphErrors EresGraph(std::vector<resolPlot> *histos, bool central, bool rms = false)
565{
566 Int_t bin;
567 Int_t count = 0;
568 TGraphErrors gr = TGraphErrors(Nbins/2);
569 Double_t sig = 0;
570 Double_t sigErr = 0;
571 for (bin = 0; bin < Nbins; bin++)
572 {
573 if (central == true && histos->at(bin).cenResolHist->GetEntries() > 100)
574 {
575 std::pair<Double_t, Double_t> sigvalues = GausFit(histos->at(bin).cenResolHist);
576 if (rms == true)
577 {
578 gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.second);
579 gr.SetPointError(count,0, sigvalues.second); // to correct
580 }
581 else
582 {
583 gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.first);
584 gr.SetPointError(count,0, sigvalues.second);
585 }
586 count++;
587 }
588
589 else if (central == false && histos->at(bin).fwdResolHist->GetEntries() > 10)
590 {
591 std::pair<Double_t, Double_t> sigvalues = GausFit(histos->at(bin).fwdResolHist);
592 if (rms == true)
593 {
594 gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.second);
595 gr.SetPointError(count,0, sigvalues.second); // to correct
596 }
597 else
598 {
599 gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.first);
600 gr.SetPointError(count,0, sigvalues.second);
601 }
602 count++;
603 }
604
605 }
606 return gr;
607}
608
609
610//------------------------------------------------------------------------------
611
612
613// type 1 : object, 2 : track, 3 : tower
614
615void addGraph(TMultiGraph *mg, TGraphErrors *gr, TLegend *leg, int type)
616{
617
618 gr->SetLineWidth(2);
619
620 switch ( type )
621 {
622 case 1:
623 gr->SetLineColor(objColor);
624 gr->SetLineStyle(objStyle);
625 std::cout << "Adding " << gr->GetName() << std::endl;
626 mg->Add(gr);
627 leg->AddEntry(gr,"Reco","l");
628 break;
629
630 case 2:
631 gr->SetLineColor(trackColor);
632 gr->SetLineStyle(trackStyle);
633 mg->Add(gr);
634 leg->AddEntry(gr,"Track","l");
635 break;
636
637 case 3:
638 gr->SetLineColor(towerColor);
639 gr->SetLineStyle(towerStyle);
640 mg->Add(gr);
641 leg->AddEntry(gr,"Tower","l");
642 break;
643
644 case 0:
645 gr->SetLineColor(objColor);
646 gr->SetLineStyle(objStyle);
647 mg->Add(gr);
648 break;
649
650 default:
651 std::cout << "wrong type, possibles choices are Object, Track and Tower" << std::endl;
652 break;
653 }
654}
655
656void addHist(TH1D *h, TLegend *leg, int type)
657{
658 h->SetLineWidth(2);
659
660 switch ( type )
661 {
662 case 1:
663 h->SetLineColor(objColor);
664 h->SetLineStyle(objStyle);
665 leg->AddEntry(h,"Reco","l");
666 break;
667
668 case 2:
669 h->SetLineColor(trackColor);
670 h->SetLineStyle(trackStyle);
671 leg->AddEntry(h,"Track","l");
672 break;
673
674 case 3:
675 h->SetLineColor(towerColor);
676 h->SetLineStyle(towerStyle);
677 leg->AddEntry(h,"Tower","l");
678 break;
679
680 case 0:
681 h->SetLineColor(objColor);
682 h->SetLineStyle(objStyle);
683 break;
684
685 default:
686 std::cout << "wrong type, possibles choices are Object, Track and Tower" << std::endl;
687 break;
688 }
689}
690
691void DrawAxis(TMultiGraph *mg, TLegend *leg, double max, int type = 0)
692{
693 mg->SetMinimum(0.);
694 mg->SetMaximum(max);
695 mg->GetXaxis()->SetLimits(ptrangemin,ptrangemax);
696 mg->GetYaxis()->SetTitle("resolution");
697 if (type == 0) mg->GetXaxis()->SetTitle("E [GeV]");
698 else mg->GetXaxis()->SetTitle("p_{T} [GeV]");
699 mg->GetYaxis()->SetTitleSize(0.07);
700 mg->GetXaxis()->SetTitleSize(0.07);
701 mg->GetYaxis()->SetLabelSize(0.06);
702 mg->GetXaxis()->SetLabelSize(0.06);
703 mg->GetYaxis()->SetLabelOffset(0.03);
704 mg->GetYaxis()->SetTitleOffset(1.4);
705 mg->GetXaxis()->SetTitleOffset(1.4);
706
707 mg->GetYaxis()->SetNdivisions(505);
708
709 leg->SetBorderSize(0);
710 leg->SetShadowColor(0);
711 leg->SetFillColor(0);
712 leg->SetFillStyle(0);
713
714 gStyle->SetOptTitle(0);
715 gPad->SetLogx();
716 gPad->SetBottomMargin(0.2);
717 gPad->SetLeftMargin(0.2);
718 gPad->Modified();
719 gPad->Update();
720
721}
722
723void DrawAxis(TH1D *h, TLegend *leg, int type = 0)
724{
725
726 h->GetYaxis()->SetRangeUser(0,1.0);
727 if (type == 0) h->GetXaxis()->SetTitle("E [GeV]");
728 else h->GetXaxis()->SetTitle("p_{T} [GeV]");
729 h->GetYaxis()->SetTitle("efficiency");
730 h->GetYaxis()->SetTitleSize(0.07);
731 h->GetXaxis()->SetTitleSize(0.07);
732 h->GetYaxis()->SetLabelSize(0.06);
733 h->GetXaxis()->SetLabelSize(0.06);
734 h->GetYaxis()->SetLabelOffset(0.03);
735 h->GetYaxis()->SetTitleOffset(1.3);
736 h->GetXaxis()->SetTitleOffset(1.4);
737
738 h->GetYaxis()->SetNdivisions(505);
739
740 leg->SetBorderSize(0);
741 leg->SetShadowColor(0);
742 leg->SetFillColor(0);
743 leg->SetFillStyle(0);
744
745 gStyle->SetOptTitle(0);
746 gStyle->SetOptStat(0);
747 gPad->SetBottomMargin(0.2);
748 gPad->SetLeftMargin(0.2);
749
750 gPad->Modified();
751 gPad->Update();
752
753}
754
755
756void Validation(const char *inputFileElectron, const char *inputFileMuon, const char *inputFilePhoton, const char *inputFileJet, const char *inputFileBJet, const char *inputFileTauJet, const char *outputFile)
757{
758 TChain *chainElectron = new TChain("Delphes");
759 chainElectron->Add(inputFileElectron);
760 ExRootTreeReader *treeReaderElectron = new ExRootTreeReader(chainElectron);
761
762 TChain *chainMuon = new TChain("Delphes");
763 chainMuon->Add(inputFileMuon);
764 ExRootTreeReader *treeReaderMuon = new ExRootTreeReader(chainMuon);
765
766 TChain *chainPhoton = new TChain("Delphes");
767 chainPhoton->Add(inputFilePhoton);
768 ExRootTreeReader *treeReaderPhoton = new ExRootTreeReader(chainPhoton);
769
770 TChain *chainJet = new TChain("Delphes");
771 chainJet->Add(inputFileJet);
772 ExRootTreeReader *treeReaderJet = new ExRootTreeReader(chainJet);
773
774 TChain *chainBJet = new TChain("Delphes");
775 chainBJet->Add(inputFileBJet);
776 ExRootTreeReader *treeReaderBJet = new ExRootTreeReader(chainBJet);
777
778 TChain *chainTauJet = new TChain("Delphes");
779 chainTauJet->Add(inputFileTauJet);
780 ExRootTreeReader *treeReaderTauJet = new ExRootTreeReader(chainTauJet);
781
782 TClonesArray *branchParticleElectron = treeReaderElectron->UseBranch("Particle");
783 TClonesArray *branchTrackElectron = treeReaderElectron->UseBranch("Track");
784 TClonesArray *branchTowerElectron = treeReaderElectron->UseBranch("Tower");
785 TClonesArray *branchElectron = treeReaderElectron->UseBranch("Electron");
786
787 TClonesArray *branchParticleMuon = treeReaderMuon->UseBranch("Particle");
788 TClonesArray *branchTrackMuon = treeReaderMuon->UseBranch("Track");
789 TClonesArray *branchMuon = treeReaderMuon->UseBranch("Muon");
790
791 TClonesArray *branchParticlePhoton = treeReaderPhoton->UseBranch("Particle");
792 TClonesArray *branchTowerPhoton = treeReaderPhoton->UseBranch("Tower");
793 TClonesArray *branchPhoton = treeReaderPhoton->UseBranch("Photon");
794
795 TClonesArray *branchGenJet = treeReaderJet->UseBranch("GenJet");
796 TClonesArray *branchPFJet = treeReaderJet->UseBranch("Jet");
797 TClonesArray *branchCaloJet = treeReaderJet->UseBranch("CaloJet");
798
799 TClonesArray *branchParticleBJet = treeReaderBJet->UseBranch("Particle");
800 TClonesArray *branchPFBJet = treeReaderBJet->UseBranch("Jet");
801
802 TClonesArray *branchParticleTauJet = treeReaderTauJet->UseBranch("Particle");
803 TClonesArray *branchPFTauJet = treeReaderTauJet->UseBranch("Jet");
804
805 TClonesArray *branchScalarHT = treeReaderJet->UseBranch("ScalarHT");
806 TClonesArray *branchMet = treeReaderJet->UseBranch("MissingET");
807
808 ///////////////
809 // Electrons //
810 ///////////////
811
812 // Reconstruction efficiency
813 TString elecs = "Electron";
814 int elID = 11;
815 std::pair<TH1D*,TH1D*> histos_el = GetEff<Electron>(branchElectron, branchParticleElectron, "Electron", elID, treeReaderElectron);
816
817 // tracking reconstruction efficiency
818 std::pair <TH1D*,TH1D*> histos_eltrack = GetEff<Track>(branchTrackElectron, branchParticleElectron, "electronTrack", elID, treeReaderElectron);
819
820 // Tower reconstruction efficiency
821 std::pair <TH1D*,TH1D*> histos_eltower = GetEff<Tower>(branchTowerElectron, branchParticleElectron, "electronTower", elID, treeReaderElectron);
822
823 // Electron Energy Resolution
824 std::vector<resolPlot> plots_el;
825 HistogramsCollection(&plots_el, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "electrons");
826 GetEres<Electron>(&plots_el, branchElectron, branchParticleElectron, elID, treeReaderElectron);
827 TGraphErrors gr_el = EresGraph(&plots_el, true);
828 TGraphErrors gr_elFwd = EresGraph(&plots_el, false);
829 gr_el.SetName("Electron");
830 gr_elFwd.SetName("ElectronFwd");
831
832 // Electron Track Energy Resolution
833 std::vector<resolPlot> plots_eltrack;
834 HistogramsCollection(&plots_eltrack, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "electronsTracks");
835 GetEres<Track>(&plots_eltrack, branchTrackElectron, branchParticleElectron, elID, treeReaderElectron);
836 TGraphErrors gr_eltrack = EresGraph(&plots_eltrack, true);
837 TGraphErrors gr_eltrackFwd = EresGraph(&plots_eltrack, false);
838 gr_eltrack.SetName("ElectronTracks");
839 gr_eltrackFwd.SetName("ElectronTracksFwd");
840
841 // Electron Tower Energy Resolution
842 std::vector<resolPlot> plots_eltower;
843 HistogramsCollection(&plots_eltower, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "electronsTower");
844 GetEres<Tower>(&plots_eltower, branchTowerElectron, branchParticleElectron, elID, treeReaderElectron);
845 TGraphErrors gr_eltower = EresGraph(&plots_eltower, true);
846 TGraphErrors gr_eltowerFwd = EresGraph(&plots_eltower, false);
847 gr_eltower.SetName("ElectronTower");
848 gr_eltrackFwd.SetName("ElectronTracksFwd");
849
850 // Canvases
851 TString elEff = "electronEff";
852 TCanvas *C_el1 = new TCanvas(elEff,elEff, 1600, 600);
853 C_el1->Divide(2);
854 C_el1->cd(1);
855 TLegend *leg_el1 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
856 leg_el1->SetHeader("#splitline{electrons}{|#eta| < 1.5}");
857 leg_el1->AddEntry("","","");
858
859 gPad->SetLogx();
860 histos_eltrack.first->Draw("][");
861 addHist(histos_eltrack.first, leg_el1, 2);
862 histos_el.first->Draw("same ][");
863 addHist(histos_el.first, leg_el1, 1);
864 DrawAxis(histos_eltrack.first, leg_el1,1);
865
866 leg_el1->Draw();
867
868 C_el1->cd(2);
869 TLegend *leg_el2 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
870 leg_el2->SetHeader("#splitline{electrons}{1.5 < |#eta| < 2.5}");
871 leg_el2->AddEntry("","","");
872
873 gPad->SetLogx();
874 histos_eltrack.second->Draw("][");
875 addHist(histos_eltrack.second, leg_el2, 2);
876 histos_el.second->Draw("same ][");
877 addHist(histos_el.second, leg_el2, 1);
878
879 DrawAxis(histos_eltrack.second, leg_el2, 1);
880 leg_el2->Draw();
881
882 TString elRes = "electronERes";
883 TString elResFwd = "electronEResForward";
884 TCanvas *C_el2 = new TCanvas(elRes,elRes, 1600, 600);
885 C_el2->Divide(2);
886 C_el2->cd(1);
887 TMultiGraph *mg_el = new TMultiGraph(elRes,elRes);
888 TLegend *leg_el = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
889 leg_el->SetHeader("#splitline{electrons}{|#eta| < 1.5}");
890 leg_el->AddEntry("","","");
891
892 addGraph(mg_el, &gr_eltower, leg_el, 3);
893 addGraph(mg_el, &gr_eltrack, leg_el, 2);
894 addGraph(mg_el, &gr_el, leg_el, 1);
895
896 mg_el->Draw("ACX");
897 leg_el->Draw();
898
899 DrawAxis(mg_el, leg_el, 0.1);
900
901 C_el2->cd(2);
902 TMultiGraph *mg_elFwd = new TMultiGraph(elResFwd,elResFwd);
903 TLegend *leg_elFwd = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
904 leg_elFwd->SetHeader("#splitline{electrons}{1.5 < |#eta| < 2.5}");
905 leg_elFwd->AddEntry("","","");
906
907 addGraph(mg_elFwd, &gr_eltowerFwd, leg_elFwd, 3);
908 addGraph(mg_elFwd, &gr_eltrackFwd, leg_elFwd, 2);
909 addGraph(mg_elFwd, &gr_elFwd, leg_elFwd, 1);
910
911 mg_elFwd->Draw("ACX");
912 leg_elFwd->Draw();
913
914 DrawAxis(mg_elFwd, leg_elFwd, 0.2);
915
916 TString pdfOutput(outputFile);
917 pdfOutput.ReplaceAll(".root", ".pdf");
918
919 C_el1->Print(pdfOutput+"(","pdf");
920 C_el2->Print(pdfOutput,"pdf");
921
922 gDirectory->cd(0);
923
924 ///////////
925 // Muons //
926 ///////////
927
928 // Reconstruction efficiency
929 int muID = 13;
930 std::pair<TH1D*,TH1D*> histos_mu = GetEff<Muon>(branchMuon, branchParticleMuon,"Muon", muID, treeReaderMuon);
931
932 // muon tracking reconstruction efficiency
933 std::pair <TH1D*,TH1D*> histos_mutrack = GetEff<Track>(branchTrackMuon, branchParticleMuon, "muonTrack", muID, treeReaderMuon);
934
935 // Muon Pt Resolution
936 std::vector<resolPlot> plots_mu;
937 HistogramsCollection(&plots_mu, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "muons");
938 GetPtres<Muon>(&plots_mu, branchMuon, branchParticleMuon, muID, treeReaderMuon);
939 TGraphErrors gr_mu = EresGraph(&plots_mu, true);
940 TGraphErrors gr_muFwd = EresGraph(&plots_mu, false);
941 gr_mu.SetName("Muon");
942 gr_muFwd.SetName("MuonFwd");
943
944 // Muon Track Energy Resolution
945 std::vector<resolPlot> plots_mutrack;
946 HistogramsCollection(&plots_mutrack, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "muonsTracks");
947 GetPtres<Track>(&plots_mutrack, branchTrackMuon, branchParticleMuon, muID, treeReaderMuon);
948 TGraphErrors gr_mutrack = EresGraph(&plots_mutrack, true);
949 TGraphErrors gr_mutrackFwd = EresGraph(&plots_mutrack, false);
950 gr_mutrackFwd.SetName("MuonTracksFwd");
951
952 // Canvas
953
954 TString muEff = "muonEff";
955 TCanvas *C_mu1 = new TCanvas(muEff,muEff, 1600, 600);
956 C_mu1->Divide(2);
957 C_mu1->cd(1);
958 TLegend *leg_mu1 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
959 leg_mu1->SetHeader("#splitline{muons}{|#eta| < 1.5}");
960 leg_mu1->AddEntry("","","");
961
962
963 gPad->SetLogx();
964 histos_mutrack.first->Draw("][");
965 addHist(histos_mutrack.first, leg_mu1, 2);
966 histos_mu.first->Draw("same ][");
967 addHist(histos_mu.first, leg_mu1, 1);
968
969 DrawAxis(histos_mutrack.first, leg_mu1, 1);
970
971 leg_mu1->Draw();
972
973 C_mu1->cd(2);
974 TLegend *leg_mu2 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
975 leg_mu2->SetHeader("#splitline{muons}{1.5 < |#eta| < 2.5}");
976 leg_mu2->AddEntry("","","");
977
978 gPad->SetLogx();
979 histos_mutrack.second->Draw("][");
980 addHist(histos_mutrack.second, leg_mu2, 2);
981 histos_mu.second->Draw("same ][");
982 addHist(histos_mu.second, leg_mu2, 1);
983
984 DrawAxis(histos_mutrack.second, leg_mu2, 1);
985 leg_mu2->Draw();
986
987 TString muRes = "muonERes";
988 TString muResFwd = "muonEResFwd";
989
990 TCanvas *C_mu = new TCanvas(muRes,muRes, 1600, 600);
991 C_mu->Divide(2);
992 C_mu->cd(1);
993 TMultiGraph *mg_mu = new TMultiGraph(muRes,muRes);
994 TLegend *leg_mu = new TLegend(topLeftLegXmin,topLeftLegYmin,topLeftLegXmax,topLeftLegYmax);
995 leg_mu->SetHeader("#splitline{muons}{|#eta| < 1.5}");
996 leg_mu->AddEntry("","","");
997
998 addGraph(mg_mu, &gr_mutrack, leg_mu, 2);
999 addGraph(mg_mu, &gr_mu, leg_mu, 1);
1000
1001 mg_mu->Draw("ACX");
1002 leg_mu->Draw();
1003
1004 DrawAxis(mg_mu, leg_mu, 0.3, 1);
1005
1006 C_mu->cd(2);
1007 TMultiGraph *mg_muFwd = new TMultiGraph(muResFwd,muResFwd);
1008 TLegend *leg_muFwd = new TLegend(topLeftLegXmin,topLeftLegYmin,topLeftLegXmax,topLeftLegYmax);
1009 leg_muFwd->SetHeader("#splitline{muons}{1.5 < |#eta| < 2.5}");
1010 leg_muFwd->AddEntry("","","");
1011
1012 addGraph(mg_muFwd, &gr_mutrackFwd, leg_muFwd, 2);
1013 addGraph(mg_muFwd, &gr_muFwd, leg_muFwd, 1);
1014
1015 mg_muFwd->Draw("ACX");
1016 leg_muFwd->Draw();
1017
1018 DrawAxis(mg_muFwd, leg_muFwd, 0.3, 1);
1019
1020 C_mu1->Print(pdfOutput,"pdf");
1021 C_mu->Print(pdfOutput,"pdf");
1022
1023 gDirectory->cd(0);
1024
1025 /////////////
1026 // Photons //
1027 /////////////
1028
1029 // Reconstruction efficiency
1030 int phID = 22;
1031 std::pair<TH1D*,TH1D*> histos_ph = GetEff<Electron>(branchPhoton, branchParticlePhoton, "Photon", phID, treeReaderPhoton);
1032 std::pair<TH1D*,TH1D*> histos_phtower = GetEff<Electron>(branchTowerPhoton, branchParticlePhoton, "Photon", phID, treeReaderPhoton);
1033
1034 // Photon Energy Resolution
1035 std::vector<resolPlot> plots_ph;
1036 HistogramsCollection(&plots_ph, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "photons");
1037 GetEres<Photon>(&plots_ph, branchPhoton, branchParticlePhoton, phID, treeReaderPhoton);
1038 TGraphErrors gr_ph = EresGraph(&plots_ph, true);
1039 TGraphErrors gr_phFwd = EresGraph(&plots_ph, false);
1040 gr_ph.SetName("Photon");
1041 gr_phFwd.SetName("PhotonFwd");
1042
1043
1044 // Photon Tower Energy Resolution
1045 std::vector<resolPlot> plots_phtower;
1046 HistogramsCollection(&plots_phtower, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "photonsTower");
1047 GetEres<Tower>(&plots_phtower, branchTowerPhoton, branchParticlePhoton, phID, treeReaderPhoton);
1048 TGraphErrors gr_phtower = EresGraph(&plots_phtower, true);
1049 TGraphErrors gr_phtowerFwd = EresGraph(&plots_phtower, false);
1050 gr_phtower.SetName("PhotonTower");
1051 gr_phtowerFwd.SetName("PhotonTowerFwd");
1052
1053 // Canvas
1054
1055 TString phEff = "photonEff";
1056 TCanvas *C_ph1 = new TCanvas(phEff,phEff, 1600, 600);
1057 C_ph1->Divide(2);
1058 C_ph1->cd(1);
1059 TLegend *leg_ph1 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
1060 leg_ph1->SetHeader("#splitline{photons}{|#eta| < 1.5}");
1061 leg_ph1->AddEntry("","","");
1062
1063
1064 gPad->SetLogx();
1065 histos_phtower.first->Draw("][");
1066 addHist(histos_phtower.first, leg_ph1, 3);
1067 histos_ph.first->Draw("same ][");
1068 addHist(histos_ph.first, leg_ph1, 1);
1069
1070 DrawAxis(histos_phtower.first, leg_ph1, 1);
1071 leg_ph1->Draw();
1072
1073 C_ph1->cd(2);
1074 TLegend *leg_ph2 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
1075 leg_ph2->SetHeader("#splitline{photons}{1.5 < |#eta| < 2.5}");
1076 leg_ph2->AddEntry("","","");
1077
1078
1079 gPad->SetLogx();
1080 histos_phtower.second->Draw("][");
1081 addHist(histos_phtower.second, leg_ph2, 3);
1082 histos_ph.second->Draw("same ][");
1083 addHist(histos_ph.second, leg_ph2, 1);
1084
1085 DrawAxis(histos_phtower.second, leg_ph2, 1);
1086 leg_ph2->Draw();
1087
1088 TString phRes = "phERes";
1089 TString phResFwd = "phEResFwd";
1090
1091 TCanvas *C_ph = new TCanvas(phRes,phRes, 1600, 600);
1092 C_ph->Divide(2);
1093 C_ph->cd(1);
1094 TMultiGraph *mg_ph = new TMultiGraph(phRes,phRes);
1095 TLegend *leg_ph = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
1096 leg_ph->SetHeader("#splitline{photons}{|#eta| < 1.5}");
1097 leg_ph->AddEntry("","","");
1098
1099 addGraph(mg_ph, &gr_phtower, leg_ph, 3);
1100 addGraph(mg_ph, &gr_ph, leg_ph, 1);
1101
1102 mg_ph->Draw("ACX");
1103 leg_ph->Draw();
1104
1105 DrawAxis(mg_ph, leg_ph, 0.1);
1106
1107 C_ph->cd(2);
1108 TMultiGraph *mg_phFwd = new TMultiGraph(phResFwd,phResFwd);
1109 TLegend *leg_phFwd = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
1110 leg_phFwd->SetHeader("#splitline{photons}{1.5 < |#eta| < 2.5}");
1111 leg_phFwd->AddEntry("","","");
1112
1113 addGraph(mg_phFwd, &gr_phtowerFwd, leg_phFwd, 3);
1114 addGraph(mg_phFwd, &gr_phFwd, leg_phFwd, 1);
1115
1116 mg_phFwd->Draw("ACX");
1117 leg_phFwd->Draw();
1118
1119 DrawAxis(mg_phFwd, leg_phFwd, 0.1);
1120
1121 C_ph1->Print(pdfOutput,"pdf");
1122 C_ph->Print(pdfOutput,"pdf");
1123
1124 gDirectory->cd(0);
1125
1126 //////////
1127 // Jets //
1128 //////////
1129
1130 // BJets Reconstruction efficiency
1131 int bID = 5;
1132 std::pair<TH1D*,TH1D*> histos_btag = GetEff<Jet>(branchPFBJet, branchParticleBJet,"BTag", bID, treeReaderBJet);
1133
1134 // TauJets Reconstruction efficiency
1135 int tauID = 15;
1136 std::pair<TH1D*,TH1D*> histos_tautag = GetEff<Jet>(branchPFTauJet, branchParticleTauJet,"TauTag", tauID, treeReaderTauJet);
1137
1138 // PFJets Energy Resolution
1139 std::vector<resolPlot> plots_pfjets;
1140 HistogramsCollection(&plots_pfjets, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "PFJet");
1141 GetJetsEres(&plots_pfjets, branchPFJet, branchGenJet, treeReaderJet);
1142 TGraphErrors gr_pfjets = EresGraph(&plots_pfjets, true);
1143 TGraphErrors gr_pfjetsFwd = EresGraph(&plots_pfjets, false);
1144 gr_pfjets.SetName("pfJet");
1145 gr_pfjetsFwd.SetName("pfJetFwd");
1146
1147 // CaloJets Energy Resolution
1148 std::vector<resolPlot> plots_calojets;
1149 HistogramsCollection(&plots_calojets, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "CaloJet");
1150 GetJetsEres(&plots_calojets, branchCaloJet, branchGenJet, treeReaderJet);
1151 TGraphErrors gr_calojets = EresGraph(&plots_calojets, true);
1152 TGraphErrors gr_calojetsFwd = EresGraph(&plots_calojets, false);
1153 gr_calojets.SetName("caloJet");
1154 gr_calojetsFwd.SetName("caloJetFwd");
1155
1156 // MET Resolution vs HT
1157 std::vector<resolPlot> plots_met;
1158 HistogramsCollection(&plots_met, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "MET", -500, 500);
1159 GetMetres(&plots_met, branchScalarHT, branchMet, branchPFJet, treeReaderJet);
1160 TGraphErrors gr_met = EresGraph(&plots_met, true);
1161 gr_calojets.SetName("MET");
1162
1163 // Canvas
1164 TString btagEff = "btagEff";
1165 TCanvas *C_btag1 = new TCanvas(btagEff,btagEff, 1600, 600);
1166 C_btag1->Divide(2);
1167 C_btag1->cd(1);
1168 TLegend *leg_btag1 = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
1169 leg_btag1->SetHeader("#splitline{B-tagging}{|#eta| < 1.5}");
1170 leg_btag1->AddEntry("","","");
1171
1172 gPad->SetLogx();
1173 histos_btag.first->Draw();
1174 addHist(histos_btag.first, leg_btag1, 0);
1175
1176 DrawAxis(histos_btag.first, leg_btag1, 1);
1177 leg_btag1->Draw();
1178
1179 C_btag1->cd(2);
1180 TLegend *leg_btag2 = new TLegend(resLegXmin,resLegYmin,resLegXmax+0.05,resLegYmax);
1181 leg_btag2->SetHeader("#splitline{B-tagging}{1.5 < |#eta| < 2.5}");
1182 leg_btag2->AddEntry("","","");
1183
1184 gPad->SetLogx();
1185 histos_btag.second->Draw();
1186 addHist(histos_btag.second, leg_btag2, 0);
1187
1188 DrawAxis(histos_btag.second, leg_btag2, 1);
1189 leg_btag2->Draw();
1190 C_btag1->cd(0);
1191
1192 TString tautagEff = "tautagEff";
1193 TCanvas *C_tautag1 = new TCanvas(tautagEff,tautagEff, 1600, 600);
1194 C_tautag1->Divide(2);
1195 C_tautag1->cd(1);
1196 TLegend *leg_tautag1 = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
1197 leg_tautag1->SetHeader("#splitline{#tau-tagging}{|#eta| < 1.5}");
1198 leg_tautag1->AddEntry("","","");
1199
1200 gPad->SetLogx();
1201 histos_tautag.first->Draw();
1202 addHist(histos_tautag.first, leg_tautag1, 0);
1203
1204 DrawAxis(histos_tautag.first, leg_tautag1, 1);
1205 leg_tautag1->Draw();
1206
1207 C_tautag1->cd(2);
1208 TLegend *leg_tautag2 = new TLegend(resLegXmin,resLegYmin,resLegXmax+0.05,resLegYmax);
1209 leg_tautag2->SetHeader("#splitline{#tau-tagging}{1.5 < |#eta| < 2.5}");
1210 leg_tautag2->AddEntry("","","");
1211
1212 gPad->SetLogx();
1213 histos_tautag.second->Draw();
1214 addHist(histos_tautag.second, leg_tautag2, 0);
1215
1216 DrawAxis(histos_tautag.second, leg_tautag2, 1);
1217 leg_tautag2->Draw();
1218 C_tautag1->cd(0);
1219
1220 TString jetRes = "jetERes";
1221 TString jetResFwd = "jetEResFwd";
1222 TCanvas *C_jet = new TCanvas(jetRes,jetRes, 1600, 600);
1223 C_jet->Divide(2);
1224
1225 C_jet->cd(1);
1226 TMultiGraph *mg_jet = new TMultiGraph(jetRes,jetRes);
1227 TLegend *leg_jet = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
1228 leg_jet->SetHeader("#splitline{jets}{|#eta| < 1.5}");
1229 leg_jet->AddEntry("","","");
1230
1231 addGraph(mg_jet, &gr_calojets, leg_jet, 3);
1232 addGraph(mg_jet, &gr_pfjets, leg_jet, 1);
1233
1234 mg_jet->Draw("ALX");
1235 leg_jet->Draw();
1236
1237 DrawAxis(mg_jet, leg_jet, 0.5);
1238
1239 C_jet->cd(2);
1240 TMultiGraph *mg_jetFwd = new TMultiGraph(jetResFwd,jetResFwd);
1241 TLegend *leg_jetFwd = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
1242 leg_jetFwd->SetHeader("#splitline{jets}{1.5 < |#eta| < 2.5}");
1243 leg_jetFwd->AddEntry("","","");
1244
1245 addGraph(mg_jetFwd, &gr_calojetsFwd, leg_jetFwd, 3);
1246 addGraph(mg_jetFwd, &gr_pfjetsFwd, leg_jetFwd, 1);
1247
1248 mg_jetFwd->Draw("ALX");
1249 leg_jetFwd->Draw();
1250
1251 DrawAxis(mg_jetFwd, leg_jetFwd, 0.5);
1252
1253 TString metRes = "MetRes";
1254 TCanvas *C_met = new TCanvas(metRes,metRes, 800, 600);
1255
1256 TMultiGraph *mg_met = new TMultiGraph(metRes,metRes);
1257 TLegend *leg_met = new TLegend(topLeftLegXmin,topLeftLegYmin+0.2,topLeftLegXmax-0.2,topLeftLegYmax);
1258 leg_met->SetHeader("E_{T}^{miss}");
1259 leg_met->AddEntry("","","");
1260
1261
1262 addGraph(mg_met, &gr_met, leg_met, 0);
1263
1264 mg_met->Draw("ACX");
1265 leg_met->Draw();
1266
1267 DrawAxis(mg_met, leg_met, 300);
1268
1269 mg_met->GetXaxis()->SetTitle("H_{T} [GeV]");
1270 mg_met->GetYaxis()->SetTitle("#sigma(ME_{x}) [GeV]");
1271
1272 C_jet->Print(pdfOutput,"pdf");
1273 C_btag1->Print(pdfOutput,"pdf");
1274 C_tautag1->Print(pdfOutput,"pdf");
1275 C_met->Print(pdfOutput+")","pdf");
1276
1277 TFile *fout = new TFile(outputFile,"recreate");
1278
1279 for (int bin = 0; bin < Nbins; bin++)
1280 {
1281 plots_el.at(bin).cenResolHist->Write();
1282 plots_eltrack.at(bin).cenResolHist->Write();
1283 plots_eltower.at(bin).cenResolHist->Write();
1284 plots_el.at(bin).fwdResolHist->Write();
1285 plots_eltrack.at(bin).fwdResolHist->Write();
1286 plots_eltower.at(bin).fwdResolHist->Write();
1287 }
1288
1289 histos_el.first->Write();
1290 histos_el.second->Write();
1291 histos_eltrack.first->Write();
1292 histos_eltrack.second->Write();
1293 histos_eltower.first->Write();
1294 histos_eltower.second->Write();
1295 C_el1->Write();
1296 C_el2->Write();
1297
1298 for (int bin = 0; bin < Nbins; bin++)
1299 {
1300 plots_mu.at(bin).cenResolHist->Write();
1301 plots_mutrack.at(bin).cenResolHist->Write();
1302 plots_mu.at(bin).fwdResolHist->Write();
1303 plots_mutrack.at(bin).fwdResolHist->Write();
1304 }
1305
1306 histos_mu.first->Write();
1307 histos_mu.second->Write();
1308 histos_mutrack.first->Write();
1309 histos_mutrack.second->Write();
1310 C_mu1->Write();
1311 C_mu->Write();
1312
1313 histos_ph.first->Write();
1314 histos_ph.second->Write();
1315 C_ph1->Write();
1316 C_ph->Write();
1317
1318 for (int bin = 0; bin < Nbins; bin++)
1319 {
1320 plots_pfjets.at(bin).cenResolHist->Write();
1321 plots_pfjets.at(bin).fwdResolHist->Write();
1322 plots_calojets.at(bin).cenResolHist->Write();
1323 plots_calojets.at(bin).fwdResolHist->Write();
1324 plots_met.at(bin).cenResolHist->Write();
1325 }
1326 histos_btag.first->Write();
1327 histos_btag.second->Write();
1328 histos_tautag.first->Write();
1329 histos_tautag.second->Write();
1330 C_btag1->Write();
1331 C_tautag1->Write();
1332 C_jet->Write();
1333 C_met->Write();
1334
1335 fout->Write();
1336
1337 cout << "** Exiting..." << endl;
1338
1339 delete treeReaderElectron;
1340 delete treeReaderMuon;
1341 delete treeReaderPhoton;
1342 delete treeReaderJet;
1343 delete treeReaderBJet;
1344 delete treeReaderTauJet;
1345 delete chainElectron;
1346 delete chainMuon;
1347 delete chainPhoton;
1348 delete chainJet;
1349 delete chainBJet;
1350 delete chainTauJet;
1351}
1352
1353//------------------------------------------------------------------------------
1354
1355int main(int argc, char *argv[])
1356{
1357 char *appName = "Validation";
1358
1359 if(argc != 8)
1360 {
1361 cout << " Usage: " << appName << " input_file_electron input_file_muon input_file_photon input_file_jet input_file_bjet input_file_taujet output_file" << endl;
1362 cout << " input_file_electron - input file in ROOT format ('Delphes' tree)," << endl;
1363 cout << " input_file_muon - input file in ROOT format ('Delphes' tree)," << endl;
1364 cout << " input_file_photon - input file in ROOT format ('Delphes' tree)," << endl;
1365 cout << " input_file_jet - input file in ROOT format ('Delphes' tree)," << endl;
1366 cout << " input_file_bjet - input file in ROOT format ('Delphes' tree)," << endl;
1367 cout << " input_file_taujet - input file in ROOT format ('Delphes' tree)," << endl;
1368 cout << " output_file - output file in ROOT format" << endl;
1369 return 1;
1370 }
1371
1372 gROOT->SetBatch();
1373
1374 int appargc = 1;
1375 char *appargv[] = {appName};
1376 TApplication app(appName, &appargc, appargv);
1377
1378 Validation(argv[1], argv[2], argv[3], argv[4], argv[5], argv[6], argv[7]);
1379}
1380
1381
Note: See TracBrowser for help on using the repository browser.