ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/SchieferD/jetcalib/mcresponse_x.cpp

Comparing UserCode/SchieferD/jetcalib/mcresponse_x.cpp (file contents):
Revision 1.2 by schiefer, Tue Aug 21 06:42:45 2007 UTC vs.
Revision 1.7 by schiefer, Sat Sep 1 12:16:09 2007 UTC

#	Line 17 | Line 17
17		#include <TEventList.h>
18		#include <TH1F.h>
19		#include <TF1.h>
20	<	#include <TGraphErrors.h>
21	<	#include <TMultiGraph.h>
22	<	#include <TCanvas.h>
23	<	#include <TLegend.h>
20	>	#include <TMath.h>
21
22		#include <iostream>
23		#include <iomanip>
27	–	#include <fstream>
24		#include <cmath>
25		#include <string>
26		#include <vector>
#	Line 34 | Line 30 | using namespace std;
30
31
32		////////////////////////////////////////////////////////////////////////////////
37	–	// defaults
38	–	////////////////////////////////////////////////////////////////////////////////
39	–	string dflt_ptbins ="7,9,12,15,20,25,30,45,65,90,120";
40	–	string dflt_etabins=".2,.4,.7,1.,1.3,1.5,1.75,2.,2.3,2.5,2.7,2.9,3.2";
41	–	//string dflt_etabins=".3,.6,.9,1.2,1.5,1.8,2.1,2.4,2.7";
42	–	string dflt_phibins="-2.75,-2.25,-1.75.,-1.25,-.75,-.25,.25,.75,1.25,1.75,2.25,2.75";
43	–
44	–
45	–	////////////////////////////////////////////////////////////////////////////////
33		// declare global functions
34		////////////////////////////////////////////////////////////////////////////////
35	+	void     replaceHistos(int nbins,
36	+	TH1F**& hVar,TH1F**& hEtCorr,TH1F**& hAbsRsp,TH1F**& hRelRsp,
37	+	TTree**& tVar);
38	+	void     replaceHistos(int nbins1,int nbins2,
39	+	TH1F***&hVar,TH1F***&hEtCorr,TH1F***&hAbsRsp,TH1F***&hRelRsp,
40	+	TTree***& tVar);
41	+	void     replaceHistos(int nbins1,int nbins2,int nbins3,
42	+	TH1F****&hVar,TH1F****&hEtCorr,
43	+	TH1F****&hAbsRsp,TH1F****&hRelRsp,
44	+	TTree****& tVar);
45	+	void     fitHistos    (TH1F**    histos,int n,double nsigma);
46	+	void     fitHistos    (TH1F***   histos,int n1,int n2,double nsigma);
47	+	void     fitHistos    (TH1F****  histos,int n1,int n2,int n3,double nsigma);
48	+	void     saveHistos   (TH1F**&   histos,int n);
49	+	void     saveHistos   (TH1F***&  histos,int n1,int n2);
50	+	void     saveHistos   (TH1F****& histos,int n1,int n2,int n3);
51
49	–	void     fillRspAndSgm(TGraphErrors* gRsp,TGraphErrors* gSgm,float nsigma,
50	–	const vector<float>& bins,TH1F** hVar,TH1F** hPt,TH1F** hRsp);
51	–	TCanvas* plotVarHistos(const vector<float>& bins,TH1F** hVar,bool logy=false);
52	–	TCanvas* plotRspHistos(const vector<float>& bins,TH1F** hRsp);
53	–	TCanvas* plotResponse (const vector<string>& jetalgs,short applyjes,
54	–	const vector<TGraphErrors*> gRsp,
55	–	const string& varname,const string& xtitle);
56	–	TCanvas* plotResponse (const string& jetalg,short applyjes,bool fitcorr,
57	–	const vector<float>& bins1,TGraphErrors** gRsp,
58	–	const string& varname1,const string vartitle1,
59	–	const string& varname2,const string& xtitle);
60	–	TCanvas* plotResolution(const vector<string>& jetalgs,
61	–	const vector<TGraphErrors*> gSgm,
62	–	const string& varname,const string& xtitle);
63	–	TCanvas* plotResolution(const string& jetalg,const vector<float>& bins1,
64	–	TGraphErrors** gSgm,
65	–	const string& varname1,const string vartitle1,
66	–	const string& varname2,const string& xtitle);
67	–
68	–	float  get_xmax(TGraph* g);
52
53
54		////////////////////////////////////////////////////////////////////////////////
#	Line 79 | Line 62 | int main(int argc,char**argv)
62		cl.parse(argc,argv);
63
64		vector<string> input     = cl.get_vector<string>("input");
65	+	string         output    = cl.get_value<string> ("output", "mcrsp.root");
66		string         selection = cl.get_value<string> ("selection",   "njt>0");
67		float          drmax     = cl.get_value<float>  ("drmax",           0.3);
68	+	short          applyjes  = cl.get_value<short>  ("applyjes",          0);
69		string         treename  = cl.get_value<string> ("treename",        "t");
70		float          nsigma    = cl.get_value<float>  ("nsigma",          3.0);
71	<	short          applyjes  = cl.get_value<short>  ("applyjes",          0);
72	<	bool           fitcorr   = cl.get_value<bool>   ("fitcorr",       false);
73	<	vector<float>  ptbins    = cl.get_vector<float> ("ptbins",  dflt_ptbins);
74	<	vector<float>  etabins   = cl.get_vector<float> ("etabins",dflt_etabins);
90	<	vector<float>  phibins   = cl.get_vector<float> ("phibins",dflt_phibins);
71	>	vector<float>  etbins    = cl.get_vector<float> ("etbins",           "");
72	>	vector<float>  etabins   = cl.get_vector<float> ("etabins",          "");
73	>	vector<float>  phibins   = cl.get_vector<float> ("phibins",          "");
74	>	vector<float>  emfbins   = cl.get_vector<float> ("emfbins",          "");
75		bool           abseta    = cl.get_value<bool>   ("abseta",         true);
76
77		if (!cl.check()) return 0;
78		cl.print();
79
80	<
80	>	int netbins  = (etbins.size() >0) ? etbins.size() +1 : 0;
81	>	int netabins = (etabins.size()>0) ? etabins.size()+1 : 0;
82	>	int nphibins = (phibins.size()>0) ? phibins.size()+1 : 0;
83	>	int nemfbins = (emfbins.size()>0) ? emfbins.size()+1 : 0;
84	>
85	>	if (netbins+netabins+nphibins+nemfbins==0) {
86	>	cout<<"Must specify bins: etbins, etabins, phibins, AND/OR emfbins!"<<endl;
87	>	return 0;
88	>	}
89	>
90		// etabins
91	<	if (!abseta) {
91	>	if (netabins>0&&!abseta) {
92		int neta=(int)etabins.size();
93		std::reverse(etabins.begin(),etabins.end());
94		for (int ieta=neta-1;ieta>=0;ieta--) etabins.push_back(etabins[ieta]);
#	Line 116 | Line 109 | int main(int argc,char**argv)
109		// prepare branch variables and histograms / graphs
110		float weight;
111		char  njt;
112	<	float jtpt[100];
112	>	float jtet[100];
113		float jteta[100];
114		float jtphi[100];
115	<	float jtgenpt[100];
115	>	float jtemf[100];
116	>	float jtgenet[100];
117		float jtgendr[100];
118		float jtjes[100][3];
119
120	<	int nptbins  = ptbins.size()+1;
121	<	int netabins = etabins.size()+1;
128	<	int nphibins = phibins.size()+1;
129	<
130	<	vector<TH1F**>         hGenPt;
131	<	vector<TH1F**>         hGenPtPtCorr;
132	<	vector<TH1F**>         hRspVsGenPt;
133	<	vector<TH1F**>         hPt;
134	<	vector<TH1F**>         hPtCorr;
135	<	vector<TH1F**>         hRspVsPt;
136	<	vector<TH1F**>         hEta;
137	<	vector<TH1F**>         hEtaPtCorr;
138	<	vector<TH1F**>         hRspVsEta;
139	<	vector<TH1F**>         hPhi;
140	<	vector<TH1F**>         hPhiPtCorr;
141	<	vector<TH1F**>         hRspVsPhi;
142	<
143	<	vector<TH1F***>        hEtaPtPt;
144	<	vector<TH1F***>        hEtaPtPtCorr;
145	<	vector<TH1F***>        hRspVsEtaPt;
146	<
147	<	vector<TGraphErrors*>  gRspVsGenPt;
148	<	vector<TGraphErrors*>  gSgmVsGenPt;
149	<	vector<TGraphErrors*>  gRspVsPt;
150	<	vector<TGraphErrors*>  gSgmVsPt;
151	<	vector<TGraphErrors*>  gRspVsEta;
152	<	vector<TGraphErrors*>  gSgmVsEta;
153	<	vector<TGraphErrors*>  gRspVsPhi;
154	<	vector<TGraphErrors*>  gSgmVsPhi;
155	<
156	<	vector<TGraphErrors**> gRspVsEtaPt;
157	<	vector<TGraphErrors**> gSgmVsEtaPt;
158	<
159	<	argc=1; //argv[1]="-b";
160	<	TRint* app = new TRint(argv[0],&argc,argv);
120	>	argc=3; argv[1]="-l"; argv[2]="-b";
121	>	TRint* app = new TRint(argv[0],&argc,argv); app->Argc();
122
123
124		//
125		// loop over all input files
126		//
127	<	TFile* plotfile = new TFile("mcresponse.root","RECREATE");
127	>	TFile* plotfile = new TFile(output.c_str(),"UPDATE");
128
129		for (unsigned int i=0;i<input.size();++i) {
130	<	TFile* f = new TFile(input[i].c_str(),"READ"); if (!f->IsOpen()) return 0;
130	>	TFile* f = new TFile(input[i].c_str(),"UPDATE"); if (!f->IsOpen()) return 0;
131		TTree* t = (TTree*)f->Get(treename.c_str());
132
133		TEventList* el = new TEventList("sel","sel");
#	Line 174 | Line 135 | int main(int argc,char**argv)
135
136		t->SetBranchAddress("weight", &weight);
137		t->SetBranchAddress("njt",    &njt);
138	<	t->SetBranchAddress("jtpt",   jtpt);
138	>	t->SetBranchAddress("jtet",   jtet);
139		t->SetBranchAddress("jteta",  jteta);
140		t->SetBranchAddress("jtphi",  jtphi);
141	<	t->SetBranchAddress("jtgenpt",jtgenpt);
141	>	t->SetBranchAddress("jtemf",  jtemf);
142	>	t->SetBranchAddress("jtgenet",jtgenet);
143		t->SetBranchAddress("jtgendr",jtgendr);
144
145		if (t->FindBranch("jtjes")) t->SetBranchAddress("jtjes", jtjes);
#	Line 185 | Line 147 | int main(int argc,char**argv)
147		for (int i1=0;i1<100;i1++)
148		for (int i2=0;i2<3;i2++) jtjes[i1][i2]=1.0;
149
150	+	string absrsp_xtitle = "E_{T}^{gen}-E_{T} [GeV]";
151	+	string relrsp_xtitle = "E_{T}^{rec}/E_{T}";
152	+
153
154	<	string rsp_xtitle = "E_{T}^{gen}-E_{T}^{rec} [GeV]";
154	>	TH1F** hGenEt(0);
155	>	TH1F** hEtCorrGenEt(0);
156	>	TH1F** hAbsRspGenEt(0);
157	>	TH1F** hRelRspGenEt(0);
158	>
159	>	if (netbins>0) {
160	>	hGenEt       = hist::initHistos("GenEt",jetalgs[i],100,
161	>	"jet E_{T}^{gen} [GeV]",
162	>	"GenEt",etbins);
163	>	hEtCorrGenEt = hist::initHistos("EtCorrGenEt",jetalgs[i],100,
164	>	"jet E_{T}^{gen} [GeV]",
165	>	"GenEt",etbins);
166	>	hAbsRspGenEt = hist::initHistos("AbsRspGenEt",jetalgs[i],
167	>	100,-50,150,
168	>	absrsp_xtitle,
169	>	"GenEt",etbins);
170	>	hRelRspGenEt = hist::initHistos("RelRspGenEt",jetalgs[i],
171	>	100,0,2,
172	>	relrsp_xtitle,
173	>	"GenEt",etbins);
174	>	}
175
176	<	hGenPt      .push_back(hist::initHistos("JetGenPt",jetalgs[i],
177	<	4000,0,4000,"jet p_{T}^{gen} [GeV]",
178	<	"GenPt",ptbins));
179	<	hGenPtPtCorr.push_back(hist::initHistos("JetPtCorr",jetalgs[i],
180	<	4000,0,4000,"jet p_{T}^{gen} [GeV]",
181	<	"GenPt",ptbins));
182	<	hRspVsGenPt .push_back(hist::initHistos("JetRspGenPt",jetalgs[i],
183	<	100,-100,100,rsp_xtitle,
184	<	"GenPt",ptbins));
185	<
186	<	hPt         .push_back(hist::initHistos("JetPt",jetalgs[i],
187	<	4000,0,4000,"jet p_{T} [GeV]",
188	<	"Pt",ptbins));
189	<	hPtCorr     .push_back(hist::initHistos("JetPtCorr",jetalgs[i],
190	<	4000,0,4000,"jet p_{T} [GeV]",
191	<	"Pt",ptbins));
207	<	hRspVsPt    .push_back(hist::initHistos("JetRspPt",jetalgs[i],
208	<	100,-100,100,rsp_xtitle,
209	<	"Pt",ptbins));
210	<
211	<	hEta        .push_back(hist::initHistos("JetEta",jetalgs[i],
212	<	100,0,5,"jet #eta","Eta",etabins));
213	<	hEtaPtCorr  .push_back(hist::initHistos("JetPtCorr",jetalgs[i],
214	<	4000,0,4000,"jet p_{T} [GeV]",
215	<	"Eta",etabins));
216	<	hRspVsEta   .push_back(hist::initHistos("JetRspEta",jetalgs[i],
217	<	100,-100,100,rsp_xtitle,"Eta",etabins));
218	<
219	<	hPhi        .push_back(hist::initHistos("JetPhi",jetalgs[i],
220	<	100,-3.2,3.2,"jet #phi","Phi",phibins));
221	<	hPhiPtCorr  .push_back(hist::initHistos("JetPtCorr",jetalgs[i],
222	<	4000,0,4000,"jet p_{T} [Gev]",
223	<	"Phi",phibins));
224	<	hRspVsPhi   .push_back(hist::initHistos("JetRspPhi",jetalgs[i],
225	<	100,-100,100,rsp_xtitle,"Phi",phibins));
226	<
227	<	hEtaPtPt    .push_back(hist::initHistos("JetPt",jetalgs[i],
228	<	4000,0,4000,"jet p_{T} [GeV]",
229	<	"Eta",etabins,"Pt",ptbins));
230	<	hEtaPtPtCorr.push_back(hist::initHistos("JetPtCorr",jetalgs[i],
231	<	4000,0,4000,"jet p_{T} [GeV]",
232	<	"Eta",etabins,"Pt",ptbins));
233	<	hRspVsEtaPt .push_back(hist::initHistos("JetRspEtaPt",jetalgs[i],
234	<	100,-100,100,rsp_xtitle,
235	<	"Eta",etabins,"Pt",ptbins));
176	>
177	>	TH1F** hEt(0);
178	>	TH1F** hEtCorrEt(0);
179	>	TH1F** hAbsRspEt(0);
180	>	TH1F** hRelRspEt(0);
181	>
182	>	if (netbins>0) {
183	>	hEt       = hist::initHistos("Et",jetalgs[i],100,
184	>	"jet E_{T} [GeV]","Et",etbins);
185	>	hEtCorrEt = hist::initHistos("EtCorrEt",jetalgs[i],100,
186	>	"jet E_{T} [GeV]","Et",etbins);
187	>	hAbsRspEt = hist::initHistos("AbsRspEt",jetalgs[i],100,-50,150,
188	>	absrsp_xtitle,"Et",etbins);
189	>	hRelRspEt = hist::initHistos("RelRspEt",jetalgs[i],100,0,2,
190	>	relrsp_xtitle,"Et",etbins);
191	>	}
192
193	+
194	+	TH1F** hEta(0);
195	+	TH1F** hEtCorrEta(0);
196	+	TH1F** hAbsRspEta(0);
197	+	TH1F** hRelRspEta(0);
198	+
199	+	if (netabins>0) {
200	+	hEta       = hist::initHistos("Eta",jetalgs[i],100,
201	+	"jet #eta","Eta",etabins);
202	+	hEtCorrEta = hist::initHistos("EtCorrEta",jetalgs[i],100,
203	+	"jet E_{T} [GeV]","Eta",etabins);
204	+	hAbsRspEta = hist::initHistos("AbsRspEta",jetalgs[i],100,-50,150,
205	+	absrsp_xtitle,"Eta",etabins);
206	+	hRelRspEta = hist::initHistos("RelRspEta",jetalgs[i],100,0,2,
207	+	relrsp_xtitle,"Eta",etabins);
208	+	}
209	+
210	+
211	+	TH1F** hPhi(0);
212	+	TH1F** hEtCorrPhi(0);
213	+	TH1F** hAbsRspPhi(0);
214	+	TH1F** hRelRspPhi(0);
215	+
216	+	if (nphibins>0) {
217	+	hPhi       = hist::initHistos("Phi",jetalgs[i],100,
218	+	"jet #phi","Phi",phibins);
219	+	hEtCorrPhi = hist::initHistos("EtCorrPhi",jetalgs[i],100,
220	+	"jet E_{T} [Gev]","Phi",phibins);
221	+	hAbsRspPhi = hist::initHistos("AbsRspPhi",jetalgs[i],100,-50,150,
222	+	absrsp_xtitle,"Phi",phibins);
223	+	hRelRspPhi = hist::initHistos("RelRspPhi",jetalgs[i],100,0,2,
224	+	relrsp_xtitle,"Phi",phibins);
225	+	}
226	+
227	+
228	+	TH1F** hEmf(0);
229	+	TH1F** hEtCorrEmf(0);
230	+	TH1F** hAbsRspEmf(0);
231	+	TH1F** hRelRspEmf(0);
232	+
233	+	if (nemfbins>0) {
234	+	hEmf       = hist::initHistos("Emf",jetalgs[i],100,
235	+	"jet emf","Emf",emfbins);
236	+	hEtCorrEmf = hist::initHistos("EtCorrEmf",jetalgs[i],100,
237	+	"jet E_{T} [Gev]","Emf",emfbins);
238	+	hAbsRspEmf = hist::initHistos("AbsRspEmf",jetalgs[i],100,-50,150,
239	+	absrsp_xtitle,"Emf",emfbins);
240	+	hRelRspEmf = hist::initHistos("RelRspEmf",jetalgs[i],100,0,2,
241	+	relrsp_xtitle,"Emf",emfbins);
242	+	}
243	+
244	+
245	+	TH1F*** hEtEtaEt(0);
246	+	TH1F*** hEtCorrEtaEt(0);
247	+	TH1F*** hAbsRspEtaEt(0);
248	+	TH1F*** hRelRspEtaEt(0);
249	+
250	+	if (netbins>0&&netabins>0) {
251	+	hEtEtaEt     = hist::initHistos("Et",jetalgs[i],100,"jet E_{T} [GeV]",
252	+	"Eta",etabins,"Et",etbins);
253	+	hEtCorrEtaEt = hist::initHistos("EtCorrEtaEt",jetalgs[i],100,
254	+	"jet E_{T} [GeV]",
255	+	"Eta",etabins,"Et",etbins);
256	+	hAbsRspEtaEt = hist::initHistos("AbsRspEtaEt",jetalgs[i],100,-50,150,
257	+	absrsp_xtitle,
258	+	"Eta",etabins,"Et",etbins);
259	+	hRelRspEtaEt = hist::initHistos("RelRspEtaEt",jetalgs[i],100,0,2,
260	+	relrsp_xtitle,
261	+	"Eta",etabins,"Et",etbins);
262	+	}
263	+
264	+
265	+	TH1F*** hEtEmfEt(0);
266	+	TH1F*** hEtCorrEmfEt(0);
267	+	TH1F*** hAbsRspEmfEt(0);
268	+	TH1F*** hRelRspEmfEt(0);
269	+
270	+	if (netbins>0&&nemfbins>0) {
271	+	hEtEmfEt     = hist::initHistos("Et",jetalgs[i],100,"jet E_{T} [GeV]",
272	+	"Emf",emfbins,"Et",etbins);
273	+	hEtCorrEmfEt = hist::initHistos("EtCorrEmfEt",jetalgs[i],100,
274	+	"jet E_{T} [GeV]",
275	+	"Emf",emfbins,"Et",etbins);
276	+	hAbsRspEmfEt = hist::initHistos("AbsRspEmfEt",jetalgs[i],100,-50,150,
277	+	absrsp_xtitle,
278	+	"Emf",emfbins,"Et",etbins);
279	+	hRelRspEmfEt = hist::initHistos("RelRspEmfEt",jetalgs[i],100,0,2,
280	+	relrsp_xtitle,
281	+	"Emf",emfbins,"Et",etbins);
282	+	}
283	+
284	+
285	+	TH1F**** hEtEmfEtaEt(0);
286	+	TH1F**** hEtCorrEmfEtaEt(0);
287	+	TH1F**** hAbsRspEmfEtaEt(0);
288	+	TH1F**** hRelRspEmfEtaEt(0);
289	+
290	+	if (netbins>0&&netabins>0&&nemfbins>0) {
291	+	hEtEmfEtaEt     = hist::initHistos("Et",jetalgs[i],100,"jet E_{T} [GeV]",
292	+	"Emf",emfbins,
293	+	"Eta",etabins,
294	+	"Et", etbins);
295	+	hEtCorrEmfEtaEt = hist::initHistos("EtCorrEmfEtaEt",jetalgs[i],100,
296	+	"jet E_{T} [GeV]",
297	+	"Emf",emfbins,
298	+	"Eta",etabins,
299	+	"Et", etbins);
300	+	hAbsRspEmfEtaEt = hist::initHistos("AbsRspEmfEtaEt",jetalgs[i],
301	+	100,-50,150,absrsp_xtitle,
302	+	"Emf",emfbins,
303	+	"Eta",etabins,
304	+	"Et", etbins);
305	+	hRelRspEmfEtaEt = hist::initHistos("RelRspEmfEtaEt",jetalgs[i],
306	+	100,0,2,relrsp_xtitle,
307	+	"Emf",emfbins,
308	+	"Eta",etabins,
309	+	"Et", etbins);
310	+	}
311
312		int nevts = el->GetN();
313		cout<<jetalgs[i]<<": "<<nevts<<" events selected."<<endl;
314
315	+	// prepare intermediate trees
316	+	float genet,et,eta,phi,emf,etcorr,absrsp,relrsp;
317	+
318	+	TTree**   tGenEt(0);
319	+	TTree**   tEt(0);
320	+	TTree**   tEta(0);
321	+	TTree**   tPhi(0);
322	+	TTree**   tEmf(0);
323	+	TTree***  tEtEtaEt(0);
324	+	TTree***  tEtEmfEt(0);
325	+	TTree**** tEtEmfEtaEt(0);
326	+
327	+	if (netbins>0)  tGenEt = new TTree*[netbins];
328	+	if (netbins>0)  tEt    = new TTree*[netbins];
329	+	if (netabins>0) tEta   = new TTree*[netabins];
330	+	if (nphibins>0) tPhi   = new TTree*[nphibins];
331	+	if (nemfbins>0) tEmf   = new TTree*[nemfbins];
332	+	if (netbins>0&&netabins>0) tEtEtaEt = new TTree**[netabins];
333	+	if (netbins>0&&nemfbins>0) tEtEmfEt = new TTree**[nemfbins];
334	+	if (netbins>0&&netabins>0&&nemfbins>0) tEtEmfEtaEt = new TTree***[nemfbins];
335	+
336	+	for (int iet=0;iet<netbins;iet++) {
337	+
338	+	tGenEt[iet]=new TTree("tGenEt","tGenEt");
339	+	tGenEt[iet]->Branch("val",   &genet, "val/F");
340	+	tGenEt[iet]->Branch("etcorr",&etcorr,"etcorr/F");
341	+	tGenEt[iet]->Branch("absrsp",&absrsp,"absrsp/F");
342	+	tGenEt[iet]->Branch("relrsp",&relrsp,"relrsp/F");
343	+	tGenEt[iet]->Branch("weight",&weight,"weight/F");
344	+
345	+	tEt[iet]=new TTree("tEt","tEt");
346	+	tEt[iet]->Branch("val",   &et,    "val/F");
347	+	tEt[iet]->Branch("etcorr",&etcorr,"etcorr/F");
348	+	tEt[iet]->Branch("absrsp",&absrsp,"absrsp/F");
349	+	tEt[iet]->Branch("relrsp",&relrsp,"relrsp/F");
350	+	tEt[iet]->Branch("weight",&weight,"weight/F");
351	+	}
352	+	for (int ieta=0;ieta<netabins;ieta++) {
353	+	tEta[ieta]=new TTree("tEta","tEta");
354	+	tEta[ieta]->Branch("val",   &eta,   "val/F");
355	+	tEta[ieta]->Branch("etcorr",&etcorr,"etcorr/F");
356	+	tEta[ieta]->Branch("absrsp",&absrsp,"absrsp/F");
357	+	tEta[ieta]->Branch("relrsp",&relrsp,"relrsp/F");
358	+	tEta[ieta]->Branch("weight",&weight,"weight/F");
359	+	}
360	+	for (int iphi=0;iphi<nphibins;iphi++) {
361	+	tPhi[iphi]=new TTree("tPhi","tPhi");
362	+	tPhi[iphi]->Branch("val",   &phi,   "val/F");
363	+	tPhi[iphi]->Branch("etcorr",&etcorr,"etcorr/F");
364	+	tPhi[iphi]->Branch("absrsp",&absrsp,"absrsp/F");
365	+	tPhi[iphi]->Branch("relrsp",&relrsp,"relrsp/F");
366	+	tPhi[iphi]->Branch("weight",&weight,"weight/F");
367	+	}
368	+	for (int iemf=0;iemf<nemfbins;iemf++) {
369	+	tEmf[iemf]=new TTree("tEmf","tEmf");
370	+	tEmf[iemf]->Branch("val",   &emf,   "val/F");
371	+	tEmf[iemf]->Branch("etcorr",&etcorr,"etcorr/F");
372	+	tEmf[iemf]->Branch("absrsp",&absrsp,"absrsp/F");
373	+	tEmf[iemf]->Branch("relrsp",&relrsp,"relrsp/F");
374	+	tEmf[iemf]->Branch("weight",&weight,"weight/F");
375	+	}
376	+	for (int ieta=0;ieta<netabins;ieta++) {
377	+	tEtEtaEt[ieta] = new TTree*[netbins];
378	+	for (int iet=0;iet<netbins;iet++) {
379	+	tEtEtaEt[ieta][iet] = new TTree("tEtEtaEt","tEtEtaEt");
380	+	tEtEtaEt[ieta][iet]->Branch("val",   &et,    "val/F");
381	+	tEtEtaEt[ieta][iet]->Branch("etcorr",&etcorr,"etcorr/F");
382	+	tEtEtaEt[ieta][iet]->Branch("absrsp",&absrsp,"absrsp/F");
383	+	tEtEtaEt[ieta][iet]->Branch("relrsp",&relrsp,"relrsp/F");
384	+	tEtEtaEt[ieta][iet]->Branch("weight",&weight,"weight/F");
385	+	}
386	+	}
387	+	for (int iemf=0;iemf<nemfbins;iemf++) {
388	+	tEtEmfEt[iemf] = new TTree*[netbins];
389	+	for (int iet=0;iet<netbins;iet++) {
390	+	tEtEmfEt[iemf][iet] = new TTree("tEtEmfEt","tEtEmfEt");
391	+	tEtEmfEt[iemf][iet]->Branch("val",   &et,    "val/F");
392	+	tEtEmfEt[iemf][iet]->Branch("etcorr",&etcorr,"etcorr/F");
393	+	tEtEmfEt[iemf][iet]->Branch("absrsp",&absrsp,"absrsp/F");
394	+	tEtEmfEt[iemf][iet]->Branch("relrsp",&relrsp,"relrsp/F");
395	+	tEtEmfEt[iemf][iet]->Branch("weight",&weight,"weight/F");
396	+	}
397	+	}
398	+	for (int iemf=0;iemf<nemfbins;iemf++) {
399	+	tEtEmfEtaEt[iemf] = new TTree**[netabins];
400	+	for (int ieta=0;ieta<netabins;ieta++) {
401	+	tEtEmfEtaEt[iemf][ieta] = new TTree*[netbins];
402	+	for (int iet=0;iet<netbins;iet++) {
403	+	tEtEmfEtaEt[iemf][ieta][iet] = new TTree("tEtEmfEtaEt","tEtEmfEtaEt");
404	+	tEtEmfEtaEt[iemf][ieta][iet]->Branch("val",   &et,    "val/F");
405	+	tEtEmfEtaEt[iemf][ieta][iet]->Branch("etcorr",&etcorr,"etcorr/F");
406	+	tEtEmfEtaEt[iemf][ieta][iet]->Branch("absrsp",&absrsp,"absrsp/F");
407	+	tEtEmfEtaEt[iemf][ieta][iet]->Branch("relrsp",&relrsp,"relrsp/F");
408	+	tEtEmfEtaEt[iemf][ieta][iet]->Branch("weight",&weight,"weight/F");
409	+	}
410	+	}
411	+	}
412	+
413	+	cout<<jetalgs[i]<<": trees created."<<endl;
414	+
415
416		// loop over all events
417		for (int ievt=0;ievt<nevts;ievt++) {
#	Line 249 | Line 423 | int main(int argc,char**argv)
423
424		if (jtgendr[ijt]<0.0||jtgendr[ijt]>drmax) continue;
425
426	<	if (ijt>0) continue;
426	>	//if (ijt>ijtmax) continue;
427	>
428	>	et      = jtet[ijt];
429	>	genet   = jtgenet[ijt];
430	>	eta     = (abseta) ? std::abs(jteta[ijt]) : jteta[ijt];
431	>	phi     = jtphi[ijt];
432	>	emf     = jtemf[ijt];
433
434	<	float pt      = jtpt[ijt];
255	<	float genpt   = jtgenpt[ijt];
256	<	float eta     = (abseta) ? std::abs(jteta[ijt]) : jteta[ijt];
257	<	float phi     = jtphi[ijt];
434	>	etcorr  = (applyjes==0) ? et : et*jtjes[ijt][applyjes-1];
435
436	<	float ptcorr  = (applyjes==0) ? pt : pt*jtjes[ijt][applyjes-1];
437	<	float deltapt = genpt - ptcorr;
436	>	absrsp  = genet - etcorr;
437	>	relrsp  = etcorr/genet;
438
439	<	// genpt
263	<	int igenpt = hist::get_ibin(genpt,ptbins);
264	<	hGenPt[i][igenpt]      ->Fill(genpt,weight);
265	<	//hGenPtPtCorr[i][igenpt]->Fill(ptcorr,weight);
266	<	hGenPtPtCorr[i][igenpt]->Fill(genpt,weight);
267	<	hRspVsGenPt[i][igenpt] ->Fill(deltapt,weight);
439	>	if (TMath::IsNaN(emf)) emf = 0.0;
440
441	<	// pt
442	<	int ipt = hist::get_ibin(pt,ptbins);
443	<	hPt[i][ipt]     ->Fill(pt,weight);
444	<	hPtCorr[i][ipt] ->Fill(ptcorr,weight);
445	<	hRspVsPt[i][ipt]->Fill(deltapt,weight);
441	>	int igenet = (netbins>0)  ? hist::get_ibin(genet,etbins) : -1;
442	>	int iet    = (netbins>0)  ? hist::get_ibin(et,etbins)    : -1;
443	>	int ieta   = (netabins>0) ? hist::get_ibin(eta,etabins)  : -1;
444	>	int iphi   = (nphibins>0) ? hist::get_ibin(phi,phibins)  : -1;
445	>	int iemf   = (nemfbins>0) ? hist::get_ibin(emf,emfbins)  : -1;
446
447	<	// eta
448	<	int ieta = hist::get_ibin(eta,etabins);
449	<	hEta[i][ieta]      ->Fill(eta,    weight);
450	<	hEtaPtCorr[i][ieta]->Fill(ptcorr, weight);
451	<	hRspVsEta[i][ieta] ->Fill(deltapt,weight);
452	<
453	<	// phi
454	<	int iphi = hist::get_ibin(phi,phibins);
455	<	hPhi[i][iphi]      ->Fill(phi,    weight);
456	<	hPhiPtCorr[i][iphi]->Fill(ptcorr, weight);
457	<	hRspVsPhi[i][iphi] ->Fill(deltapt,weight);
286	<
287	<	// pt-eta
288	<	hEtaPtPt[i][ieta][ipt]    ->Fill(pt,     weight);
289	<	hEtaPtPtCorr[i][ieta][ipt]->Fill(ptcorr, weight);
290	<	hRspVsEtaPt[i][ieta][ipt] ->Fill(deltapt,weight);
291	<
447	>	if (netbins>0)  tGenEt[igenet]->Fill();
448	>	if (netbins>0)  tEt[iet]->Fill();
449	>	if (netabins>0) tEta[ieta]->Fill();
450	>	if (nphibins>0) tPhi[iphi]->Fill();
451	>	if (nemfbins>0) tEmf[iemf]->Fill();
452	>
453	>	if (netbins>0&&netabins>0) tEtEtaEt[ieta][iet]->Fill();
454	>	if (netbins>0&&nemfbins>0) tEtEmfEt[iemf][iet]->Fill();
455	>
456	>	if (netbins>0&&netabins>0&&nemfbins>0)
457	>	tEtEmfEtaEt[iemf][ieta][iet]->Fill();
458
459		} // jets
460		} // evts
461
462		cout<<jetalgs[i]<<": all events processed / histograms filled."<<endl;
297	–
463
464	<	TDirectory* d = plotfile->mkdir(jetalgs[i].c_str());
464	>
465	>	// replace histograms
466	>
467	>	replaceHistos(netbins,
468	>	hGenEt,hEtCorrGenEt,hAbsRspGenEt,hRelRspGenEt,tGenEt);
469	>	replaceHistos(netbins,hEt,hEtCorrEt,hAbsRspEt,hRelRspEt,tEt);
470	>	replaceHistos(netabins,hEta,hEtCorrEta,hAbsRspEta,hRelRspEta,tEta);
471	>	replaceHistos(nphibins,hPhi,hEtCorrPhi,hAbsRspPhi,hRelRspPhi,tPhi);
472	>	replaceHistos(nemfbins,hEmf,hEtCorrEmf,hAbsRspEmf,hRelRspEmf,tEmf);
473	>	replaceHistos(netabins,netbins,
474	>	hEtEtaEt,hEtCorrEtaEt,hAbsRspEtaEt,hRelRspEtaEt,
475	>	tEtEtaEt);
476	>	replaceHistos(nemfbins,netbins,
477	>	hEtEmfEt,hEtCorrEmfEt,hAbsRspEmfEt,hRelRspEmfEt,
478	>	tEtEmfEt);
479	>	replaceHistos(nemfbins,netabins,netbins,
480	>	hEtEmfEtaEt,hEtCorrEmfEtaEt,
481	>	hAbsRspEmfEtaEt,hRelRspEmfEtaEt,
482	>	tEtEmfEtaEt);
483	>
484	>	cout<<jetalgs[i]<<": all histograms replaced, and tree deleted."<<endl;
485	>
486	>
487	>	TDirectory* d = plotfile->GetDirectory(jetalgs[i].c_str());
488	>	if (0!=d) plotfile->rmdir(jetalgs[i].c_str());
489	>	d = plotfile->mkdir(jetalgs[i].c_str());
490		d->cd();
491
492	<
493	<	// response / resolution vs *genpt*
494	<	gRspVsGenPt.push_back(new TGraphErrors(nptbins));
305	<	gSgmVsGenPt.push_back(new TGraphErrors(nptbins));
306	<	gRspVsGenPt.back()->SetName(("RspVsGenPt_"+jetalgs[i]).c_str());
307	<	gSgmVsGenPt.back()->SetName(("SgmVsGenPt_"+jetalgs[i]).c_str());
308	<	fillRspAndSgm(gRspVsGenPt.back(),gSgmVsGenPt.back(),nsigma,
309	<	ptbins,hGenPt[i],hGenPtPtCorr[i],hRspVsGenPt[i]);
310	<
311	<	// response / resolution vs *pt*
312	<	gRspVsPt.push_back(new TGraphErrors(nptbins));
313	<	gSgmVsPt.push_back(new TGraphErrors(nptbins));
314	<	gRspVsPt.back()->SetName(("RspVsPt_"+jetalgs[i]).c_str());
315	<	gSgmVsPt.back()->SetName(("SgmVsPt_"+jetalgs[i]).c_str());
316	<	fillRspAndSgm(gRspVsPt.back(),gSgmVsPt.back(),nsigma,
317	<	ptbins,hPt[i],hPtCorr[i],hRspVsPt[i]);
318	<
319	<	// response / resolution vs *eta*
320	<	gRspVsEta.push_back(new TGraphErrors(netabins));
321	<	gSgmVsEta.push_back(new TGraphErrors(netabins));
322	<	gRspVsEta.back()->SetName(("RspVsEta_"+jetalgs[i]).c_str());
323	<	gSgmVsEta.back()->SetName(("SgmVsEta_"+jetalgs[i]).c_str());
324	<	fillRspAndSgm(gRspVsEta.back(),gSgmVsEta.back(),nsigma,
325	<	etabins,hEta[i],hEtaPtCorr[i],hRspVsEta[i]);
326	<
327	<	// response / resolution vs *phi*
328	<	gRspVsPhi.push_back(new TGraphErrors(nphibins));
329	<	gSgmVsPhi.push_back(new TGraphErrors(nphibins));
330	<	gRspVsPhi.back()->SetName(("RspVsPhi_"+jetalgs[i]).c_str());
331	<	gSgmVsPhi.back()->SetName(("SgmVsPhi_"+jetalgs[i]).c_str());
332	<	fillRspAndSgm(gRspVsPhi.back(),gSgmVsPhi.back(),nsigma,
333	<	phibins,hPhi[i],hPhiPtCorr[i],hRspVsPhi[i]);
334	<
335	<	// response / resolution vs *pt-eta*
336	<	TGraphErrors** g;
337	<	g=new TGraphErrors*[netabins];
338	<	for (int ieta=0;ieta<netabins;ieta++) g[ieta] = new TGraphErrors(nptbins);
339	<	gRspVsEtaPt.push_back(g);
340	<	g=new TGraphErrors*[netabins];
341	<	for (int ieta=0;ieta<netabins;ieta++) g[ieta] = new TGraphErrors(nptbins);
342	<	gSgmVsEtaPt.push_back(g);
343	<	for (int ieta=0;ieta<netabins;ieta++) {
344	<	gRspVsEtaPt[i][ieta]->SetName(("RspVsPhi_"+
345	<	hist::get_bin_name("Eta",ieta,etabins)+
346	<	"_"+jetalgs[i]).c_str());
347	<	gSgmVsEtaPt[i][ieta]->SetName(("SgmVsPhi_"+
348	<	hist::get_bin_name("Eta",ieta,etabins)+
349	<	"_"+jetalgs[i]).c_str());
350	<	fillRspAndSgm(gRspVsEtaPt[i][ieta],gSgmVsEtaPt[i][ieta],nsigma,ptbins,
351	<	hEtaPtPt[i][ieta],hEtaPtPtCorr[i][ieta],hRspVsEtaPt[i][ieta]);
352	<	}
353	<	}
354	<
355	<
356	<	TCanvas* c = (TCanvas*)gROOT->FindObject("c1");
357	<	if (0!=c) delete c;
358	<
492	>	// fit response histos
493	>	fitHistos(hAbsRspGenEt,netbins,nsigma);
494	>	fitHistos(hRelRspGenEt,netbins,nsigma);
495
496	<	//
497	<	// make plots
362	<	//
363	<	for (unsigned int ialg=0;ialg<jetalgs.size();ialg++) {
364	<
365	<	TDirectory* d = plotfile->GetDirectory(jetalgs[ialg].c_str());
366	<	d->cd();
496	>	fitHistos(hAbsRspEt,   netbins,nsigma);
497	>	fitHistos(hRelRspEt,   netbins,nsigma);
498
499	<	// genpt
500	<	TCanvas* cGenPt       = plotVarHistos(ptbins,hGenPt[ialg],true);
501	<	TCanvas* cGenPtPtCorr = plotVarHistos(ptbins,hGenPtPtCorr[ialg],true);
502	<	TCanvas* cRspGenPt    = plotRspHistos(ptbins,hRspVsGenPt[ialg]);
503	<
373	<	delete cGenPt;
374	<	delete cGenPtPtCorr;
375	<	delete cRspGenPt;
376	<
377	<	// pt
378	<	TCanvas* cPt     = plotVarHistos(ptbins,hPt[ialg],true);
379	<	TCanvas* cPtCorr = plotVarHistos(ptbins,hPtCorr[ialg],true);
380	<	TCanvas* cRspPt  = plotRspHistos(ptbins,hRspVsPt[ialg]);
381	<
382	<	delete cPt;
383	<	delete cPtCorr;
384	<	delete cRspPt;
385	<
386	<	// eta
387	<	TCanvas* cEta       = plotVarHistos(etabins,hEta[ialg]);
388	<	TCanvas* cEtaPtCorr = plotVarHistos(etabins,hEtaPtCorr[ialg],true);
389	<	TCanvas* cRspEta    = plotRspHistos(etabins,hRspVsEta[ialg]);
390	<
391	<	delete cEta;
392	<	delete cEtaPtCorr;
393	<	delete cRspEta;
394	<
395	<	// phi
396	<	TCanvas* cPhi       = plotVarHistos(phibins,hPhi[ialg]);
397	<	TCanvas* cPhiPtCorr = plotVarHistos(phibins,hPhiPtCorr[ialg],true);
398	<	TCanvas* cRspPhi    = plotRspHistos(phibins,hRspVsPhi[ialg]);
399	<
400	<	delete cPhi;
401	<	delete cPhiPtCorr;
402	<	delete cRspPhi;
499	>	fitHistos(hAbsRspEta,  netabins,nsigma);
500	>	fitHistos(hRelRspEta,  netabins,nsigma);
501	>
502	>	fitHistos(hAbsRspPhi,  nphibins,nsigma);
503	>	fitHistos(hRelRspPhi,  nphibins,nsigma);
504
505	<	// eta-pt
506	<	for (int ieta=0;ieta<netabins;ieta++) {
406	<	TCanvas* cEtaPtPt     = plotVarHistos(ptbins,hEtaPtPt[ialg][ieta],true);
407	<	TCanvas* cEtaPtPtCorr = plotVarHistos(ptbins,hEtaPtPtCorr[ialg][ieta],true);
408	<	TCanvas* cRspEtaPt    = plotRspHistos(ptbins,hRspVsEtaPt[ialg][ieta]);
409	<
410	<	delete cEtaPtPt;
411	<	delete cEtaPtPtCorr;
412	<	delete cRspEtaPt;
413	<	}
505	>	fitHistos(hAbsRspEmf,  nemfbins,nsigma);
506	>	fitHistos(hRelRspEmf,  nemfbins,nsigma);
507
508	<	}
509	<
510	<	plotfile->cd();
508	>	fitHistos(hAbsRspEtaEt,netabins,netbins,nsigma);
509	>	fitHistos(hRelRspEtaEt,netabins,netbins,nsigma);
510	>
511	>	fitHistos(hAbsRspEmfEt,nemfbins,netbins,nsigma);
512	>	fitHistos(hRelRspEmfEt,nemfbins,netbins,nsigma);
513
514	<	// response & resolution vs genpT
515	<	TCanvas* cRspVsGenPt=plotResponse(jetalgs,applyjes,gRspVsGenPt,
516	<	"GenPt","jet p_{T}^{gen} [GeV]");
517	<	TCanvas* cSgmVsGenPt=plotResolution(jetalgs,gSgmVsGenPt,
423	<	"GenPt","jet p_{T}^{gen} [GeV]");
424	<	cRspVsGenPt->cd();
425	<	cSgmVsGenPt->cd();
426	<
427	<	// response & resolution vs pT
428	<	TCanvas* cRspVsPt  = plotResponse(jetalgs,applyjes,gRspVsPt,"Pt","jet p_{T} [GeV]");
429	<	TCanvas* cSgmVsPt  = plotResolution(jetalgs,gSgmVsPt,"Pt","jet p_{T} [GeV]");
430	<	cRspVsPt->cd();
431	<	cSgmVsPt->cd();
432	<
433	<	// response & resolution vs eta
434	<	TCanvas* cRspVsEta = plotResponse(jetalgs,applyjes,gRspVsEta,"Eta","jet #eta");
435	<	TCanvas* cSgmVsEta = plotResolution(jetalgs,gSgmVsEta,"Eta","jet #eta");
436	<	cRspVsEta->cd();
437	<	cSgmVsEta->cd();
438	<
439	<	// response & resolution vs phi
440	<	TCanvas* cRspVsPhi = plotResponse(jetalgs,applyjes,gRspVsPhi,"Phi","jet #phi");
441	<	TCanvas* cSgmVsPhi = plotResolution(jetalgs,gSgmVsPhi,"Phi","jet #phi");
442	<	cRspVsPhi->cd();
443	<	cSgmVsPhi->cd();
514	>	fitHistos(hAbsRspEmfEtaEt,nemfbins,netabins,netbins,nsigma);
515	>	fitHistos(hRelRspEmfEtaEt,nemfbins,netabins,netbins,nsigma);
516	>
517	>	cout<<jetalgs[i]<<": fits done."<<endl;
518
445	–	// response & resolution vs eta/pT
446	–	for (unsigned int ialg=0;ialg<jetalgs.size();ialg++) {
519
520	<	TDirectory* d = plotfile->GetDirectory(jetalgs[ialg].c_str());
521	<	d->cd();
520	>	// save histograms
521	>	saveHistos(hGenEt,      netbins);
522	>	saveHistos(hEtCorrGenEt,netbins);
523	>	saveHistos(hAbsRspGenEt,netbins);
524	>	saveHistos(hRelRspGenEt,netbins);
525	>
526	>	saveHistos(hEt,         netbins);
527	>	saveHistos(hEtCorrEt,   netbins);
528	>	saveHistos(hAbsRspEt,   netbins);
529	>	saveHistos(hRelRspEt,   netbins);
530	>
531	>	saveHistos(hEta,        netabins);
532	>	saveHistos(hEtCorrEta,  netabins);
533	>	saveHistos(hAbsRspEta,  netabins);
534	>	saveHistos(hRelRspEta,  netabins);
535	>
536	>	saveHistos(hPhi,        nphibins);
537	>	saveHistos(hEtCorrPhi,  nphibins);
538	>	saveHistos(hAbsRspPhi,  nphibins);
539	>	saveHistos(hRelRspPhi,  nphibins);
540	>
541	>	saveHistos(hEmf,        nemfbins);
542	>	saveHistos(hEtCorrEmf,  nemfbins);
543	>	saveHistos(hAbsRspEmf,  nemfbins);
544	>	saveHistos(hRelRspEmf,  nemfbins);
545
546	<	TCanvas* cRspVsEtaPt = plotResponse(jetalgs[ialg],applyjes,fitcorr,
547	<	etabins,gRspVsEtaPt[ialg],
548	<	"Eta","#eta","Pt","jet p_{T} [GeV]");
549	<	TCanvas* cSgmVsEtaPt = plotResolution(jetalgs[ialg],etabins,gSgmVsEtaPt[ialg],
550	<	"Eta","#eta","Pt","jet p_{T} [GeV]");
546	>	saveHistos(hEtEtaEt,    netabins,netbins);
547	>	saveHistos(hEtCorrEtaEt,netabins,netbins);
548	>	saveHistos(hAbsRspEtaEt,netabins,netbins);
549	>	saveHistos(hRelRspEtaEt,netabins,netbins);
550	>
551	>	saveHistos(hEtEmfEt,    nemfbins,netbins);
552	>	saveHistos(hEtCorrEmfEt,nemfbins,netbins);
553	>	saveHistos(hAbsRspEmfEt,nemfbins,netbins);
554	>	saveHistos(hRelRspEmfEt,nemfbins,netbins);
555
556	<	cRspVsEtaPt->cd();
557	<	cSgmVsEtaPt->cd();
556	>	saveHistos(hEtEmfEtaEt,    nemfbins,netabins,netbins);
557	>	saveHistos(hEtCorrEmfEtaEt,nemfbins,netabins,netbins);
558	>	saveHistos(hAbsRspEmfEtaEt,nemfbins,netabins,netbins);
559	>	saveHistos(hRelRspEmfEtaEt,nemfbins,netabins,netbins);
560	>
561	>	plotfile->Write();
562	>
563	>	cout<<jetalgs[i]<<": histograms saved."<<endl;
564		}
565
566		plotfile->Write();
567		plotfile->Close();
568		delete plotfile;
569
465	–	app->Run();
466	–
570		return 0;
571		}
572
#	Line 474 | Line 577 | int main(int argc,char**argv)
577		////////////////////////////////////////////////////////////////////////////////
578
579		//______________________________________________________________________________
580	<	void fillRspAndSgm(TGraphErrors* gRsp,TGraphErrors* gSgm,float nsigma,
581	<	const vector<float>& bins,TH1F** hVar,TH1F** hPt,TH1F** hRsp)
580	>	void replaceHistos(int nbins,
581	>	TH1F**& hVar,TH1F**& hEtCorr,TH1F**& hAbsRsp,TH1F**& hRelRsp,
582	>	TTree**& tVar)
583		{
584	<	int nbins = bins.size()+1;
584	>	if (nbins==0) return;
585	>
586	>	TH1::SetDefaultSumw2();
587		for (int ibin=0;ibin<nbins;ibin++) {
588
589	<	float var = hVar[ibin]->GetMean();
590	<	float pt  = hPt[ibin]->GetMean();
591	<	TH1F* h   = hRsp[ibin];
592	<
593	<	float de       = h->GetMean();
594	<	float errde    = h->GetMeanError();
595	<	float sgmde    = h->GetRMS();
596	<	float errsgmde = h->GetRMSError();
597	<
598	<	if (h->GetEntries()>50) {
599	<	string fname = "fitfnc" + (string)h->GetName();
600	<	TF1* f = new TF1(fname.c_str(),"gaus",de-nsigma*sgmde,de+nsigma*sgmde);
601	<	f->SetLineColor(kRed);
602	<	f->SetLineWidth(1);
603	<	f->SetParameter(1,de);
604	<	f->SetParameter(2,sgmde);
605	<	h->Fit(f,"QR");
606	<
607	<	de       = f->GetParameter(1);
608	<	errde    = f->GetParError(1);
609	<	sgmde    = f->GetParameter(2);
610	<	errsgmde = f->GetParError(2);
611	<	}
612	<
613	<	float rsp    = pt/(pt+de);
614	<	float errrsp = pt/(pt+de)/(pt+de)*errde;
615	<	float sgm    = pt/(pt+de)/(pt+de)*sgmde;
616	<	float errsgm = 2*pt*sgmde/(pt+de)/(pt+de)/(pt+de)*errde;
617	<
618	<	sgm    /= rsp;
619	<	errsgm /= rsp;
620	<
621	<	gRsp->SetPoint(ibin,var,rsp);
622	<	gRsp->SetPointError(ibin,0.0,errrsp);
623	<	gSgm->SetPoint(ibin,var,sgm);
624	<	gSgm->SetPointError(ibin,0.0,errsgm);
589	>	TH1F* h(0);
590	>
591	>	if (tVar[ibin]->GetEntries()==0) continue;
592	>
593	>	tVar[ibin]->Project("h(50)","val","weight*(1)");
594	>	h = (TH1F*)gROOT->FindObject("h");
595	>	h->SetNameTitle(hVar[ibin]->GetName(),hVar[ibin]->GetTitle());
596	>	h->SetXTitle(hVar[ibin]->GetXaxis()->GetTitle());
597	>	h->SetYTitle(hVar[ibin]->GetYaxis()->GetTitle());
598	>	delete hVar[ibin];
599	>	hVar[ibin] = h;
600	>
601	>	tVar[ibin]->Project("h(50)","etcorr","weight*(1)");
602	>	h = (TH1F*)gROOT->FindObject("h");
603	>	h->SetNameTitle(hEtCorr[ibin]->GetName(),hEtCorr[ibin]->GetTitle());
604	>	h->SetXTitle(hEtCorr[ibin]->GetXaxis()->GetTitle());
605	>	h->SetYTitle(hEtCorr[ibin]->GetYaxis()->GetTitle());
606	>	delete hEtCorr[ibin];
607	>	hEtCorr[ibin] = h;
608	>
609	>	tVar[ibin]->Project("h(50)","absrsp","weight*(1)");
610	>	h = (TH1F*)gROOT->FindObject("h");
611	>	h->SetNameTitle(hAbsRsp[ibin]->GetName(),hAbsRsp[ibin]->GetTitle());
612	>	h->SetXTitle(hAbsRsp[ibin]->GetXaxis()->GetTitle());
613	>	h->SetYTitle(hAbsRsp[ibin]->GetYaxis()->GetTitle());
614	>	delete hAbsRsp[ibin];
615	>	hAbsRsp[ibin] = h;
616	>
617	>	tVar[ibin]->Project("h(50)","relrsp","weight*(1)");
618	>	h = (TH1F*)gROOT->FindObject("h");
619	>	h->SetNameTitle(hRelRsp[ibin]->GetName(),hRelRsp[ibin]->GetTitle());
620	>	h->SetXTitle(hRelRsp[ibin]->GetXaxis()->GetTitle());
621	>	h->SetYTitle(hRelRsp[ibin]->GetYaxis()->GetTitle());
622	>	delete hRelRsp[ibin];
623	>	hRelRsp[ibin] = h;
624	>
625	>	delete tVar[ibin];
626		}
627
628	<	gRsp->Write();
522	<	gSgm->Write();
628	>	delete [] tVar;
629		}
630
631
632		//______________________________________________________________________________
633	<	TCanvas* plotVarHistos(const vector<float>& bins,TH1F** hVar,bool logy)
633	>	void replaceHistos(int nbins1,int nbins2,
634	>	TH1F***& hVar,TH1F***& hEtCorr,TH1F***& hAbsRsp,TH1F***& hRelRsp,
635	>	TTree***& tVar)
636		{
637	<	string::size_type pos;
530	<	string cname,jetalg;
531	<
532	<	cname  = hVar[0]->GetName();
533	<	pos    = cname.find_last_of('_');
534	<	jetalg = cname.substr(pos+1);
535	<	cname  = cname.substr(0,pos);
536	<	pos    = cname.find_last_of(':');
537	<	cname  = cname.substr(0,pos);
538	<	cname += "_"+jetalg;
539	<
540	<	TCanvas* c = new TCanvas(cname.c_str(),cname.c_str());
541	<
542	<	int nbins = bins.size()+1;
543	<	int ncx= int(std::sqrt((float)nbins));
544	<	int ncy= ncx;
545	<	if (ncx*ncy<nbins) ncx++;
546	<	if (ncx*ncy<nbins) ncy++;
547	<
548	<	c->Divide(ncx,ncy);
549	<
550	<	for (int ibin=0;ibin<nbins;ibin++) {
551	<	c->cd(ibin+1);
552	<	gPad->SetLogy(logy);
553	<	hVar[ibin]->Draw("EHIST");
554	<	hVar[ibin]->Write();
555	<	}
556	<
557	<	//c->Write();
637	>	if (nbins1==0) return;
638
639	<	return c;
639	>	for (int i1=0;i1<nbins1;i1++)
640	>	replaceHistos(nbins2,hVar[i1],hEtCorr[i1],hAbsRsp[i1],hRelRsp[i1],tVar[i1]);
641	>	delete [] tVar;
642		}
643
644
645		//______________________________________________________________________________
646	<	TCanvas* plotRspHistos(const vector<float>& bins,TH1F** hRsp)
646	>	void replaceHistos(int nbins1,int nbins2,int nbins3,
647	>	TH1F****& hVar,TH1F****& hEtCorr,TH1F****&hAbsRsp,TH1F****&hRelRsp,
648	>	TTree****& tVar)
649		{
650	<	string::size_type pos;
651	<	string cname,jetalg;
652	<
653	<
654	<	cname  = hRsp[0]->GetName();
571	<	pos    = cname.find_last_of('_');
572	<	jetalg = cname.substr(pos+1);
573	<	cname  = cname.substr(0,pos);
574	<	pos    = cname.find_last_of('_');
575	<	cname  = cname.substr(0,pos);
576	<	cname += "_"+jetalg;
577	<
578	<	TCanvas* c = new TCanvas(cname.c_str(),cname.c_str());
579	<
580	<	int nbins = bins.size()+1;
581	<	int ncx= int(std::sqrt((float)nbins));
582	<	int ncy= ncx;
583	<	if (ncx*ncy<nbins) ncx++;
584	<	if (ncx*ncy<nbins) ncy++;
585	<
586	<	c->Divide(ncx,ncy);
587	<
588	<	for (int ibin=0;ibin<nbins;ibin++) {
589	<	c->cd(ibin+1);
590	<	if (hRsp[ibin]->GetEntries()>0) gPad->SetLogy();
591	<	hRsp[ibin]->Draw("EHIST");
592	<	string fname = "fitfnc" + (string)hRsp[ibin]->GetName();
593	<	TF1* fitfnc=hRsp[ibin]->GetFunction(fname.c_str());
594	<	if (0!=fitfnc) fitfnc->Draw("SAME");
595	<	hRsp[ibin]->Write();
596	<	}
597	<
598	<	c->Write();
599	<
600	<	return c;
650	>	if (nbins1==0) return;
651	>	for (int i1=0;i1<nbins1;i1++)
652	>	replaceHistos(nbins2,nbins3,
653	>	hVar[i1],hEtCorr[i1],hAbsRsp[i1],hRelRsp[i1],tVar[i1]);
654	>	delete [] tVar;
655		}
656
657
658		//______________________________________________________________________________
659	<	TCanvas* plotResponse(const vector<string>& jetalgs,short applyjes,
606	<	const vector<TGraphErrors*> gRsp,
607	<	const string& varname,const string& xtitle)
659	>	void fitHistos(TH1F** histos,int n,double nsigma)
660		{
661	<	assert(jetalgs.size()==gRsp.size());
662	<
663	<	string   cname = "RspVs"+varname;
664	<	TCanvas* c     = new TCanvas(cname.c_str(),cname.c_str());
665	<	TLegend* l     = new TLegend(0.7,0.25+jetalgs.size()*0.045,0.85,0.25);
666	<
667	<	l->SetFillColor(10);
668	<	Color_t color=kBlack;
669	<	TMultiGraph* mg = new TMultiGraph();
670	<
671	<	for (unsigned int i=0;i<gRsp.size();i++) {
620	<	TGraphErrors* g = gRsp[i];
621	<	l->AddEntry(g,jetalgs[i].c_str(),"l");
622	<	while (color==kYellow||color==kWhite||color==10) color++;
623	<	g->SetMarkerColor(color);
624	<	g->SetLineColor(color);
625	<	g->SetMarkerStyle(kFullCircle);
626	<	g->SetMarkerSize(0.5);
627	<	mg->Add(g,"P");
628	<
629	<	if (varname=="Pt") {
630	<	string fname = "fitfnc"+cname;
631	<	string fitfncAsString =
632	<	//(applyjes==0) ? "[0]+[1]*log(x)+[2]*log(x)*log(x)+[3]*x+[4]*x*x" : "pol0";
633	<	(applyjes==0) ? "[0]*(pow(log(x),[1])-1/x)" : "pol0";
634	<
635	<	TF1* fitfnc = new TF1(fname.c_str(),fitfncAsString.c_str(),1,get_xmax(g));
636	<
637	<	if (applyjes==0) {
638	<	fitfnc->SetParameter(0,0.4);
639	<	fitfnc->SetParameter(1,0.4);
640	<	}
641	<	else {
642	<	fitfnc->SetParameter(0,1.0);
643	<	}
644	<
645	<	fitfnc->SetLineWidth(1);
646	<	fitfnc->SetLineColor(color);
647	<	g->Fit(fitfnc,"QR");
648	<	}
649	<
650	<	g->Write();
651	<
652	<	color++;
661	>	for (int i=0;i<n;i++) {
662	>	float nrm  = histos[i]->Integral();
663	>	float mean = histos[i]->GetMean();
664	>	float rms  = histos[i]->GetRMS();
665	>	TF1* f = new TF1("gaus","gaus",mean-nsigma*rms,mean+nsigma*rms);
666	>	f->SetParameter(0,nrm);
667	>	f->SetParameter(1,mean);
668	>	f->SetParameter(2,rms);
669	>	f->SetLineWidth(1);
670	>	f->SetLineColor(kRed);
671	>	histos[i]->Fit(f,"QR");
672		}
654	–
655	–	string mgtitle = "Jet Response Vs "+varname;
656	–	mg->Draw("a");
657	–	mg->SetTitle(mgtitle.c_str());
658	–	mg->GetXaxis()->SetTitle(xtitle.c_str());
659	–	mg->GetYaxis()->SetTitle("E_{T}^{rec}/E_{T}^{gen}");
660	–	if (varname=="Pt") mg->GetXaxis()->SetLimits(0.0,mg->GetXaxis()->GetXmax());
661	–	mg->SetMinimum(0.0);
662	–	mg->SetMaximum(1.1);
663	–	l->Draw();
664	–
665	–	c->Write();
666	–
667	–	return c;
673		}
674
675
676		//______________________________________________________________________________
677	<	TCanvas* plotResponse(const string& jetalg,short applyjes,bool fitcorr,
673	<	const vector<float>& bins1,TGraphErrors** gRsp,
674	<	const string& varname1,const string vartitle1,
675	<	const string& varname2,const string& xtitle)
676	<
677	>	void fitHistos(TH1F*** histos,int n1,int n2,double nsigma)
678		{
679	<	int      nbins = bins1.size()+1;
679	<
680	<	string   fitfncAsString = "pol0";
681	<	int      fitfncNbParams = 1;
682	<
683	<	ofstream fout;
684	<
685	<	if (applyjes==0) {
686	<	//fitfncAsString = "[0]+[1]*log(x)+[2]*log(x)*log(x)+[3]*x+[4]*x*x";
687	<	//fitfncNbParams = 5;
688	<	fitfncAsString = "[0]*(pow(log(x),[1])-1/x)";
689	<	fitfncNbParams = 2;
690	<	}
691	<	if (fitcorr) {
692	<	fout.open((jetalg+".txt").c_str()); assert(fout.is_open());
693	<	fout<<"#\n# jet correction parameters for "<<jetalg<<" algorithm\n#\n\n"
694	<	<<"# number of eta bins\n"<<nbins
695	<	<<"\n\n# fit function\n"<<fitfncAsString<<"\n"<<fitfncNbParams
696	<	<<endl<<endl;
697	<	}
698	<
699	<	string   cname = "RspVs"+varname1+varname2+"_"+jetalg;
700	<	TCanvas* c = new TCanvas(cname.c_str(),cname.c_str());
701	<	TLegend* l = new TLegend(0.7,0.25+nbins*0.045,0.85,0.25);
702	<
703	<	l->SetFillColor(10);
704	<	Color_t color=kBlack;
705	<	TMultiGraph* mg = new TMultiGraph();
706	<
707	<	for (int ibin=0;ibin<nbins;ibin++) {
708	<
709	<	TGraphErrors* g = gRsp[ibin];
710	<	l->AddEntry(g,hist::get_bin_title(vartitle1,ibin,bins1).c_str(),"l");
711	<	while (color==kYellow||color==kWhite||color==10) color++;
712	<	g->SetMarkerColor(color);
713	<	g->SetLineColor(color);
714	<	g->SetMarkerStyle(kFullCircle);
715	<	g->SetMarkerSize(0.5);
716	<	mg->Add(g,"P");
717	<
718	<	if (varname2=="Pt") {
719	<	string fname = "fitfnc"+cname;
720	<	TF1* fitfnc = new TF1(fname.c_str(),fitfncAsString.c_str(),1,get_xmax(g));
721	<
722	<	if (applyjes==0) {
723	<	fitfnc->SetParameter(0,0.4);
724	<	fitfnc->SetParameter(1,0.4);
725	<	}
726	<	else {
727	<	fitfnc->SetParameter(0,1.0);
728	<	}
729	<
730	<	fitfnc->SetLineWidth(1);
731	<	fitfnc->SetLineColor(color);
732	<	g->Fit(fitfnc,"QR");
733	<
734	<	if (fitcorr==0) {
735	<	if (ibin==nbins-1)
736	<	fout<<"# eta > "<<bins1.back()<<"\n999.9"<<endl;
737	<	else
738	<	fout<<"# eta < "<<bins1[ibin]<<"\n"<<bins1[ibin]<<endl;
739	<	for (int ipar=0;ipar<fitfnc->GetNpar();ipar++)
740	<	fout<<setw(12)<<fitfnc->GetParameter(ipar)<<" "
741	<	<<setw(12)<<fitfnc->GetParError(ipar)<<endl;
742	<	fout<<endl;
743	<	}
744	<
745	<	}
746	<
747	<	g->Write();
748	<
749	<	color++;
750	<	}
751	<
752	<	if (fout.is_open()) fout.close();
753	<
754	<	string mgtitle = jetalg+": Jet Response vs "+varname1+" and "+varname2;
755	<	mg->Draw("a");
756	<	mg->SetTitle(mgtitle.c_str());
757	<	mg->GetXaxis()->SetTitle(xtitle.c_str());
758	<	mg->GetYaxis()->SetTitle("E_{T}^{rec}/E_{T}^{gen}");
759	<	if (varname2=="Pt") mg->GetXaxis()->SetLimits(0.0,mg->GetXaxis()->GetXmax());
760	<	mg->SetMinimum(0.0);
761	<	mg->SetMaximum(1.1);
762	<	l->Draw();
763	<
764	<	c->Write();
765	<
766	<	return c;
679	>	for (int i1=0;i1<n1;i1++) fitHistos(histos[i1],n2,nsigma);
680		}
681
682
683		//______________________________________________________________________________
684	<	TCanvas* plotResolution(const vector<string>& jetalgs,
772	<	const vector<TGraphErrors*> gSgm,
773	<	const string& varname,const string& xtitle)
684	>	void fitHistos(TH1F**** histos,int n1,int n2,int n3,double nsigma)
685		{
686	<	string   cname = "SgmVs"+varname;
687	<	TCanvas* c     = new TCanvas(cname.c_str(),cname.c_str());
777	<	TLegend* l     = new TLegend(0.7,0.85-jetalgs.size()*0.045,0.85,0.85);
778	<
779	<	l->SetFillColor(10);
780	<	Color_t color=kBlack;
781	<	TMultiGraph* mg = new TMultiGraph();
782	<
783	<	for (unsigned int i=0;i<gSgm.size();i++) {
784	<	TGraphErrors* g = gSgm[i];
785	<	l->AddEntry(g,jetalgs[i].c_str(),"l");
786	<	while (color==kYellow||color==kWhite||color==10) color++;
787	<	g->SetMarkerColor(color);
788	<	g->SetLineColor(color);
789	<	g->SetMarkerStyle(kFullCircle);
790	<	g->SetMarkerSize(0.5);
791	<	mg->Add(g,"P");
792	<
793	<	if (varname=="Pt") {
794	<	string fname = "fitfnc"+cname;
795	<	TF1* fitfnc = new TF1(fname.c_str(),"sqrt([0]*[0]/x*x+[1]*[1]/x+[2]*[2])");
796	<	fitfnc->SetParameter(0,0.0);
797	<	fitfnc->SetParameter(1,0.5);
798	<	fitfnc->SetParameter(2,0.1);
799	<	fitfnc->SetLineWidth(1);
800	<	fitfnc->SetLineColor(color);
801	<	fitfnc->SetRange(2.0,get_xmax(g));
802	<	g->Fit(fitfnc,"QR");
803	<	}
686	>	for (int i1=0;i1<n1;i1++) fitHistos(histos[i1],n2,n3,nsigma);
687	>	}
688
689	<	g->Write();
690	<
691	<	color++;
689	>
690	>	//______________________________________________________________________________
691	>	void saveHistos(TH1F**& histos,int n)
692	>	{
693	>	if (n==0) return;
694	>	for (int i=0;i<n;i++) {
695	>	histos[i]->Write();
696	>	delete histos[i];
697		}
698	<
810	<	string mgtitle = "Jet Resolution Vs "+varname;
811	<	mg->Draw("a");
812	<	mg->SetTitle(mgtitle.c_str());
813	<	mg->GetXaxis()->SetTitle(xtitle.c_str());
814	<	mg->GetYaxis()
815	<	->SetTitle("#sigma(E_{T}^{rec}/E_{T}^{gen})/<E_{T}^{rec}/E_{T}^{gen}>");
816	<	if (varname=="Pt") mg->GetXaxis()->SetLimits(0.0,mg->GetXaxis()->GetXmax());
817	<	mg->SetMinimum(0.0);
818	<	mg->SetMaximum(0.5);
819	<	l->Draw();
820	<
821	<	c->Write();
822	<
823	<	return c;
698	>	delete [] histos;
699		}
700
701
702		//______________________________________________________________________________
703	<	TCanvas* plotResolution(const string& jetalg,const vector<float>& bins1,
829	<	TGraphErrors** gSgm,
830	<	const string& varname1,const string vartitle1,
831	<	const string& varname2,const string& xtitle)
703	>	void saveHistos(TH1F***& histos,int n1,int n2)
704		{
705	<	int      nbins = bins1.size()+1;
706	<
707	<	string   cname = "SgmVs"+varname1+varname2+"_"+jetalg;
836	<	TCanvas* c     = new TCanvas(cname.c_str(),cname.c_str());
837	<	TLegend* l     = new TLegend(0.7,0.85-nbins*0.045,0.85,0.85);
838	<
839	<	l->SetFillColor(10);
840	<	Color_t color=kBlack;
841	<	TMultiGraph* mg = new TMultiGraph();
842	<
843	<	for (int ibin=0;ibin<nbins;ibin++) {
844	<	TGraphErrors* g = gSgm[ibin];
845	<	l->AddEntry(g,hist::get_bin_title(vartitle1,ibin,bins1).c_str(),"l");
846	<	while (color==kYellow||color==kWhite||color==10) color++;
847	<	g->SetMarkerColor(color);
848	<	g->SetLineColor(color);
849	<	g->SetMarkerStyle(kFullCircle);
850	<	g->SetMarkerSize(0.5);
851	<	mg->Add(g,"P");
852	<
853	<	if (varname2=="Pt") {
854	<	string fname = "fitfnc"+cname;
855	<	TF1* fitfnc = new TF1(fname.c_str(),"sqrt([0]*[0]/x*x+[1]*[1]/x+[2]*[2])");
856	<	fitfnc->SetParameter(0,0.0);
857	<	fitfnc->SetParameter(1,0.5);
858	<	fitfnc->SetParameter(2,0.1);
859	<	fitfnc->SetLineWidth(1);
860	<	fitfnc->SetLineColor(color);
861	<	fitfnc->SetRange(2.0,get_xmax(g));
862	<	g->Fit(fitfnc,"QR");
863	<	}
864	<
865	<	g->Write();
866	<
867	<	color++;
868	<	}
869	<
870	<	string mgtitle=jetalg+": Jet Resolution Vs "+varname1+" and "+varname2;
871	<	mg->Draw("a");
872	<	mg->SetTitle(mgtitle.c_str());
873	<	mg->GetXaxis()->SetTitle(xtitle.c_str());
874	<	mg->GetYaxis()
875	<	->SetTitle("#sigma(E_{T}^{rec}/E_{T}^{gen})/<E_{T}^{rec}/E_{T}^{gen}>");
876	<	if (varname2=="Pt") mg->GetXaxis()->SetLimits(0.0,mg->GetXaxis()->GetXmax());
877	<	mg->SetMinimum(0.0);
878	<	mg->SetMaximum(0.5);
879	<	l->Draw();
880	<
881	<	c->Write();
882	<
883	<	return c;
705	>	if (n1==0) return;
706	>	for (int i1=0;i1<n1;i1++) saveHistos(histos[i1],n2);
707	>	delete [] histos;
708		}
709
710
711		//______________________________________________________________________________
712	<	float get_xmax(TGraph* g)
712	>	void saveHistos(TH1F****& histos,int n1,int n2,int n3)
713		{
714	<	float result(-1e200);
715	<	for (int i=0;i<g->GetN();++i) if (g->GetX()[i]>result) result=g->GetX()[i];
716	<	return result;
714	>	if (n1==0) return;
715	>	for (int i1=0;i1<n1;i1++) saveHistos(histos[i1],n2,n3);
716	>	delete [] histos;
717		}

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines