Fork me on GitHub

source: git/examples/Validation.cpp@ 8f0b34c

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 8f0b34c was 32677fb, checked in by Alexandre Mertens <alexandre.mertens@…>, 8 years ago

style changes
  • 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 C_el1->Print("delphes_validation.pdf(","pdf");
917 C_el2->Print("delphes_validation.pdf","pdf");
918
919 gDirectory->cd(0);
920
921 ///////////
922 // Muons //
923 ///////////
924
925 // Reconstruction efficiency
926 int muID = 13;
927 std::pair<TH1D*,TH1D*> histos_mu = GetEff<Muon>(branchMuon, branchParticleMuon,"Muon", muID, treeReaderMuon);
928
929 // muon tracking reconstruction efficiency
930 std::pair <TH1D*,TH1D*> histos_mutrack = GetEff<Track>(branchTrackMuon, branchParticleMuon, "muonTrack", muID, treeReaderMuon);
931
932 // Muon Pt Resolution
933 std::vector<resolPlot> plots_mu;
934 HistogramsCollection(&plots_mu, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "muons");
935 GetPtres<Muon>(&plots_mu, branchMuon, branchParticleMuon, muID, treeReaderMuon);
936 TGraphErrors gr_mu = EresGraph(&plots_mu, true);
937 TGraphErrors gr_muFwd = EresGraph(&plots_mu, false);
938 gr_mu.SetName("Muon");
939 gr_muFwd.SetName("MuonFwd");
940
941 // Muon Track Energy Resolution
942 std::vector<resolPlot> plots_mutrack;
943 HistogramsCollection(&plots_mutrack, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "muonsTracks");
944 GetPtres<Track>(&plots_mutrack, branchTrackMuon, branchParticleMuon, muID, treeReaderMuon);
945 TGraphErrors gr_mutrack = EresGraph(&plots_mutrack, true);
946 TGraphErrors gr_mutrackFwd = EresGraph(&plots_mutrack, false);
947 gr_mutrackFwd.SetName("MuonTracksFwd");
948
949 // Canvas
950
951 TString muEff = "muonEff";
952 TCanvas *C_mu1 = new TCanvas(muEff,muEff, 1600, 600);
953 C_mu1->Divide(2);
954 C_mu1->cd(1);
955 TLegend *leg_mu1 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
956 leg_mu1->SetHeader("#splitline{muons}{|#eta| < 1.5}");
957 leg_mu1->AddEntry("","","");
958
959
960 gPad->SetLogx();
961 histos_mutrack.first->Draw("][");
962 addHist(histos_mutrack.first, leg_mu1, 2);
963 histos_mu.first->Draw("same ][");
964 addHist(histos_mu.first, leg_mu1, 1);
965
966 DrawAxis(histos_mutrack.first, leg_mu1, 1);
967
968 leg_mu1->Draw();
969
970 C_mu1->cd(2);
971 TLegend *leg_mu2 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
972 leg_mu2->SetHeader("#splitline{muons}{1.5 < |#eta| < 2.5}");
973 leg_mu2->AddEntry("","","");
974
975 gPad->SetLogx();
976 histos_mutrack.second->Draw("][");
977 addHist(histos_mutrack.second, leg_mu2, 2);
978 histos_mu.second->Draw("same ][");
979 addHist(histos_mu.second, leg_mu2, 1);
980
981 DrawAxis(histos_mutrack.second, leg_mu2, 1);
982 leg_mu2->Draw();
983
984 TString muRes = "muonERes";
985 TString muResFwd = "muonEResFwd";
986
987 TCanvas *C_mu = new TCanvas(muRes,muRes, 1600, 600);
988 C_mu->Divide(2);
989 C_mu->cd(1);
990 TMultiGraph *mg_mu = new TMultiGraph(muRes,muRes);
991 TLegend *leg_mu = new TLegend(topLeftLegXmin,topLeftLegYmin,topLeftLegXmax,topLeftLegYmax);
992 leg_mu->SetHeader("#splitline{muons}{|#eta| < 1.5}");
993 leg_mu->AddEntry("","","");
994
995 addGraph(mg_mu, &gr_mutrack, leg_mu, 2);
996 addGraph(mg_mu, &gr_mu, leg_mu, 1);
997
998 mg_mu->Draw("ACX");
999 leg_mu->Draw();
1000
1001 DrawAxis(mg_mu, leg_mu, 0.3, 1);
1002
1003 C_mu->cd(2);
1004 TMultiGraph *mg_muFwd = new TMultiGraph(muResFwd,muResFwd);
1005 TLegend *leg_muFwd = new TLegend(topLeftLegXmin,topLeftLegYmin,topLeftLegXmax,topLeftLegYmax);
1006 leg_muFwd->SetHeader("#splitline{muons}{1.5 < |#eta| < 2.5}");
1007 leg_muFwd->AddEntry("","","");
1008
1009 addGraph(mg_muFwd, &gr_mutrackFwd, leg_muFwd, 2);
1010 addGraph(mg_muFwd, &gr_muFwd, leg_muFwd, 1);
1011
1012 mg_muFwd->Draw("ACX");
1013 leg_muFwd->Draw();
1014
1015 DrawAxis(mg_muFwd, leg_muFwd, 0.3, 1);
1016
1017 C_mu1->Print("delphes_validation.pdf","pdf");
1018 C_mu->Print("delphes_validation.pdf","pdf");
1019
1020 gDirectory->cd(0);
1021
1022 /////////////
1023 // Photons //
1024 /////////////
1025
1026 // Reconstruction efficiency
1027 int phID = 22;
1028 std::pair<TH1D*,TH1D*> histos_ph = GetEff<Electron>(branchPhoton, branchParticlePhoton, "Photon", phID, treeReaderPhoton);
1029 std::pair<TH1D*,TH1D*> histos_phtower = GetEff<Electron>(branchTowerPhoton, branchParticlePhoton, "Photon", phID, treeReaderPhoton);
1030
1031 // Photon Energy Resolution
1032 std::vector<resolPlot> plots_ph;
1033 HistogramsCollection(&plots_ph, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "photons");
1034 GetEres<Photon>(&plots_ph, branchPhoton, branchParticlePhoton, phID, treeReaderPhoton);
1035 TGraphErrors gr_ph = EresGraph(&plots_ph, true);
1036 TGraphErrors gr_phFwd = EresGraph(&plots_ph, false);
1037 gr_ph.SetName("Photon");
1038 gr_phFwd.SetName("PhotonFwd");
1039
1040
1041 // Photon Tower Energy Resolution
1042 std::vector<resolPlot> plots_phtower;
1043 HistogramsCollection(&plots_phtower, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "photonsTower");
1044 GetEres<Tower>(&plots_phtower, branchTowerPhoton, branchParticlePhoton, phID, treeReaderPhoton);
1045 TGraphErrors gr_phtower = EresGraph(&plots_phtower, true);
1046 TGraphErrors gr_phtowerFwd = EresGraph(&plots_phtower, false);
1047 gr_phtower.SetName("PhotonTower");
1048 gr_phtowerFwd.SetName("PhotonTowerFwd");
1049
1050 // Canvas
1051
1052 TString phEff = "photonEff";
1053 TCanvas *C_ph1 = new TCanvas(phEff,phEff, 1600, 600);
1054 C_ph1->Divide(2);
1055 C_ph1->cd(1);
1056 TLegend *leg_ph1 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
1057 leg_ph1->SetHeader("#splitline{photons}{|#eta| < 1.5}");
1058 leg_ph1->AddEntry("","","");
1059
1060
1061 gPad->SetLogx();
1062 histos_phtower.first->Draw("][");
1063 addHist(histos_phtower.first, leg_ph1, 3);
1064 histos_ph.first->Draw("same ][");
1065 addHist(histos_ph.first, leg_ph1, 1);
1066
1067 DrawAxis(histos_phtower.first, leg_ph1, 1);
1068 leg_ph1->Draw();
1069
1070 C_ph1->cd(2);
1071 TLegend *leg_ph2 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
1072 leg_ph2->SetHeader("#splitline{photons}{1.5 < |#eta| < 2.5}");
1073 leg_ph2->AddEntry("","","");
1074
1075
1076 gPad->SetLogx();
1077 histos_phtower.second->Draw("][");
1078 addHist(histos_phtower.second, leg_ph2, 3);
1079 histos_ph.second->Draw("same ][");
1080 addHist(histos_ph.second, leg_ph2, 1);
1081
1082 DrawAxis(histos_phtower.second, leg_ph2, 1);
1083 leg_ph2->Draw();
1084
1085 TString phRes = "phERes";
1086 TString phResFwd = "phEResFwd";
1087
1088 TCanvas *C_ph = new TCanvas(phRes,phRes, 1600, 600);
1089 C_ph->Divide(2);
1090 C_ph->cd(1);
1091 TMultiGraph *mg_ph = new TMultiGraph(phRes,phRes);
1092 TLegend *leg_ph = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
1093 leg_ph->SetHeader("#splitline{photons}{|#eta| < 1.5}");
1094 leg_ph->AddEntry("","","");
1095
1096 addGraph(mg_ph, &gr_phtower, leg_ph, 3);
1097 addGraph(mg_ph, &gr_ph, leg_ph, 1);
1098
1099 mg_ph->Draw("ACX");
1100 leg_ph->Draw();
1101
1102 DrawAxis(mg_ph, leg_ph, 0.1);
1103
1104 C_ph->cd(2);
1105 TMultiGraph *mg_phFwd = new TMultiGraph(phResFwd,phResFwd);
1106 TLegend *leg_phFwd = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
1107 leg_phFwd->SetHeader("#splitline{photons}{1.5 < |#eta| < 2.5}");
1108 leg_phFwd->AddEntry("","","");
1109
1110 addGraph(mg_phFwd, &gr_phtowerFwd, leg_phFwd, 3);
1111 addGraph(mg_phFwd, &gr_phFwd, leg_phFwd, 1);
1112
1113 mg_phFwd->Draw("ACX");
1114 leg_phFwd->Draw();
1115
1116 DrawAxis(mg_phFwd, leg_phFwd, 0.1);
1117
1118 C_ph1->Print("delphes_validation.pdf","pdf");
1119 C_ph->Print("delphes_validation.pdf","pdf");
1120
1121 gDirectory->cd(0);
1122
1123 //////////
1124 // Jets //
1125 //////////
1126
1127 // BJets Reconstruction efficiency
1128 int bID = 5;
1129 std::pair<TH1D*,TH1D*> histos_btag = GetEff<Jet>(branchPFBJet, branchParticleBJet,"BTag", bID, treeReaderBJet);
1130
1131 // TauJets Reconstruction efficiency
1132 int tauID = 15;
1133 std::pair<TH1D*,TH1D*> histos_tautag = GetEff<Jet>(branchPFTauJet, branchParticleTauJet,"TauTag", tauID, treeReaderTauJet);
1134
1135 // PFJets Energy Resolution
1136 std::vector<resolPlot> plots_pfjets;
1137 HistogramsCollection(&plots_pfjets, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "PFJet");
1138 GetJetsEres(&plots_pfjets, branchPFJet, branchGenJet, treeReaderJet);
1139 TGraphErrors gr_pfjets = EresGraph(&plots_pfjets, true);
1140 TGraphErrors gr_pfjetsFwd = EresGraph(&plots_pfjets, false);
1141 gr_pfjets.SetName("pfJet");
1142 gr_pfjetsFwd.SetName("pfJetFwd");
1143
1144 // CaloJets Energy Resolution
1145 std::vector<resolPlot> plots_calojets;
1146 HistogramsCollection(&plots_calojets, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "CaloJet");
1147 GetJetsEres(&plots_calojets, branchCaloJet, branchGenJet, treeReaderJet);
1148 TGraphErrors gr_calojets = EresGraph(&plots_calojets, true);
1149 TGraphErrors gr_calojetsFwd = EresGraph(&plots_calojets, false);
1150 gr_calojets.SetName("caloJet");
1151 gr_calojetsFwd.SetName("caloJetFwd");
1152
1153 // MET Resolution vs HT
1154 std::vector<resolPlot> plots_met;
1155 HistogramsCollection(&plots_met, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "MET", -500, 500);
1156 GetMetres(&plots_met, branchScalarHT, branchMet, branchPFJet, treeReaderJet);
1157 TGraphErrors gr_met = EresGraph(&plots_met, true);
1158 gr_calojets.SetName("MET");
1159
1160 // Canvas
1161 TString btagEff = "btagEff";
1162 TCanvas *C_btag1 = new TCanvas(btagEff,btagEff, 1600, 600);
1163 C_btag1->Divide(2);
1164 C_btag1->cd(1);
1165 TLegend *leg_btag1 = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
1166 leg_btag1->SetHeader("#splitline{B-tagging}{|#eta| < 1.5}");
1167 leg_btag1->AddEntry("","","");
1168
1169 gPad->SetLogx();
1170 histos_btag.first->Draw();
1171 addHist(histos_btag.first, leg_btag1, 0);
1172
1173 DrawAxis(histos_btag.first, leg_btag1, 1);
1174 leg_btag1->Draw();
1175
1176 C_btag1->cd(2);
1177 TLegend *leg_btag2 = new TLegend(resLegXmin,resLegYmin,resLegXmax+0.05,resLegYmax);
1178 leg_btag2->SetHeader("#splitline{B-tagging}{1.5 < |#eta| < 2.5}");
1179 leg_btag2->AddEntry("","","");
1180
1181 gPad->SetLogx();
1182 histos_btag.second->Draw();
1183 addHist(histos_btag.second, leg_btag2, 0);
1184
1185 DrawAxis(histos_btag.second, leg_btag2, 1);
1186 leg_btag2->Draw();
1187 C_btag1->cd(0);
1188
1189 TString tautagEff = "tautagEff";
1190 TCanvas *C_tautag1 = new TCanvas(tautagEff,tautagEff, 1600, 600);
1191 C_tautag1->Divide(2);
1192 C_tautag1->cd(1);
1193 TLegend *leg_tautag1 = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
1194 leg_tautag1->SetHeader("#splitline{#tau-tagging}{|#eta| < 1.5}");
1195 leg_tautag1->AddEntry("","","");
1196
1197 gPad->SetLogx();
1198 histos_tautag.first->Draw();
1199 addHist(histos_tautag.first, leg_tautag1, 0);
1200
1201 DrawAxis(histos_tautag.first, leg_tautag1, 1);
1202 leg_tautag1->Draw();
1203
1204 C_tautag1->cd(2);
1205 TLegend *leg_tautag2 = new TLegend(resLegXmin,resLegYmin,resLegXmax+0.05,resLegYmax);
1206 leg_tautag2->SetHeader("#splitline{#tau-tagging}{1.5 < |#eta| < 2.5}");
1207 leg_tautag2->AddEntry("","","");
1208
1209 gPad->SetLogx();
1210 histos_tautag.second->Draw();
1211 addHist(histos_tautag.second, leg_tautag2, 0);
1212
1213 DrawAxis(histos_tautag.second, leg_tautag2, 1);
1214 leg_tautag2->Draw();
1215 C_tautag1->cd(0);
1216
1217 TString jetRes = "jetERes";
1218 TString jetResFwd = "jetEResFwd";
1219 TCanvas *C_jet = new TCanvas(jetRes,jetRes, 1600, 600);
1220 C_jet->Divide(2);
1221
1222 C_jet->cd(1);
1223 TMultiGraph *mg_jet = new TMultiGraph(jetRes,jetRes);
1224 TLegend *leg_jet = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
1225 leg_jet->SetHeader("#splitline{jets}{|#eta| < 1.5}");
1226 leg_jet->AddEntry("","","");
1227
1228 addGraph(mg_jet, &gr_calojets, leg_jet, 3);
1229 addGraph(mg_jet, &gr_pfjets, leg_jet, 1);
1230
1231 mg_jet->Draw("ALX");
1232 leg_jet->Draw();
1233
1234 DrawAxis(mg_jet, leg_jet, 0.5);
1235
1236 C_jet->cd(2);
1237 TMultiGraph *mg_jetFwd = new TMultiGraph(jetResFwd,jetResFwd);
1238 TLegend *leg_jetFwd = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
1239 leg_jetFwd->SetHeader("#splitline{jets}{1.5 < |#eta| < 2.5}");
1240 leg_jetFwd->AddEntry("","","");
1241
1242 addGraph(mg_jetFwd, &gr_calojetsFwd, leg_jetFwd, 3);
1243 addGraph(mg_jetFwd, &gr_pfjetsFwd, leg_jetFwd, 1);
1244
1245 mg_jetFwd->Draw("ALX");
1246 leg_jetFwd->Draw();
1247
1248 DrawAxis(mg_jetFwd, leg_jetFwd, 0.5);
1249
1250 TString metRes = "MetRes";
1251 TCanvas *C_met = new TCanvas(metRes,metRes, 800, 600);
1252
1253 TMultiGraph *mg_met = new TMultiGraph(metRes,metRes);
1254 TLegend *leg_met = new TLegend(topLeftLegXmin,topLeftLegYmin+0.2,topLeftLegXmax-0.2,topLeftLegYmax);
1255 leg_met->SetHeader("E_{T}^{miss}");
1256 leg_met->AddEntry("","","");
1257
1258
1259 addGraph(mg_met, &gr_met, leg_met, 0);
1260
1261 mg_met->Draw("ACX");
1262 leg_met->Draw();
1263
1264 DrawAxis(mg_met, leg_met, 300);
1265
1266 mg_met->GetXaxis()->SetTitle("H_{T} [GeV]");
1267 mg_met->GetYaxis()->SetTitle("#sigma(ME_{x}) [GeV]");
1268
1269 C_jet->Print("delphes_validation.pdf","pdf");
1270 C_btag1->Print("delphes_validation.pdf","pdf");
1271 C_tautag1->Print("delphes_validation.pdf","pdf");
1272 C_met->Print("delphes_validation.pdf)","pdf");
1273
1274 TFile *fout = new TFile(outputFile,"recreate");
1275
1276 for (int bin = 0; bin < Nbins; bin++)
1277 {
1278 plots_el.at(bin).cenResolHist->Write();
1279 plots_eltrack.at(bin).cenResolHist->Write();
1280 plots_eltower.at(bin).cenResolHist->Write();
1281 plots_el.at(bin).fwdResolHist->Write();
1282 plots_eltrack.at(bin).fwdResolHist->Write();
1283 plots_eltower.at(bin).fwdResolHist->Write();
1284 }
1285
1286 histos_el.first->Write();
1287 histos_el.second->Write();
1288 histos_eltrack.first->Write();
1289 histos_eltrack.second->Write();
1290 histos_eltower.first->Write();
1291 histos_eltower.second->Write();
1292 C_el1->Write();
1293 C_el2->Write();
1294
1295 for (int bin = 0; bin < Nbins; bin++)
1296 {
1297 plots_mu.at(bin).cenResolHist->Write();
1298 plots_mutrack.at(bin).cenResolHist->Write();
1299 plots_mu.at(bin).fwdResolHist->Write();
1300 plots_mutrack.at(bin).fwdResolHist->Write();
1301 }
1302
1303 histos_mu.first->Write();
1304 histos_mu.second->Write();
1305 histos_mutrack.first->Write();
1306 histos_mutrack.second->Write();
1307 C_mu1->Write();
1308 C_mu->Write();
1309
1310 histos_ph.first->Write();
1311 histos_ph.second->Write();
1312 C_ph1->Write();
1313 C_ph->Write();
1314
1315 for (int bin = 0; bin < Nbins; bin++)
1316 {
1317 plots_pfjets.at(bin).cenResolHist->Write();
1318 plots_pfjets.at(bin).fwdResolHist->Write();
1319 plots_calojets.at(bin).cenResolHist->Write();
1320 plots_calojets.at(bin).fwdResolHist->Write();
1321 plots_met.at(bin).cenResolHist->Write();
1322 }
1323 histos_btag.first->Write();
1324 histos_btag.second->Write();
1325 histos_tautag.first->Write();
1326 histos_tautag.second->Write();
1327 C_btag1->Write();
1328 C_tautag1->Write();
1329 C_jet->Write();
1330 C_met->Write();
1331
1332 fout->Write();
1333
1334 cout << "** Exiting..." << endl;
1335
1336 delete treeReaderElectron;
1337 delete treeReaderMuon;
1338 delete treeReaderPhoton;
1339 delete treeReaderJet;
1340 delete treeReaderBJet;
1341 delete treeReaderTauJet;
1342 delete chainElectron;
1343 delete chainMuon;
1344 delete chainPhoton;
1345 delete chainJet;
1346 delete chainBJet;
1347 delete chainTauJet;
1348}
1349
1350//------------------------------------------------------------------------------
1351
1352int main(int argc, char *argv[])
1353{
1354 char *appName = "Validation";
1355
1356 if(argc != 8)
1357 {
1358 cout << " Usage: " << appName << " input_file_electron input_file_muon input_file_photon input_file_jet input_file_bjet input_file_taujet output_file" << endl;
1359 cout << " input_file_electron - input file in ROOT format ('Delphes' tree)," << endl;
1360 cout << " input_file_muon - input file in ROOT format ('Delphes' tree)," << endl;
1361 cout << " input_file_photon - input file in ROOT format ('Delphes' tree)," << endl;
1362 cout << " input_file_jet - input file in ROOT format ('Delphes' tree)," << endl;
1363 cout << " input_file_bjet - input file in ROOT format ('Delphes' tree)," << endl;
1364 cout << " input_file_taujet - input file in ROOT format ('Delphes' tree)," << endl;
1365 cout << " output_file - output file in ROOT format" << endl;
1366 return 1;
1367 }
1368
1369 gROOT->SetBatch();
1370
1371 int appargc = 1;
1372 char *appargv[] = {appName};
1373 TApplication app(appName, &appargc, appargv);
1374
1375 Validation(argv[1], argv[2], argv[3], argv[4], argv[5], argv[6], argv[7]);
1376}
1377
1378
Note: See TracBrowser for help on using the repository browser.