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.1 by schiefer, Wed Aug 15 17:20:06 2007 UTC vs.
Revision 1.13 by schiefer, Mon Sep 3 17:20:56 2007 UTC

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

Diff Legend

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