SchieferD/jetcalib/mcrecoeff_x.cpp

////////////////////////////////////////////////////////////////////////////////
//
// mcrecoeff_x
// -----------
//
//            06/16/2007 Philipp Schieferdecker <philipp.schieferdecker@cern.ch>
////////////////////////////////////////////////////////////////////////////////


#include "utils/cmdline.h"

#include <TROOT.h>
#include <TStyle.h>
#include <TRint.h>
#include <TFile.h>
#include <TTree.h>
#include <TH1F.h>
#include <TCanvas.h>
#include <TLegend.h>

#include <iostream>
#include <iomanip>
#include <cmath>
#include <string>
#include <vector>


using namespace std;


////////////////////////////////////////////////////////////////////////////////
// declare global functions
////////////////////////////////////////////////////////////////////////////////

//______________________________________________________________________________
void   setBinomialErrors(TH1F* hEff,const TH1F* hRec, const TH1F* hAll);
string getLegLabel(const string& jetalg);


////////////////////////////////////////////////////////////////////////////////
// main
////////////////////////////////////////////////////////////////////////////////

//______________________________________________________________________________
int main(int argc,char**argv)
{
  cmdline cl;
  cl.parse(argc,argv);

  vector<string> input    = cl.get_vector<string>("input");
  
  if (!cl.check()) return 0;
  cl.print();

  
  // jetalgs
  vector<string> jetalgs;
  for (vector<string>::const_iterator it=input.begin();it!=input.end();++it) {
    string jetalg = *it;
    string::size_type pos;
    while ((pos=jetalg.find('/'))!=string::npos) jetalg = jetalg.substr(pos+1);
    jetalg = jetalg.substr(0,jetalg.find("."));
    jetalgs.push_back(jetalg);
  }
  
  // loop over all input samples
  vector<TH1F*> etHistos;
  vector<TH1F*> ptHistos;
  vector<TH1F*> etaHistos;
  vector<TH1F*> phiHistos;
  
  for (unsigned int i=0;i<input.size();++i) {
    TFile* f = new TFile(input[i].c_str(),"READ");
    string hname;
    
    hname = "PtEff_"+jetalgs[i];
    TH1F* hPtAll = (TH1F*)f->Get("GenJetPtAll");
    TH1F* hPtRec = (TH1F*)f->Get("GenJetPtRec");
    TH1F* hPtEff = new TH1F(*hPtAll);
    hPtEff->SetNameTitle(hname.c_str(),hname.c_str());
    hPtEff->Divide(hPtRec,hPtAll,1.0,1.0,"B");
    setBinomialErrors(hPtEff,hPtRec,hPtAll);
    ptHistos.push_back(hPtEff);

    hname = "EtaEff_"+jetalgs[i];
    TH1F* hEtaAll = (TH1F*)f->Get("GenJetEtaAll");
    TH1F* hEtaRec = (TH1F*)f->Get("GenJetEtaRec");
    TH1F* hEtaEff = new TH1F(*hEtaAll);
    hEtaEff->SetNameTitle(hname.c_str(),hname.c_str());
    hEtaEff->Divide(hEtaRec,hEtaAll,1.0,1.0,"B");
    setBinomialErrors(hEtaEff,hEtaRec,hEtaAll);
    etaHistos.push_back(hEtaEff);

    hname = "PhiEff_"+jetalgs[i];
    TH1F* hPhiAll = (TH1F*)f->Get("GenJetPhiAll");
    TH1F* hPhiRec = (TH1F*)f->Get("GenJetPhiRec");
    TH1F* hPhiEff = new TH1F(*hPhiAll);
    hPhiEff->SetNameTitle(hname.c_str(),hname.c_str());
    hPhiEff->Divide(hPhiRec,hPhiAll,1.0,1.0,"B");
    setBinomialErrors(hPhiEff,hPhiRec,hPhiAll);
    phiHistos.push_back(hPhiEff);
  }
  
  
  // make plots
  argc=1;
  TRint* app = new TRint(argv[0],&argc,argv);
  
  gStyle->SetOptStat(0);
  
  Color_t color;
  
  // phi canvas
  TCanvas* cPhiEff  = new TCanvas("EffVsPhi", "EffVsPhi",510,400,500,400);
  TLegend* lPhiEff  = new TLegend(0.65,0.85-phiHistos.size()*0.05,0.86,0.85);
  cPhiEff->cd();
  lPhiEff->SetFillColor(10);
  
  color = kBlack;
  for (unsigned int i=0;i<phiHistos.size();i++) {
    TH1F* h = phiHistos[i];
    lPhiEff->AddEntry(h,getLegLabel(jetalgs[i]).c_str(),"l");
    h->SetMarkerColor(color);
    h->SetLineColor(color++);
    h->SetMarkerStyle(kFullCircle);
    h->SetMarkerSize(0.5);
    if (i==0) {
      h->SetTitle("#epsilon_{RECO} vs #phi");
      h->SetXTitle("genjet #phi");
      h->SetYTitle("#epsilon_{RECO}");
      h->Draw("E");
      h->SetMinimum(0.0);
      h->SetMaximum(1.05);
    }
    else h->Draw("ESAME");
  }
  lPhiEff->Draw();

  // eta canvas
  TCanvas* cEtaEff  = new TCanvas("EffVsEta", "EffVsEta",10,400,500,400);
  TLegend* lEtaEff  = new TLegend(0.65,0.85-etaHistos.size()*0.05,0.86,0.85);
  cEtaEff->cd();
  lEtaEff->SetFillColor(10);
  
  color = kBlack;
  for (unsigned int i=0;i<etaHistos.size();i++) {
    TH1F* h = etaHistos[i];
    lEtaEff->AddEntry(h,getLegLabel(jetalgs[i]).c_str(),"l");
    h->SetMarkerColor(color);
    h->SetLineColor(color++);
    h->SetMarkerStyle(kFullCircle);
    h->SetMarkerSize(0.5);
    if (i==0) {
      h->SetTitle("#epsilon_{RECO} vs #eta");
      h->SetXTitle("genjet #eta");
      h->SetYTitle("#epsilon_{RECO}");
      h->Draw("E");
      h->SetMinimum(0.0);
      h->SetMaximum(1.05);
    }
    else h->Draw("ESAME");
    
  }
  lEtaEff->Draw();

  // pt canvas
  TCanvas* cPtEff  = new TCanvas("EffVsPt", "EffVsPt",10,10,500,400);
  TLegend* lPtEff  = new TLegend(0.65,0.22+ptHistos.size()*0.05,0.86,0.22);
  cPtEff->cd();
  lPtEff->SetFillColor(10);
  
  color = kBlack;
  for (unsigned int i=0;i<ptHistos.size();i++) {
    TH1F* h = ptHistos[i];
    lPtEff->AddEntry(h,getLegLabel(jetalgs[i]).c_str(),"l");
    h->SetMarkerColor(color);
    h->SetLineColor(color++);
    h->SetMarkerStyle(kFullCircle);
    h->SetMarkerSize(0.5);
    if (i==0) {
      h->SetTitle("#epsilon_{RECO} vs p_{T}");
      h->SetXTitle("genjet p_{T} [GeV]");
      h->SetYTitle("#epsilon_{RECO}");
      h->Draw("E");
      h->SetMinimum(0.0);
      h->SetMaximum(1.05);
    }
    else h->Draw("ESAME");
  }
  lPtEff->Draw();

  app->Run();
  
  return 0;
}


////////////////////////////////////////////////////////////////////////////////
// implementation of global function
////////////////////////////////////////////////////////////////////////////////

//______________________________________________________________________________
void setBinomialErrors(TH1F* hEff,const TH1F* hRec, const TH1F* hAll)
{
  for (int i=1;i<hEff->GetNbinsX();i++) {
    double nrec = hRec->GetBinContent(i);
    double nall = hAll->GetBinContent(i);
    double eeff = (nall>0.0) ? std::sqrt(nrec/(nall*nall)*(1-nrec/nall)) : 0.0;
    hEff->SetBinError(i,eeff);
  }
}


//______________________________________________________________________________
string getLegLabel(const string& jetalg)
{
  string result;
  string tmp(jetalg);
  
  if (tmp.find("kt")==0) {
    result = "fast k_{T}, D=";
    tmp=tmp.substr(2);
  }
  else if (tmp.find("sc")==0) {
    result = "SISCone, R=";
    tmp = tmp.substr(2);
  }
  else if (tmp.find("ic")==0) {
    result = "Iterative Cone, R=";
    tmp = tmp.substr(2);
  }
  else if (tmp.find("mc")==0) {
    result = "Midpoint, R=";
    tmp = tmp.substr(2);
  }
  else return result;

  stringstream ssparam;
  ssparam<<tmp;
  float param;
  ssparam>>param;
  param = 0.1*param;
  stringstream ssparam2;
  ssparam2<<param;
  result += ssparam2.str();

  return result;
}
Revision:	1.4
Committed:	Fri Nov 23 13:34:44 2007 UTC (17 years, 5 months ago) by schiefer
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.3:	+0 -36 lines
Error occurred while calculating annotation data.
Log Message:	* empty log message *
#	Content
1	////////////////////////////////////////////////////////////////////////////////
2	//
3	// mcrecoeff_x
4	// -----------
5	//
6	// 06/16/2007 Philipp Schieferdecker <philipp.schieferdecker@cern.ch>
7	////////////////////////////////////////////////////////////////////////////////
8
9
10	#include "utils/cmdline.h"
11
12	#include <TROOT.h>
13	#include <TStyle.h>
14	#include <TRint.h>
15	#include <TFile.h>
16	#include <TTree.h>
17	#include <TH1F.h>
18	#include <TCanvas.h>
19	#include <TLegend.h>
20
21	#include <iostream>
22	#include <iomanip>
23	#include <cmath>
24	#include <string>
25	#include <vector>
26
27
28	using namespace std;
29
30
31	////////////////////////////////////////////////////////////////////////////////
32	// declare global functions
33	////////////////////////////////////////////////////////////////////////////////
34
35	//______________________________________________________________________________
36	void setBinomialErrors(TH1F* hEff,const TH1F* hRec, const TH1F* hAll);
37	string getLegLabel(const string& jetalg);
38
39
40	////////////////////////////////////////////////////////////////////////////////
41	// main
42	////////////////////////////////////////////////////////////////////////////////
43
44	//______________________________________________________________________________
45	int main(int argc,char**argv)
46	{
47	cmdline cl;
48	cl.parse(argc,argv);
49
50	vector<string> input = cl.get_vector<string>("input");
51
52	if (!cl.check()) return 0;
53	cl.print();
54
55
56	// jetalgs
57	vector<string> jetalgs;
58	for (vector<string>::const_iterator it=input.begin();it!=input.end();++it) {
59	string jetalg = *it;
60	string::size_type pos;
61	while ((pos=jetalg.find('/'))!=string::npos) jetalg = jetalg.substr(pos+1);
62	jetalg = jetalg.substr(0,jetalg.find("."));
63	jetalgs.push_back(jetalg);
64	}
65
66	// loop over all input samples
67	vector<TH1F*> etHistos;
68	vector<TH1F*> ptHistos;
69	vector<TH1F*> etaHistos;
70	vector<TH1F*> phiHistos;
71
72	for (unsigned int i=0;i<input.size();++i) {
73	TFile* f = new TFile(input[i].c_str(),"READ");
74	string hname;
75
76	hname = "PtEff_"+jetalgs[i];
77	TH1F* hPtAll = (TH1F*)f->Get("GenJetPtAll");
78	TH1F* hPtRec = (TH1F*)f->Get("GenJetPtRec");
79	TH1F* hPtEff = new TH1F(*hPtAll);
80	hPtEff->SetNameTitle(hname.c_str(),hname.c_str());
81	hPtEff->Divide(hPtRec,hPtAll,1.0,1.0,"B");
82	setBinomialErrors(hPtEff,hPtRec,hPtAll);
83	ptHistos.push_back(hPtEff);
84
85	hname = "EtaEff_"+jetalgs[i];
86	TH1F* hEtaAll = (TH1F*)f->Get("GenJetEtaAll");
87	TH1F* hEtaRec = (TH1F*)f->Get("GenJetEtaRec");
88	TH1F* hEtaEff = new TH1F(*hEtaAll);
89	hEtaEff->SetNameTitle(hname.c_str(),hname.c_str());
90	hEtaEff->Divide(hEtaRec,hEtaAll,1.0,1.0,"B");
91	setBinomialErrors(hEtaEff,hEtaRec,hEtaAll);
92	etaHistos.push_back(hEtaEff);
93
94	hname = "PhiEff_"+jetalgs[i];
95	TH1F* hPhiAll = (TH1F*)f->Get("GenJetPhiAll");
96	TH1F* hPhiRec = (TH1F*)f->Get("GenJetPhiRec");
97	TH1F* hPhiEff = new TH1F(*hPhiAll);
98	hPhiEff->SetNameTitle(hname.c_str(),hname.c_str());
99	hPhiEff->Divide(hPhiRec,hPhiAll,1.0,1.0,"B");
100	setBinomialErrors(hPhiEff,hPhiRec,hPhiAll);
101	phiHistos.push_back(hPhiEff);
102	}
103
104
105	// make plots
106	argc=1;
107	TRint* app = new TRint(argv[0],&argc,argv);
108
109	gStyle->SetOptStat(0);
110
111	Color_t color;
112
113	// phi canvas
114	TCanvas* cPhiEff = new TCanvas("EffVsPhi", "EffVsPhi",510,400,500,400);
115	TLegend* lPhiEff = new TLegend(0.65,0.85-phiHistos.size()*0.05,0.86,0.85);
116	cPhiEff->cd();
117	lPhiEff->SetFillColor(10);
118
119	color = kBlack;
120	for (unsigned int i=0;i<phiHistos.size();i++) {
121	TH1F* h = phiHistos[i];
122	lPhiEff->AddEntry(h,getLegLabel(jetalgs[i]).c_str(),"l");
123	h->SetMarkerColor(color);
124	h->SetLineColor(color++);
125	h->SetMarkerStyle(kFullCircle);
126	h->SetMarkerSize(0.5);
127	if (i==0) {
128	h->SetTitle("#epsilon_{RECO} vs #phi");
129	h->SetXTitle("genjet #phi");
130	h->SetYTitle("#epsilon_{RECO}");
131	h->Draw("E");
132	h->SetMinimum(0.0);
133	h->SetMaximum(1.05);
134	}
135	else h->Draw("ESAME");
136	}
137	lPhiEff->Draw();
138
139	// eta canvas
140	TCanvas* cEtaEff = new TCanvas("EffVsEta", "EffVsEta",10,400,500,400);
141	TLegend* lEtaEff = new TLegend(0.65,0.85-etaHistos.size()*0.05,0.86,0.85);
142	cEtaEff->cd();
143	lEtaEff->SetFillColor(10);
144
145	color = kBlack;
146	for (unsigned int i=0;i<etaHistos.size();i++) {
147	TH1F* h = etaHistos[i];
148	lEtaEff->AddEntry(h,getLegLabel(jetalgs[i]).c_str(),"l");
149	h->SetMarkerColor(color);
150	h->SetLineColor(color++);
151	h->SetMarkerStyle(kFullCircle);
152	h->SetMarkerSize(0.5);
153	if (i==0) {
154	h->SetTitle("#epsilon_{RECO} vs #eta");
155	h->SetXTitle("genjet #eta");
156	h->SetYTitle("#epsilon_{RECO}");
157	h->Draw("E");
158	h->SetMinimum(0.0);
159	h->SetMaximum(1.05);
160	}
161	else h->Draw("ESAME");
162
163	}
164	lEtaEff->Draw();
165
166	// pt canvas
167	TCanvas* cPtEff = new TCanvas("EffVsPt", "EffVsPt",10,10,500,400);
168	TLegend* lPtEff = new TLegend(0.65,0.22+ptHistos.size()*0.05,0.86,0.22);
169	cPtEff->cd();
170	lPtEff->SetFillColor(10);
171
172	color = kBlack;
173	for (unsigned int i=0;i<ptHistos.size();i++) {
174	TH1F* h = ptHistos[i];
175	lPtEff->AddEntry(h,getLegLabel(jetalgs[i]).c_str(),"l");
176	h->SetMarkerColor(color);
177	h->SetLineColor(color++);
178	h->SetMarkerStyle(kFullCircle);
179	h->SetMarkerSize(0.5);
180	if (i==0) {
181	h->SetTitle("#epsilon_{RECO} vs p_{T}");
182	h->SetXTitle("genjet p_{T} [GeV]");
183	h->SetYTitle("#epsilon_{RECO}");
184	h->Draw("E");
185	h->SetMinimum(0.0);
186	h->SetMaximum(1.05);
187	}
188	else h->Draw("ESAME");
189	}
190	lPtEff->Draw();
191
192	app->Run();
193
194	return 0;
195	}
196
197
198	////////////////////////////////////////////////////////////////////////////////
199	// implementation of global function
200	////////////////////////////////////////////////////////////////////////////////
201
202	//______________________________________________________________________________
203	void setBinomialErrors(TH1F* hEff,const TH1F* hRec, const TH1F* hAll)
204	{
205	for (int i=1;i<hEff->GetNbinsX();i++) {
206	double nrec = hRec->GetBinContent(i);
207	double nall = hAll->GetBinContent(i);
208	double eeff = (nall>0.0) ? std::sqrt(nrec/(nallnall)(1-nrec/nall)) : 0.0;
209	hEff->SetBinError(i,eeff);
210	}
211	}
212
213
214	//______________________________________________________________________________
215	string getLegLabel(const string& jetalg)
216	{
217	string result;
218	string tmp(jetalg);
219
220	if (tmp.find("kt")==0) {
221	result = "fast k_{T}, D=";
222	tmp=tmp.substr(2);
223	}
224	else if (tmp.find("sc")==0) {
225	result = "SISCone, R=";
226	tmp = tmp.substr(2);
227	}
228	else if (tmp.find("ic")==0) {
229	result = "Iterative Cone, R=";
230	tmp = tmp.substr(2);
231	}
232	else if (tmp.find("mc")==0) {
233	result = "Midpoint, R=";
234	tmp = tmp.substr(2);
235	}
236	else return result;
237
238	stringstream ssparam;
239	ssparam<<tmp;
240	float param;
241	ssparam>>param;
242	param = 0.1*param;
243	stringstream ssparam2;
244	ssparam2<<param;
245	result += ssparam2.str();
246
247	return result;
248	}