Fork me on GitHub

source: git/examples/Validation.C@ fe3908f

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

working with calos

  • Property mode set to 100644
File size: 26.2 KB
Line 
1/*
2This macro shows how to compute jet energy scale.
3root -l examples/Example4.C'("delphes_output.root", "plots.root")'
4
5The output ROOT file contains the pT(MC)/pT(Reco) distributions for
6various pT(Reco) and |eta| bins. The peak value of such distribution is
7interpreted as the jet energy correction to be applied for that
8given pT(Reco), |eta| bin.
9
10This can be done by modifying the "ScaleFormula" input parameter to
11the JetEnergyScale module in the delphes_card_XXX.tcl
12
13e.g a smooth function:
14
15 set ScaleFormula { sqrt(3.0 - 0.1*(abs(eta)))^2 / pt + 1.0) }
16
17or a binned function:
18
19 set ScaleFormula {(abs(eta) > 0.0 && abs(eta) <= 2.5) * (pt > 20.0 && pt <= 50.0) * (1.10) +
20 (abs(eta) > 0.0 && abs(eta) <= 2.5) * (pt > 50.0 && pt <= 100.0) * (1.05) +
21 (abs(eta) > 0.0 && abs(eta) <= 2.5) * (pt > 100.0) * (1.00) +
22 (abs(eta) > 2.5 && abs(eta) <= 5.0) * (pt > 20.0 && pt <= 50.0) * (1.10) +
23 (abs(eta) > 2.5 && abs(eta) <= 5.0) * (pt > 50.0 && pt <= 100.0) * (1.05) +
24 (abs(eta) > 2.5 && abs(eta) <= 5.0) * (pt > 100.0) * (1.00)}
25
26Be aware that a binned jet energy scale can produce "steps" in the corrected
27jet pt distribution ...
28*/
29
30#ifdef __CLING__
31R__LOAD_LIBRARY(libDelphes)
32#include "classes/DelphesClasses.h"
33#include "external/ExRootAnalysis/ExRootTreeReader.h"
34#include "external/ExRootAnalysis/ExRootResult.h"
35#else
36class ExRootTreeReader;
37class ExRootResult;
38#endif
39
40#include "TCanvas.h"
41#include "TSystem.h"
42#include "TStyle.h"
43#include "TLegend.h"
44#include <TH1.h>
45#include "TString.h"
46#include "vector"
47#include <TMath.h>
48#include <iostream>
49#include "TGraph.h"
50#include "TGraphErrors.h"
51#include "TMultiGraph.h"
52#include <typeinfo>
53#include "TLorentzVector.h"
54
55//------------------------------------------------------------------------------
56
57double ptrangemin = 10;
58double ptrangemax = 10000;
59static const int Nbins = 20;
60
61int objStyle = 1;
62int trackStyle = 7;
63int towerStyle = 5;
64
65Color_t objColor = kBlack;
66Color_t trackColor = kBlack;
67Color_t towerColor = kBlack;
68
69double effLegXmin = 0.22;
70double effLegXmax = 0.5;
71double effLegYmin = 0.22;
72double effLegYmax = 0.4;
73
74struct resolPlot
75{
76 TH1 *cenResolHist;
77 TH1 *fwdResolHist;
78 double ptmin;
79 double ptmax;
80 TString obj;
81
82 resolPlot();
83 resolPlot(double ptdown, double ptup, TString object);
84 void set(double ptdown, double ptup, TString object);
85 void print(){std::cout << ptmin << std::endl;}
86};
87
88
89resolPlot::resolPlot()
90{
91}
92
93resolPlot::resolPlot(double ptdown, double ptup, TString object)
94{
95 this->set(ptdown,ptup,object);
96}
97
98void resolPlot::set(double ptdown, double ptup, TString object){
99 ptmin = ptdown;
100 ptmax = ptup;
101 obj = object;
102
103 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", 500, -1, 1);
104 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", 500, 0.4, 0.4);
105
106}
107
108void HistogramsCollection(std::vector<resolPlot> *histos, double ptmin, double ptmax, TString obj)
109{
110 double width;
111 double ptdown;
112 double ptup;
113 resolPlot ptemp;
114
115 for (int i = 0; i < Nbins; i++)
116 {
117 width = (ptmax - ptmin) / Nbins;
118 ptdown = TMath::Power(10,ptmin + i * width );
119 ptup = TMath::Power(10,ptmin + (i+1) * width );
120 ptemp.set(ptdown, ptup, obj);
121 histos->push_back(ptemp);
122 }
123}
124
125//------------------------------------------------------------------------------
126
127class ExRootResult;
128class ExRootTreeReader;
129
130//------------------------------------------------------------------------------
131
132void BinLogX(TH1*h)
133{
134
135 TAxis *axis = h->GetXaxis();
136 int bins = axis->GetNbins();
137
138 Axis_t from = axis->GetXmin();
139 Axis_t to = axis->GetXmax();
140 Axis_t width = (to - from) / bins;
141 Axis_t *new_bins = new Axis_t[bins + 1];
142
143 for (int i = 0; i <= bins; i++) {
144 new_bins[i] = TMath::Power(10, from + i * width);
145
146 }
147 axis->Set(bins, new_bins);
148 delete new_bins;
149}
150
151
152//------------------------------------------------------------------------------
153
154template<typename T>
155std::pair<TH1D*, TH1D*> GetEff(TClonesArray *branchReco, TClonesArray *branchParticle, TString name, int pdgID, ExRootTreeReader *treeReader)
156{
157
158 cout << "** Computing Efficiency of reconstructing "<< branchReco->GetName() << " induced by " << branchParticle->GetName() << " with PID " << pdgID << endl;
159
160 Long64_t allEntries = treeReader->GetEntries();
161
162 GenParticle *particle;
163 T *recoObj;
164
165 TLorentzVector recoMomentum, genMomentum, bestRecoMomentum;
166
167 Float_t deltaR;
168 Float_t pt, eta;
169 Long64_t entry;
170
171 Int_t i, j;
172
173 TH1D *histGenPt = new TH1D(name+" gen spectra Pt",name+" gen spectra Pt", Nbins, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax));
174 TH1D *histRecoPt = new TH1D(name+" reco spectra Pt",name+" reco spectra Pt", Nbins, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax));
175 TH1D *histGenEta = new TH1D(name+" gen spectra Eta",name+" gen spectra Eta", 12, -3, 3);
176 TH1D *histRecoEta = new TH1D(name+" reco spectra Eta",name+" reco spectra Eta", 12, -3, 3);
177
178
179 BinLogX(histGenPt);
180 BinLogX(histRecoPt);
181
182 // Loop over all events
183 for(entry = 0; entry < allEntries; ++entry)
184 {
185 // Load selected branches with data from specified event
186 treeReader->ReadEntry(entry);
187
188 // Loop over all generated particle in event
189 for(i = 0; i < branchParticle->GetEntriesFast(); ++i)
190 {
191
192 particle = (GenParticle*) branchParticle->At(i);
193 genMomentum = particle->P4();
194
195 deltaR = 999;
196
197 if (particle->PID == pdgID && genMomentum.Pt() > ptrangemin && genMomentum.Pt() < ptrangemax )
198 {
199
200 // Loop over all reco object in event
201 for(j = 0; j < branchReco->GetEntriesFast(); ++j)
202 {
203 recoObj = (T*)branchReco->At(j);
204 recoMomentum = recoObj->P4();
205 // this is simply to avoid warnings from initial state particle
206 // having infite rapidity ...
207 //if(Momentum.Px() == 0 && genMomentum.Py() == 0) continue;
208
209 // take the closest parton candidate
210 if(genMomentum.DeltaR(recoMomentum) < deltaR)
211 {
212 deltaR = genMomentum.DeltaR(recoMomentum);
213 bestRecoMomentum = recoMomentum;
214 }
215 }
216
217 pt = genMomentum.Pt();
218 eta = genMomentum.Eta();
219
220 histGenPt->Fill(pt);
221 histGenEta->Fill(eta);
222
223 if(deltaR < 0.3)
224 {
225 histRecoPt->Fill(pt);
226 histRecoEta->Fill(eta);
227 }
228 }
229 }
230 }
231
232
233 std::pair<TH1D*,TH1D*> histos;
234
235 histRecoPt->Divide(histGenPt);
236 histRecoEta->Divide(histGenEta);
237
238 histos.first = histRecoPt;
239 histos.second = histRecoEta;
240
241 return histos;
242}
243
244template<typename T>
245void GetEres(std::vector<resolPlot> *histos, TClonesArray *branchReco, TClonesArray *branchParticle, int pdgID, ExRootTreeReader *treeReader)
246{
247 Long64_t allEntries = treeReader->GetEntries();
248
249 cout << "** Computing resolution of " << branchReco->GetName() << " induced by " << branchParticle->GetName() << " with PID " << pdgID << endl;
250
251 GenParticle *particle;
252 T* recoObj;
253
254 TLorentzVector recoMomentum, genMomentum, bestGenMomentum;
255
256 Float_t deltaR;
257 Float_t pt, eta;
258 Long64_t entry;
259
260 Int_t i, j, bin;
261
262 // Loop over all events
263 for(entry = 0; entry < allEntries; ++entry)
264 {
265 // Load selected branches with data from specified event
266 treeReader->ReadEntry(entry);
267
268 // Loop over all reconstructed jets in event
269 for(i = 0; i < branchReco->GetEntriesFast(); ++i)
270 {
271 recoObj = (T*) branchReco->At(i);
272 recoMomentum = recoObj->P4();
273
274 deltaR = 999;
275
276 // Loop over all hard partons in event
277 for(j = 0; j < branchParticle->GetEntriesFast(); ++j)
278 {
279 particle = (GenParticle*) branchParticle->At(j);
280 if (particle->PID == pdgID && particle->Status == 1)
281 {
282 genMomentum = particle->P4();
283
284 // this is simply to avoid warnings from initial state particle
285 // having infite rapidity ...
286 if(genMomentum.Px() == 0 && genMomentum.Py() == 0) continue;
287
288 // take the closest parton candidate
289 if(genMomentum.DeltaR(recoMomentum) < deltaR)
290 {
291 deltaR = genMomentum.DeltaR(recoMomentum);
292 bestGenMomentum = genMomentum;
293 }
294 }
295 }
296
297 if(deltaR < 0.3)
298 {
299 pt = bestGenMomentum.Pt();
300 eta = TMath::Abs(bestGenMomentum.Eta());
301
302 for (bin = 0; bin < Nbins; bin++)
303 {
304 if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta > 0.0 && eta < 2.5)
305 {
306 if (eta < 1.5) {histos->at(bin).cenResolHist->Fill((bestGenMomentum.E()-recoMomentum.E())/bestGenMomentum.E());}
307 else if (eta < 2.5) {histos->at(bin).fwdResolHist->Fill((bestGenMomentum.E()-recoMomentum.E())/bestGenMomentum.E());}
308 }
309 }
310 }
311 }
312 }
313}
314void GetJetsEres(std::vector<resolPlot> *histos, TClonesArray *branchJet, TClonesArray *branchGenJet, ExRootTreeReader *treeReader)
315{
316
317 Long64_t allEntries = treeReader->GetEntries();
318
319 cout << "** Computing resolution of " << branchJet->GetName() << " induced by " << branchGenJet->GetName() << endl;
320
321 Jet *jet, *genjet;
322
323 TLorentzVector jetMomentum, genJetMomentum, bestGenJetMomentum;
324
325 Float_t deltaR;
326 Float_t pt, eta;
327 Long64_t entry;
328
329 Int_t i, j, bin;
330
331 // Loop over all events
332 for(entry = 0; entry < allEntries; ++entry)
333 {
334 // Load selected branches with data from specified event
335 treeReader->ReadEntry(entry);
336
337 if(entry%10000 == 0) cout << "Event number: "<< entry <<endl;
338
339 // Loop over all reconstructed jets in event
340 for(i = 0; i < TMath::Min(2,branchJet->GetEntriesFast()); ++i) //branchJet->GetEntriesFast(); ++i)
341 {
342
343 jet = (Jet*) branchJet->At(i);
344 jetMomentum = jet->P4();
345
346 deltaR = 999;
347
348 // Loop over all hard partons in event
349 for(j = 0; j < TMath::Min(2,branchGenJet->GetEntriesFast()); ++j)
350 {
351 genjet = (Jet*) branchGenJet->At(j);
352
353 genJetMomentum = genjet->P4();
354
355 // this is simply to avoid warnings from initial state particle
356 // having infite rapidity ...
357 if(genJetMomentum.Px() == 0 && genJetMomentum.Py() == 0) continue;
358
359 // take the closest parton candidate
360 if(genJetMomentum.DeltaR(jetMomentum) < deltaR)
361 {
362 deltaR = genJetMomentum.DeltaR(jetMomentum);
363 bestGenJetMomentum = genJetMomentum;
364 }
365 }
366
367 if(deltaR < 0.3)
368 {
369 pt = genJetMomentum.Pt();
370 eta = TMath::Abs(genJetMomentum.Eta());
371
372 for (bin = 0; bin < Nbins; bin++)
373 {
374 if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta > 0.0 && eta < 2.5)
375 {
376 histos->at(bin).cenResolHist->Fill((bestGenJetMomentum.E()-jetMomentum.E())/bestGenJetMomentum.E());
377 }
378 }
379 }
380 }
381 }
382}
383std::pair<Double_t, Double_t> GausFit(TH1* hist)
384{
385 TF1 *f1 = new TF1("f1", "gaus", hist->GetMean()-2*hist->GetRMS(), hist->GetMean()+2*hist->GetRMS());
386 hist->Fit("f1","RQ");
387 Double_t sig = f1->GetParameter(2);
388 Double_t sigErr = f1->GetParError(2);
389 delete f1;
390 return make_pair (sig, sigErr);
391 //return make_pair (hist->GetRMS(), hist->GetRMSError());
392}
393
394
395TGraphErrors EresGraph(std::vector<resolPlot> *histos, bool central)
396{
397 Int_t bin;
398 Int_t count = 0;
399 TGraphErrors gr = TGraphErrors(Nbins/2);
400 Double_t sig = 0;
401 Double_t sigErr = 0;
402 for (bin = 0; bin < Nbins; bin++)
403 {
404 if (central == true && histos->at(bin).cenResolHist->GetEntries() > 100)
405 {
406 std::cout << " pt : " << (histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0;
407 std::cout << " mean : " << histos->at(bin).cenResolHist->GetMean() << " RMS : " << histos->at(bin).cenResolHist->GetRMS();
408 std::cout << " entries : " << histos->at(bin).cenResolHist->GetEntries() << std::endl;
409 std::pair<Double_t, Double_t> sigvalues = GausFit(histos->at(bin).cenResolHist);
410 gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.first);
411 gr.SetPointError(count,0, sigvalues.second);
412 count++;
413 }
414 /*
415 else if (histos->at(bin).cenResolHist->GetEntries() > 10)
416 {
417 histos->at(bin).fwdResolHist->Fit("gaus","","", -2*histos->at(bin).fwdResolHist->GetRMS(), 2*histos->at(bin).fwdResolHist->GetRMS());
418 TF1 *f = histos->at(bin).fwdResolHist->GetFunction("gaus");
419 Double_t sig = f->GetParameter(2);
420 Double_t sigErr = f->GetParError(2);
421 gr.SetPoint(bin,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sig);
422 gr.SetPointError(bin,0, sigErr);
423 }
424 */
425 }
426 return gr;
427}
428
429
430//------------------------------------------------------------------------------
431
432
433// type 1 : object, 2 : track, 3 : tower
434
435void addGraph(TMultiGraph *mg, TGraphErrors *gr, TLegend *leg, int type)
436{
437
438 gr->SetLineWidth(3);
439
440 switch ( type )
441 {
442 case 1:
443 gr->SetLineColor(objColor);
444 gr->SetLineStyle(objStyle);
445 std::cout << "Adding " << gr->GetName() << std::endl;
446 mg->Add(gr);
447 leg->AddEntry(gr,"Reco","l");
448 break;
449
450 case 2:
451 gr->SetLineColor(trackColor);
452 gr->SetLineStyle(trackStyle);
453 mg->Add(gr);
454 leg->AddEntry(gr,"Track","l");
455 break;
456
457 case 3:
458 gr->SetLineColor(towerColor);
459 gr->SetLineStyle(towerStyle);
460 mg->Add(gr);
461 leg->AddEntry(gr,"Tower","l");
462 break;
463
464 default:
465 std::cout << "wrong type, possibles choices are Object, Track and Tower" << std::endl;
466 break;
467 }
468}
469
470void addHist(TH1D *h, TLegend *leg, int type)
471{
472 h->SetLineWidth(3);
473
474 switch ( type )
475 {
476 case 1:
477 h->SetLineColor(objColor);
478 h->SetLineStyle(objStyle);
479 leg->AddEntry(h,"Reco","l");
480 break;
481
482 case 2:
483 h->SetLineColor(trackColor);
484 h->SetLineStyle(trackStyle);
485 leg->AddEntry(h,"Track","l");
486 break;
487
488 case 3:
489 h->SetLineColor(towerColor);
490 h->SetLineStyle(towerStyle);
491 leg->AddEntry(h,"Tower","l");
492 break;
493
494 default:
495 std::cout << "wrong type, possibles choices are Object, Track and Tower" << std::endl;
496 break;
497 }
498}
499
500void DrawAxis(TMultiGraph *mg, TLegend *leg, double max)
501{
502 mg->SetMinimum(0.);
503 mg->SetMaximum(max);
504 mg->GetXaxis()->SetLimits(ptrangemin,ptrangemax);
505 mg->GetYaxis()->SetTitle("#sigma (E) / E_{gen}");
506 mg->GetXaxis()->SetTitle("E_{gen}");
507 mg->GetYaxis()->SetTitleSize(0.07);
508 mg->GetXaxis()->SetTitleSize(0.07);
509 mg->GetYaxis()->SetLabelSize(0.07);
510 mg->GetXaxis()->SetLabelSize(0.07);
511
512 leg->SetLineStyle(0);
513 leg->SetFillStyle(0);
514 leg->SetLineWidth(0);
515 leg->SetLineColor(0);
516
517 gStyle->SetOptTitle(0);
518 gPad->SetLogx();
519 gPad->SetBottomMargin(0.2);
520 gPad->SetLeftMargin(0.2);
521 gPad->Modified();
522 gPad->Update();
523
524}
525
526void DrawAxis(TH1D *h, TLegend *leg, int type)
527{
528
529 h->GetYaxis()->SetRangeUser(0,1.2);
530 if (type == 0) h->GetXaxis()->SetTitle("E_{gen}");
531 else h->GetXaxis()->SetTitle("#eta");
532 h->GetYaxis()->SetTitle("#epsilon");
533 h->GetYaxis()->SetTitleSize(0.07);
534 h->GetXaxis()->SetTitleSize(0.07);
535 h->GetYaxis()->SetLabelSize(0.07);
536 h->GetXaxis()->SetLabelSize(0.07);
537 leg->SetLineStyle(0);
538 leg->SetFillStyle(0);
539 leg->SetLineWidth(0);
540 leg->SetLineColor(0);
541
542 gStyle->SetOptTitle(0);
543 gStyle->SetOptStat(0);
544 gPad->SetBottomMargin(0.2);
545 gPad->SetLeftMargin(0.2);
546
547 gPad->Modified();
548 gPad->Update();
549
550}
551
552
553void Validation(const char *inputFile, const char *outputFile)
554{
555 //gSystem->Load("libDelphes");
556
557 std::cout << "input file : " << inputFile << " " << " , output file : " << outputFile << std::endl;
558
559 TChain *chain = new TChain("Delphes");
560 chain->Add(inputFile);
561
562 ExRootTreeReader *treeReader = new ExRootTreeReader(chain);
563 TClonesArray *branchParticle = treeReader->UseBranch("Particle");
564 TClonesArray *branchElectron = treeReader->UseBranch("Electron");
565 TClonesArray *branchMuon = treeReader->UseBranch("Muon");
566 TClonesArray *branchPhoton = treeReader->UseBranch("Photon");
567 TClonesArray *branchTrack = treeReader->UseBranch("Track");
568 TClonesArray *branchTower = treeReader->UseBranch("Tower");
569 TClonesArray *branchGenJet = treeReader->UseBranch("GenJet");
570 TClonesArray *branchPFJet = treeReader->UseBranch("Jet");
571 TClonesArray *branchCaloJet = treeReader->UseBranch("CaloJet");
572
573
574
575 ///////////////
576 // Electrons //
577 ///////////////
578
579 // Reconstruction efficiency
580 TString elecs = "Electron";
581 int elID = 11;
582 std::pair<TH1D*,TH1D*> histos_el = GetEff<Electron>(branchElectron, branchParticle, "Electron", elID, treeReader);
583
584 // tracking reconstruction efficiency
585 std::pair <TH1D*,TH1D*> histos_eltrack = GetEff<Track>(branchTrack, branchParticle, "electronTrack", elID, treeReader);
586
587 // Tower reconstruction efficiency
588 std::pair <TH1D*,TH1D*> histos_eltower = GetEff<Tower>(branchTower, branchParticle, "electronTower", elID, treeReader);
589
590 // Electron Energy Resolution
591 std::vector<resolPlot> plots_el;
592 HistogramsCollection(&plots_el, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "electrons");
593 GetEres<Electron>( &plots_el, branchElectron, branchParticle, elID, treeReader);
594 TGraphErrors gr_el = EresGraph(&plots_el, true);
595 gr_el.SetName("Electron");
596
597 // Electron Track Energy Resolution
598 std::vector<resolPlot> plots_eltrack;
599 HistogramsCollection(&plots_eltrack, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "electronsTracks");
600 GetEres<Track>( &plots_eltrack, branchTrack, branchParticle, elID, treeReader);
601 TGraphErrors gr_eltrack = EresGraph(&plots_eltrack, true);
602 gr_eltrack.SetName("ElectronTracks");
603
604 // Electron Tower Energy Resolution
605 std::vector<resolPlot> plots_eltower;
606 HistogramsCollection(&plots_eltower, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "electronsTower");
607 GetEres<Tower>( &plots_eltower, branchTower, branchParticle, elID, treeReader);
608 TGraphErrors gr_eltower = EresGraph(&plots_eltower, true);
609 gr_eltower.SetName("ElectronTower");
610
611 // Canvases
612 TString elEff = "electronEff";
613 TCanvas *C_el1 = new TCanvas(elEff,elEff, 1000, 500);
614 C_el1->Divide(2);
615 C_el1->cd(1);
616 TLegend *leg_el1 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
617
618 gPad->SetLogx();
619 histos_eltrack.first->Draw();
620 addHist(histos_eltrack.first, leg_el1, 2);
621 histos_eltower.first->Draw("same");
622 addHist(histos_eltower.first, leg_el1, 3);
623 histos_el.first->Draw("same");
624 addHist(histos_el.first, leg_el1, 1);
625
626 DrawAxis(histos_eltrack.first, leg_el1, 0);
627 leg_el1->Draw();
628
629 C_el1->cd(2);
630 TLegend *leg_el2 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
631
632 histos_eltrack.second->Draw();
633 addHist(histos_eltrack.second, leg_el2, 2);
634 histos_eltower.second->Draw("same");
635 addHist(histos_eltower.second, leg_el2, 3);
636 histos_el.second->Draw("same");
637 addHist(histos_el.second, leg_el2, 1);
638
639 DrawAxis(histos_eltrack.second, leg_el2, 1);
640 delete(leg_el2);
641 C_el1->cd(0);
642
643 TString elRes = "electronERes";
644 TCanvas *C_el2 = new TCanvas(elRes,elRes, 1000, 500);
645 TMultiGraph *mg_el = new TMultiGraph(elRes,elRes);
646 TLegend *leg_el = new TLegend(0.52,0.7,0.9,0.9);
647
648 addGraph(mg_el, &gr_eltower, leg_el, 3);
649 addGraph(mg_el, &gr_eltrack, leg_el, 2);
650 addGraph(mg_el, &gr_el, leg_el, 1);
651
652 mg_el->Draw("ACX");
653 leg_el->Draw();
654
655 DrawAxis(mg_el, leg_el, 0.1);
656
657 C_el1->SaveAs(elEff+".eps");
658 C_el2->SaveAs(elRes+".eps");
659
660 gDirectory->cd(0);
661
662/*
663 ///////////
664 // Muons //
665 ///////////
666
667 // Reconstruction efficiency
668 int muID = 13;
669 std::pair<TH1D*,TH1D*> histos_mu = GetEff<Muon>(branchMuon, branchParticle,"Muon", muID, treeReader);
670
671 // muon tracking reconstruction efficiency
672 std::pair <TH1D*,TH1D*> histos_mutrack = GetEff<Track>(branchTrack, branchParticle, "muonTrack", muID, treeReader);
673
674 // Muon Energy Resolution
675 std::vector<resolPlot> plots_mu;
676 HistogramsCollection(&plots_mu, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "muons");
677 GetEres<Muon>( &plots_mu, branchMuon, branchParticle, muID, treeReader);
678 TGraphErrors gr_mu = EresGraph(&plots_mu, true);
679 gr_mu.SetName("Muon");
680
681 // Muon Track Energy Resolution
682 std::vector<resolPlot> plots_mutrack;
683 HistogramsCollection(&plots_mutrack, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "muonsTracks");
684 GetEres<Track>( &plots_mutrack, branchTrack, branchParticle, muID, treeReader);
685 TGraphErrors gr_mutrack = EresGraph(&plots_mutrack, true);
686 gr_eltrack.SetName("MuonTracks");
687
688 // Canvas
689
690 TString muEff = "muonEff";
691 TCanvas *C_mu1 = new TCanvas(muEff,muEff, 1000, 500);
692 C_mu1->Divide(2);
693 C_mu1->cd(1);
694 TLegend *leg_mu1 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
695
696 gPad->SetLogx();
697 histos_mutrack.first->Draw();
698 addHist(histos_mutrack.first, leg_mu1, 2);
699 histos_mu.first->Draw("same");
700 addHist(histos_mu.first, leg_mu1, 1);
701
702 DrawAxis(histos_mutrack.first, leg_mu1, 0);
703 leg_mu1->Draw();
704
705 C_mu1->cd(2);
706 TLegend *leg_mu2 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
707
708 histos_mutrack.second->Draw();
709 addHist(histos_mutrack.second, leg_mu2, 2);
710 histos_mu.second->Draw("same");
711 addHist(histos_mu.second, leg_mu2, 1);
712
713 DrawAxis(histos_mutrack.second, leg_mu2, 1);
714
715 TString muRes = "muonERes";
716 TCanvas *C_mu = new TCanvas(muRes,muRes, 1000, 500);
717 TMultiGraph *mg_mu = new TMultiGraph(muRes,muRes);
718 TLegend *leg_mu = new TLegend(0.52,0.7,0.9,0.9);
719
720 addGraph(mg_mu, &gr_mutrack, leg_mu, 2);
721 addGraph(mg_mu, &gr_mu, leg_mu, 1);
722
723 mg_mu->Draw("ACX");
724 leg_mu->Draw();
725
726 DrawAxis(mg_mu, leg_mu, 0.3);
727
728 C_mu1->SaveAs(muEff+".eps");
729 C_mu->SaveAs(muRes+".eps");
730
731 /////////////
732 // Photons //
733 /////////////
734
735 // Reconstruction efficiency
736 int phID = 22;
737 std::pair<TH1D*,TH1D*> histos_ph = GetEff<Electron>(branchPhoton, branchParticle, "Photon", phID, treeReader);
738 std::pair<TH1D*,TH1D*> histos_phtower = GetEff<Electron>(branchTower, branchParticle, "Photon", phID, treeReader);
739
740 // Photon Energy Resolution
741 std::vector<resolPlot> plots_ph;
742 HistogramsCollection(&plots_ph, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "photons");
743 GetEres<Photon>( &plots_ph, branchPhoton, branchParticle, phID, treeReader);
744 TGraphErrors gr_ph = EresGraph(&plots_ph, true);
745 gr_ph.SetName("Photon");
746
747 // Photon Tower Energy Resolution
748 std::vector<resolPlot> plots_phtower;
749 HistogramsCollection(&plots_phtower, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "photonsTower");
750 GetEres<Tower>( &plots_phtower, branchTower, branchParticle, phID, treeReader);
751 TGraphErrors gr_phtower = EresGraph(&plots_phtower, true);
752 gr_phtower.SetName("PhotonTower");
753
754 // Canvas
755
756 TString phEff = "photonEff";
757 TCanvas *C_ph1 = new TCanvas(phEff,phEff, 1000, 500);
758 C_ph1->Divide(2);
759 C_ph1->cd(1);
760 TLegend *leg_ph1 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
761
762 gPad->SetLogx();
763 histos_phtower.first->Draw();
764 addHist(histos_phtower.first, leg_ph1, 3);
765 histos_ph.first->Draw("same");
766 addHist(histos_ph.first, leg_ph1, 1);
767
768 DrawAxis(histos_phtower.first, leg_ph1, 0);
769 leg_ph1->Draw();
770
771 C_ph1->cd(2);
772 TLegend *leg_ph2 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax);
773
774 histos_phtower.second->Draw("same");
775 addHist(histos_phtower.second, leg_ph2, 3);
776 histos_ph.second->Draw("same");
777 addHist(histos_ph.second, leg_ph2, 1);
778
779 DrawAxis(histos_phtower.second, leg_ph2, 1);
780
781 C_ph1->SaveAs(phEff+".eps");
782
783 TString phRes = "phERes";
784 TCanvas *C_ph = new TCanvas(phRes,phRes, 1000, 500);
785 TMultiGraph *mg_ph = new TMultiGraph(phRes,phRes);
786 TLegend *leg_ph = new TLegend(0.52,0.7,0.9,0.9);
787
788 addGraph(mg_ph, &gr_phtower, leg_ph, 3);
789 addGraph(mg_ph, &gr_ph, leg_ph, 1);
790
791 mg_ph->Draw("ACX");
792 leg_ph->Draw();
793
794 DrawAxis(mg_ph, leg_ph, 0.3);
795
796 C_ph->SaveAs(phRes+".eps");
797
798*/
799 //////////
800 // Jets //
801 //////////
802
803
804 // PFJets Energy Resolution
805 std::vector<resolPlot> plots_pfjets;
806 HistogramsCollection(&plots_pfjets, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "PFJet");
807 GetJetsEres( &plots_pfjets, branchPFJet, branchGenJet, treeReader);
808 TGraphErrors gr_pfjets = EresGraph(&plots_pfjets, true);
809 gr_pfjets.SetName("pfJet");
810
811
812 // PFJets Energy Resolution
813 std::vector<resolPlot> plots_calojets;
814 HistogramsCollection(&plots_calojets, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "CaloJet");
815 GetJetsEres( &plots_calojets, branchCaloJet, branchGenJet, treeReader);
816 TGraphErrors gr_calojets = EresGraph(&plots_calojets, true);
817 gr_calojets.SetName("caloJet");
818
819
820 TString jetRes = "jetERes";
821 TCanvas *C_jet = new TCanvas(jetRes,jetRes, 1000, 500);
822 TMultiGraph *mg_jet = new TMultiGraph(jetRes,jetRes);
823 TLegend *leg_jet = new TLegend(0.52,0.7,0.9,0.9);
824
825 addGraph(mg_jet, &gr_calojets, leg_jet, 3);
826 addGraph(mg_jet, &gr_pfjets, leg_jet, 1);
827
828 mg_jet->Draw("ACX");
829 leg_jet->Draw();
830
831 DrawAxis(mg_jet, leg_jet, 0.25);
832
833 C_jet->SaveAs(jetRes+".eps");
834
835
836 /*
837 // CaloJets Energy Resolution
838 std::vector<resolPlot> plots_calojets;
839 HistogramsCollection(&plots_calojets, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "caloJet");
840 GetJetsEres( &plots_calojets, branchCaloJet, branchGenJet, treeReader);
841 TGraphErrors gr_calojets = EresGraph(&plots_calojets, true);
842 gr_calojets.SetName("caloJet");
843 */
844
845
846
847 TFile *fout = new TFile(outputFile,"recreate");
848
849 for (int bin = 0; bin < Nbins; bin++)
850 {
851 plots_pfjets.at(bin).cenResolHist->Write();
852 plots_calojets.at(bin).cenResolHist->Write();
853 plots_el.at(bin).cenResolHist->Write();
854 plots_eltrack.at(bin).cenResolHist->Write();
855 plots_eltower.at(bin).cenResolHist->Write();
856 }
857
858/*
859 // gr.Write();
860 histos_el.first->Write();
861 //histos_el.second->Write();
862 histos_eltrack.first->Write();
863 //histos_eltrack.second->Write();
864 histos_eltower.first->Write();
865
866 histos_mu.first->Write();
867 histos_mu.second->Write();
868 histos_mutrack.first->Write();
869 histos_mutrack.second->Write();
870
871 histos_ph.first->Write();
872 histos_ph.second->Write();
873
874 //gr_el.Write();
875 //gr_eltrack.Write();
876 //gr_eltower.Write();
877*/
878
879 C_el1->Write();
880 C_el2->Write();
881 C_jet->Write();
882/*
883 C_mu->Write();
884 C_ph->Write();
885*/
886 gr_pfjets.Write();
887 gr_calojets.Write();
888
889 fout->Write();
890
891
892 cout << "** Exiting..." << endl;
893
894 //delete plots;
895 //delete result;
896 delete treeReader;
897 delete chain;
898}
899
900//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.