10 |
|
#include "utils/cmdline.h" |
11 |
|
#include "utils/hist.h" |
12 |
|
|
13 |
+ |
#include <TSystem.h> |
14 |
|
#include <TROOT.h> |
15 |
|
#include <TRint.h> |
16 |
|
#include <TFile.h> |
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> |
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 |
– |
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); |
53 |
|
|
54 |
|
|
55 |
|
//////////////////////////////////////////////////////////////////////////////// |
63 |
|
cl.parse(argc,argv); |
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 |
+ |
short applyjes = cl.get_value<short> ("applyjes", 0); |
70 |
|
string treename = cl.get_value<string> ("treename", "t"); |
71 |
|
float nsigma = cl.get_value<float> ("nsigma", 3.0); |
72 |
< |
short applyjes = cl.get_value<short> ("applyjes", 0); |
73 |
< |
bool fitcorr = cl.get_value<bool> ("fitcorr", false); |
74 |
< |
vector<float> ptbins = cl.get_vector<float> ("ptbins", dflt_ptbins); |
75 |
< |
vector<float> etabins = cl.get_vector<float> ("etabins",dflt_etabins); |
90 |
< |
vector<float> phibins = cl.get_vector<float> ("phibins",dflt_phibins); |
72 |
> |
vector<float> etbins = cl.get_vector<float> ("etbins", ""); |
73 |
> |
vector<float> etabins = cl.get_vector<float> ("etabins", ""); |
74 |
> |
vector<float> phibins = cl.get_vector<float> ("phibins", ""); |
75 |
> |
vector<float> emfbins = cl.get_vector<float> ("emfbins", ""); |
76 |
|
bool abseta = cl.get_value<bool> ("abseta", true); |
77 |
|
|
78 |
|
if (!cl.check()) return 0; |
79 |
|
cl.print(); |
80 |
|
|
81 |
< |
|
81 |
> |
int netbins = (etbins.size() >0) ? etbins.size() +1 : 0; |
82 |
> |
int netabins = (etabins.size()>0) ? etabins.size()+1 : 0; |
83 |
> |
int nphibins = (phibins.size()>0) ? phibins.size()+1 : 0; |
84 |
> |
int nemfbins = (emfbins.size()>0) ? emfbins.size()+1 : 0; |
85 |
> |
|
86 |
> |
if (netbins+netabins+nphibins+nemfbins==0) { |
87 |
> |
cout<<"Must specify bins: etbins, etabins, phibins, AND/OR emfbins!"<<endl; |
88 |
> |
return 0; |
89 |
> |
} |
90 |
> |
|
91 |
|
// etabins |
92 |
< |
if (!abseta) { |
92 |
> |
if (netabins>0&&!abseta) { |
93 |
|
int neta=(int)etabins.size(); |
94 |
|
std::reverse(etabins.begin(),etabins.end()); |
95 |
|
for (int ieta=neta-1;ieta>=0;ieta--) etabins.push_back(etabins[ieta]); |
110 |
|
// prepare branch variables and histograms / graphs |
111 |
|
float weight; |
112 |
|
char njt; |
113 |
< |
float jtpt[100]; |
113 |
> |
float jtet[100]; |
114 |
|
float jteta[100]; |
115 |
|
float jtphi[100]; |
116 |
< |
float jtgenpt[100]; |
116 |
> |
float jtemf[100]; |
117 |
> |
float jtgenet[100]; |
118 |
|
float jtgendr[100]; |
119 |
|
float jtjes[100][3]; |
120 |
|
|
121 |
< |
int nptbins = ptbins.size()+1; |
122 |
< |
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); |
121 |
> |
argc=3; argv[1]="-l"; argv[2]="-b"; |
122 |
> |
TRint* app = new TRint(argv[0],&argc,argv); app->Argc(); |
123 |
|
|
124 |
|
|
125 |
|
// |
126 |
|
// loop over all input files |
127 |
|
// |
128 |
< |
TFile* plotfile = new TFile("mcresponse.root","RECREATE"); |
128 |
> |
TFile* plotfile = new TFile(output.c_str(),"UPDATE"); |
129 |
|
|
130 |
|
for (unsigned int i=0;i<input.size();++i) { |
131 |
|
TFile* f = new TFile(input[i].c_str(),"READ"); if (!f->IsOpen()) return 0; |
136 |
|
|
137 |
|
t->SetBranchAddress("weight", &weight); |
138 |
|
t->SetBranchAddress("njt", &njt); |
139 |
< |
t->SetBranchAddress("jtpt", jtpt); |
139 |
> |
t->SetBranchAddress("jtet", jtet); |
140 |
|
t->SetBranchAddress("jteta", jteta); |
141 |
|
t->SetBranchAddress("jtphi", jtphi); |
142 |
< |
t->SetBranchAddress("jtgenpt",jtgenpt); |
142 |
> |
t->SetBranchAddress("jtemf", jtemf); |
143 |
> |
t->SetBranchAddress("jtgenet",jtgenet); |
144 |
|
t->SetBranchAddress("jtgendr",jtgendr); |
145 |
< |
|
146 |
< |
if (t->FindBranch("jtjes")) t->SetBranchAddress("jtjes", jtjes); |
147 |
< |
else |
145 |
> |
|
146 |
> |
vector<TBranch*> branches; |
147 |
> |
branches.push_back(t->GetBranch("weight")); |
148 |
> |
branches.push_back(t->GetBranch("njt")); |
149 |
> |
branches.push_back(t->GetBranch("jtet")); |
150 |
> |
branches.push_back(t->GetBranch("jteta")); |
151 |
> |
branches.push_back(t->GetBranch("jtphi")); |
152 |
> |
branches.push_back(t->GetBranch("jtemf")); |
153 |
> |
branches.push_back(t->GetBranch("jtgenet")); |
154 |
> |
branches.push_back(t->GetBranch("jtgendr")); |
155 |
> |
|
156 |
> |
if (t->FindBranch("jtjes")) { |
157 |
> |
t->SetBranchAddress("jtjes", jtjes); |
158 |
> |
branches.push_back(t->GetBranch("jtjes")); |
159 |
> |
} |
160 |
> |
else { |
161 |
|
for (int i1=0;i1<100;i1++) |
162 |
|
for (int i2=0;i2<3;i2++) jtjes[i1][i2]=1.0; |
163 |
+ |
} |
164 |
+ |
string absrsp_xtitle = "E_{T}^{gen}-E_{T} [GeV]"; |
165 |
+ |
string relrsp_xtitle = "E_{T}^{rec}/E_{T}"; |
166 |
+ |
|
167 |
+ |
|
168 |
+ |
TH1F** hGenEt(0); |
169 |
+ |
TH1F** hEtCorrGenEt(0); |
170 |
+ |
TH1F** hAbsRspGenEt(0); |
171 |
+ |
TH1F** hRelRspGenEt(0); |
172 |
+ |
|
173 |
+ |
if (netbins>0) { |
174 |
+ |
hGenEt = hist::initHistos("GenEt",jetalgs[i],100, |
175 |
+ |
"jet E_{T}^{gen} [GeV]", |
176 |
+ |
"GenEt",etbins); |
177 |
+ |
hEtCorrGenEt = hist::initHistos("EtCorrGenEt",jetalgs[i],100, |
178 |
+ |
"jet E_{T}^{gen} [GeV]", |
179 |
+ |
"GenEt",etbins); |
180 |
+ |
hAbsRspGenEt = hist::initHistos("AbsRspGenEt",jetalgs[i], |
181 |
+ |
100,-50,150, |
182 |
+ |
absrsp_xtitle, |
183 |
+ |
"GenEt",etbins); |
184 |
+ |
hRelRspGenEt = hist::initHistos("RelRspGenEt",jetalgs[i], |
185 |
+ |
100,0,2, |
186 |
+ |
relrsp_xtitle, |
187 |
+ |
"GenEt",etbins); |
188 |
+ |
} |
189 |
+ |
|
190 |
+ |
|
191 |
+ |
TH1F** hEt(0); |
192 |
+ |
TH1F** hEtCorrEt(0); |
193 |
+ |
TH1F** hAbsRspEt(0); |
194 |
+ |
TH1F** hRelRspEt(0); |
195 |
+ |
|
196 |
+ |
if (netbins>0) { |
197 |
+ |
hEt = hist::initHistos("Et",jetalgs[i],100, |
198 |
+ |
"jet E_{T} [GeV]","Et",etbins); |
199 |
+ |
hEtCorrEt = hist::initHistos("EtCorrEt",jetalgs[i],100, |
200 |
+ |
"jet E_{T} [GeV]","Et",etbins); |
201 |
+ |
hAbsRspEt = hist::initHistos("AbsRspEt",jetalgs[i],100,-50,150, |
202 |
+ |
absrsp_xtitle,"Et",etbins); |
203 |
+ |
hRelRspEt = hist::initHistos("RelRspEt",jetalgs[i],100,0,2, |
204 |
+ |
relrsp_xtitle,"Et",etbins); |
205 |
+ |
} |
206 |
+ |
|
207 |
+ |
|
208 |
+ |
TH1F** hEta(0); |
209 |
+ |
TH1F** hEtCorrEta(0); |
210 |
+ |
TH1F** hAbsRspEta(0); |
211 |
+ |
TH1F** hRelRspEta(0); |
212 |
+ |
|
213 |
+ |
if (netabins>0) { |
214 |
+ |
hEta = hist::initHistos("Eta",jetalgs[i],100, |
215 |
+ |
"jet #eta","Eta",etabins); |
216 |
+ |
hEtCorrEta = hist::initHistos("EtCorrEta",jetalgs[i],100, |
217 |
+ |
"jet E_{T} [GeV]","Eta",etabins); |
218 |
+ |
hAbsRspEta = hist::initHistos("AbsRspEta",jetalgs[i],100,-50,150, |
219 |
+ |
absrsp_xtitle,"Eta",etabins); |
220 |
+ |
hRelRspEta = hist::initHistos("RelRspEta",jetalgs[i],100,0,2, |
221 |
+ |
relrsp_xtitle,"Eta",etabins); |
222 |
+ |
} |
223 |
|
|
224 |
+ |
|
225 |
+ |
TH1F** hPhi(0); |
226 |
+ |
TH1F** hEtCorrPhi(0); |
227 |
+ |
TH1F** hAbsRspPhi(0); |
228 |
+ |
TH1F** hRelRspPhi(0); |
229 |
+ |
|
230 |
+ |
if (nphibins>0) { |
231 |
+ |
hPhi = hist::initHistos("Phi",jetalgs[i],100, |
232 |
+ |
"jet #phi","Phi",phibins); |
233 |
+ |
hEtCorrPhi = hist::initHistos("EtCorrPhi",jetalgs[i],100, |
234 |
+ |
"jet E_{T} [Gev]","Phi",phibins); |
235 |
+ |
hAbsRspPhi = hist::initHistos("AbsRspPhi",jetalgs[i],100,-50,150, |
236 |
+ |
absrsp_xtitle,"Phi",phibins); |
237 |
+ |
hRelRspPhi = hist::initHistos("RelRspPhi",jetalgs[i],100,0,2, |
238 |
+ |
relrsp_xtitle,"Phi",phibins); |
239 |
+ |
} |
240 |
+ |
|
241 |
+ |
|
242 |
+ |
TH1F** hEmf(0); |
243 |
+ |
TH1F** hEtCorrEmf(0); |
244 |
+ |
TH1F** hAbsRspEmf(0); |
245 |
+ |
TH1F** hRelRspEmf(0); |
246 |
+ |
|
247 |
+ |
if (nemfbins>0) { |
248 |
+ |
hEmf = hist::initHistos("Emf",jetalgs[i],100, |
249 |
+ |
"jet emf","Emf",emfbins); |
250 |
+ |
hEtCorrEmf = hist::initHistos("EtCorrEmf",jetalgs[i],100, |
251 |
+ |
"jet E_{T} [Gev]","Emf",emfbins); |
252 |
+ |
hAbsRspEmf = hist::initHistos("AbsRspEmf",jetalgs[i],100,-50,150, |
253 |
+ |
absrsp_xtitle,"Emf",emfbins); |
254 |
+ |
hRelRspEmf = hist::initHistos("RelRspEmf",jetalgs[i],100,0,2, |
255 |
+ |
relrsp_xtitle,"Emf",emfbins); |
256 |
+ |
} |
257 |
|
|
189 |
– |
string rsp_xtitle = "E_{T}^{gen}-E_{T}^{rec} [GeV]"; |
258 |
|
|
259 |
< |
hGenPt .push_back(hist::initHistos("JetGenPt",jetalgs[i], |
260 |
< |
4000,0,4000,"jet p_{T}^{gen} [GeV]", |
261 |
< |
"GenPt",ptbins)); |
262 |
< |
hGenPtPtCorr.push_back(hist::initHistos("JetPtCorr",jetalgs[i], |
263 |
< |
4000,0,4000,"jet p_{T}^{gen} [GeV]", |
264 |
< |
"GenPt",ptbins)); |
265 |
< |
hRspVsGenPt .push_back(hist::initHistos("JetRspGenPt",jetalgs[i], |
266 |
< |
100,-100,100,rsp_xtitle, |
267 |
< |
"GenPt",ptbins)); |
268 |
< |
|
269 |
< |
hPt .push_back(hist::initHistos("JetPt",jetalgs[i], |
270 |
< |
4000,0,4000,"jet p_{T} [GeV]", |
271 |
< |
"Pt",ptbins)); |
272 |
< |
hPtCorr .push_back(hist::initHistos("JetPtCorr",jetalgs[i], |
273 |
< |
4000,0,4000,"jet p_{T} [GeV]", |
274 |
< |
"Pt",ptbins)); |
275 |
< |
hRspVsPt .push_back(hist::initHistos("JetRspPt",jetalgs[i], |
276 |
< |
100,-100,100,rsp_xtitle, |
277 |
< |
"Pt",ptbins)); |
278 |
< |
|
279 |
< |
hEta .push_back(hist::initHistos("JetEta",jetalgs[i], |
280 |
< |
100,0,5,"jet #eta","Eta",etabins)); |
281 |
< |
hEtaPtCorr .push_back(hist::initHistos("JetPtCorr",jetalgs[i], |
282 |
< |
4000,0,4000,"jet p_{T} [GeV]", |
283 |
< |
"Eta",etabins)); |
284 |
< |
hRspVsEta .push_back(hist::initHistos("JetRspEta",jetalgs[i], |
285 |
< |
100,-100,100,rsp_xtitle,"Eta",etabins)); |
286 |
< |
|
287 |
< |
hPhi .push_back(hist::initHistos("JetPhi",jetalgs[i], |
288 |
< |
100,-3.2,3.2,"jet #phi","Phi",phibins)); |
289 |
< |
hPhiPtCorr .push_back(hist::initHistos("JetPtCorr",jetalgs[i], |
290 |
< |
4000,0,4000,"jet p_{T} [Gev]", |
291 |
< |
"Phi",phibins)); |
292 |
< |
hRspVsPhi .push_back(hist::initHistos("JetRspPhi",jetalgs[i], |
293 |
< |
100,-100,100,rsp_xtitle,"Phi",phibins)); |
294 |
< |
|
295 |
< |
hEtaPtPt .push_back(hist::initHistos("JetPt",jetalgs[i], |
296 |
< |
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)); |
259 |
> |
TH1F*** hEtEtaEt(0); |
260 |
> |
TH1F*** hEtCorrEtaEt(0); |
261 |
> |
TH1F*** hAbsRspEtaEt(0); |
262 |
> |
TH1F*** hRelRspEtaEt(0); |
263 |
> |
|
264 |
> |
if (netbins>0&&netabins>0) { |
265 |
> |
hEtEtaEt = hist::initHistos("Et",jetalgs[i],100,"jet E_{T} [GeV]", |
266 |
> |
"Eta",etabins,"Et",etbins); |
267 |
> |
hEtCorrEtaEt = hist::initHistos("EtCorrEtaEt",jetalgs[i],100, |
268 |
> |
"jet E_{T} [GeV]", |
269 |
> |
"Eta",etabins,"Et",etbins); |
270 |
> |
hAbsRspEtaEt = hist::initHistos("AbsRspEtaEt",jetalgs[i],100,-50,150, |
271 |
> |
absrsp_xtitle, |
272 |
> |
"Eta",etabins,"Et",etbins); |
273 |
> |
hRelRspEtaEt = hist::initHistos("RelRspEtaEt",jetalgs[i],100,0,2, |
274 |
> |
relrsp_xtitle, |
275 |
> |
"Eta",etabins,"Et",etbins); |
276 |
> |
} |
277 |
> |
|
278 |
> |
|
279 |
> |
TH1F*** hEtEmfEt(0); |
280 |
> |
TH1F*** hEtCorrEmfEt(0); |
281 |
> |
TH1F*** hAbsRspEmfEt(0); |
282 |
> |
TH1F*** hRelRspEmfEt(0); |
283 |
> |
|
284 |
> |
if (netbins>0&&nemfbins>0) { |
285 |
> |
hEtEmfEt = hist::initHistos("Et",jetalgs[i],100,"jet E_{T} [GeV]", |
286 |
> |
"Emf",emfbins,"Et",etbins); |
287 |
> |
hEtCorrEmfEt = hist::initHistos("EtCorrEmfEt",jetalgs[i],100, |
288 |
> |
"jet E_{T} [GeV]", |
289 |
> |
"Emf",emfbins,"Et",etbins); |
290 |
> |
hAbsRspEmfEt = hist::initHistos("AbsRspEmfEt",jetalgs[i],100,-50,150, |
291 |
> |
absrsp_xtitle, |
292 |
> |
"Emf",emfbins,"Et",etbins); |
293 |
> |
hRelRspEmfEt = hist::initHistos("RelRspEmfEt",jetalgs[i],100,0,2, |
294 |
> |
relrsp_xtitle, |
295 |
> |
"Emf",emfbins,"Et",etbins); |
296 |
> |
} |
297 |
|
|
298 |
+ |
|
299 |
+ |
TH1F**** hEtEmfEtaEt(0); |
300 |
+ |
TH1F**** hEtCorrEmfEtaEt(0); |
301 |
+ |
TH1F**** hAbsRspEmfEtaEt(0); |
302 |
+ |
TH1F**** hRelRspEmfEtaEt(0); |
303 |
+ |
|
304 |
+ |
if (netbins>0&&netabins>0&&nemfbins>0) { |
305 |
+ |
hEtEmfEtaEt = hist::initHistos("Et",jetalgs[i],100,"jet E_{T} [GeV]", |
306 |
+ |
"Emf",emfbins, |
307 |
+ |
"Eta",etabins, |
308 |
+ |
"Et", etbins); |
309 |
+ |
hEtCorrEmfEtaEt = hist::initHistos("EtCorrEmfEtaEt",jetalgs[i],100, |
310 |
+ |
"jet E_{T} [GeV]", |
311 |
+ |
"Emf",emfbins, |
312 |
+ |
"Eta",etabins, |
313 |
+ |
"Et", etbins); |
314 |
+ |
hAbsRspEmfEtaEt = hist::initHistos("AbsRspEmfEtaEt",jetalgs[i], |
315 |
+ |
100,-50,150,absrsp_xtitle, |
316 |
+ |
"Emf",emfbins, |
317 |
+ |
"Eta",etabins, |
318 |
+ |
"Et", etbins); |
319 |
+ |
hRelRspEmfEtaEt = hist::initHistos("RelRspEmfEtaEt",jetalgs[i], |
320 |
+ |
100,0,2,relrsp_xtitle, |
321 |
+ |
"Emf",emfbins, |
322 |
+ |
"Eta",etabins, |
323 |
+ |
"Et", etbins); |
324 |
+ |
} |
325 |
|
|
326 |
|
int nevts = el->GetN(); |
327 |
|
cout<<jetalgs[i]<<": "<<nevts<<" events selected."<<endl; |
328 |
|
|
329 |
+ |
// prepare intermediate trees |
330 |
+ |
float genet,et,eta,phi,emf,etcorr,absrsp,relrsp; |
331 |
+ |
|
332 |
+ |
TTree** tGenEt(0); |
333 |
+ |
TTree** tEt(0); |
334 |
+ |
TTree** tEta(0); |
335 |
+ |
TTree** tPhi(0); |
336 |
+ |
TTree** tEmf(0); |
337 |
+ |
TTree*** tEtEtaEt(0); |
338 |
+ |
TTree*** tEtEmfEt(0); |
339 |
+ |
TTree**** tEtEmfEtaEt(0); |
340 |
+ |
|
341 |
+ |
gROOT->cd(); |
342 |
+ |
|
343 |
+ |
if (netbins>0) tGenEt = new TTree*[netbins]; |
344 |
+ |
if (netbins>0) tEt = new TTree*[netbins]; |
345 |
+ |
if (netabins>0) tEta = new TTree*[netabins]; |
346 |
+ |
if (nphibins>0) tPhi = new TTree*[nphibins]; |
347 |
+ |
if (nemfbins>0) tEmf = new TTree*[nemfbins]; |
348 |
+ |
if (netbins>0&&netabins>0) tEtEtaEt = new TTree**[netabins]; |
349 |
+ |
if (netbins>0&&nemfbins>0) tEtEmfEt = new TTree**[nemfbins]; |
350 |
+ |
if (netbins>0&&netabins>0&&nemfbins>0) tEtEmfEtaEt = new TTree***[nemfbins]; |
351 |
|
|
352 |
+ |
for (int iet=0;iet<netbins;iet++) { |
353 |
+ |
|
354 |
+ |
tGenEt[iet]=new TTree("tGenEt","tGenEt"); |
355 |
+ |
tGenEt[iet]->Branch("val", &genet, "val/F"); |
356 |
+ |
tGenEt[iet]->Branch("etcorr",&etcorr,"etcorr/F"); |
357 |
+ |
tGenEt[iet]->Branch("absrsp",&absrsp,"absrsp/F"); |
358 |
+ |
tGenEt[iet]->Branch("relrsp",&relrsp,"relrsp/F"); |
359 |
+ |
tGenEt[iet]->Branch("weight",&weight,"weight/F"); |
360 |
+ |
|
361 |
+ |
tEt[iet]=new TTree("tEt","tEt"); |
362 |
+ |
tEt[iet]->Branch("val", &et, "val/F"); |
363 |
+ |
tEt[iet]->Branch("etcorr",&etcorr,"etcorr/F"); |
364 |
+ |
tEt[iet]->Branch("absrsp",&absrsp,"absrsp/F"); |
365 |
+ |
tEt[iet]->Branch("relrsp",&relrsp,"relrsp/F"); |
366 |
+ |
tEt[iet]->Branch("weight",&weight,"weight/F"); |
367 |
+ |
} |
368 |
+ |
for (int ieta=0;ieta<netabins;ieta++) { |
369 |
+ |
tEta[ieta]=new TTree("tEta","tEta"); |
370 |
+ |
tEta[ieta]->Branch("val", &eta, "val/F"); |
371 |
+ |
tEta[ieta]->Branch("etcorr",&etcorr,"etcorr/F"); |
372 |
+ |
tEta[ieta]->Branch("absrsp",&absrsp,"absrsp/F"); |
373 |
+ |
tEta[ieta]->Branch("relrsp",&relrsp,"relrsp/F"); |
374 |
+ |
tEta[ieta]->Branch("weight",&weight,"weight/F"); |
375 |
+ |
} |
376 |
+ |
for (int iphi=0;iphi<nphibins;iphi++) { |
377 |
+ |
tPhi[iphi]=new TTree("tPhi","tPhi"); |
378 |
+ |
tPhi[iphi]->Branch("val", &phi, "val/F"); |
379 |
+ |
tPhi[iphi]->Branch("etcorr",&etcorr,"etcorr/F"); |
380 |
+ |
tPhi[iphi]->Branch("absrsp",&absrsp,"absrsp/F"); |
381 |
+ |
tPhi[iphi]->Branch("relrsp",&relrsp,"relrsp/F"); |
382 |
+ |
tPhi[iphi]->Branch("weight",&weight,"weight/F"); |
383 |
+ |
} |
384 |
+ |
for (int iemf=0;iemf<nemfbins;iemf++) { |
385 |
+ |
tEmf[iemf]=new TTree("tEmf","tEmf"); |
386 |
+ |
tEmf[iemf]->Branch("val", &emf, "val/F"); |
387 |
+ |
tEmf[iemf]->Branch("etcorr",&etcorr,"etcorr/F"); |
388 |
+ |
tEmf[iemf]->Branch("absrsp",&absrsp,"absrsp/F"); |
389 |
+ |
tEmf[iemf]->Branch("relrsp",&relrsp,"relrsp/F"); |
390 |
+ |
tEmf[iemf]->Branch("weight",&weight,"weight/F"); |
391 |
+ |
} |
392 |
+ |
|
393 |
+ |
if (netbins>0) { |
394 |
+ |
for (int ieta=0;ieta<netabins;ieta++) { |
395 |
+ |
tEtEtaEt[ieta] = new TTree*[netbins]; |
396 |
+ |
for (int iet=0;iet<netbins;iet++) { |
397 |
+ |
tEtEtaEt[ieta][iet] = new TTree("tEtEtaEt","tEtEtaEt"); |
398 |
+ |
tEtEtaEt[ieta][iet]->Branch("val", &et, "val/F"); |
399 |
+ |
tEtEtaEt[ieta][iet]->Branch("etcorr",&etcorr,"etcorr/F"); |
400 |
+ |
tEtEtaEt[ieta][iet]->Branch("absrsp",&absrsp,"absrsp/F"); |
401 |
+ |
tEtEtaEt[ieta][iet]->Branch("relrsp",&relrsp,"relrsp/F"); |
402 |
+ |
tEtEtaEt[ieta][iet]->Branch("weight",&weight,"weight/F"); |
403 |
+ |
} |
404 |
+ |
} |
405 |
+ |
for (int iemf=0;iemf<nemfbins;iemf++) { |
406 |
+ |
tEtEmfEt[iemf] = new TTree*[netbins]; |
407 |
+ |
for (int iet=0;iet<netbins;iet++) { |
408 |
+ |
tEtEmfEt[iemf][iet] = new TTree("tEtEmfEt","tEtEmfEt"); |
409 |
+ |
tEtEmfEt[iemf][iet]->Branch("val", &et, "val/F"); |
410 |
+ |
tEtEmfEt[iemf][iet]->Branch("etcorr",&etcorr,"etcorr/F"); |
411 |
+ |
tEtEmfEt[iemf][iet]->Branch("absrsp",&absrsp,"absrsp/F"); |
412 |
+ |
tEtEmfEt[iemf][iet]->Branch("relrsp",&relrsp,"relrsp/F"); |
413 |
+ |
tEtEmfEt[iemf][iet]->Branch("weight",&weight,"weight/F"); |
414 |
+ |
} |
415 |
+ |
} |
416 |
+ |
if (netabins>0) { |
417 |
+ |
for (int iemf=0;iemf<nemfbins;iemf++) { |
418 |
+ |
tEtEmfEtaEt[iemf] = new TTree**[netabins]; |
419 |
+ |
for (int ieta=0;ieta<netabins;ieta++) { |
420 |
+ |
tEtEmfEtaEt[iemf][ieta] = new TTree*[netbins]; |
421 |
+ |
for (int iet=0;iet<netbins;iet++) { |
422 |
+ |
tEtEmfEtaEt[iemf][ieta][iet] = new TTree("tEtEmfEtaEt", |
423 |
+ |
"tEtEmfEtaEt"); |
424 |
+ |
tEtEmfEtaEt[iemf][ieta][iet]->Branch("val", &et, "val/F"); |
425 |
+ |
tEtEmfEtaEt[iemf][ieta][iet]->Branch("etcorr",&etcorr,"etcorr/F"); |
426 |
+ |
tEtEmfEtaEt[iemf][ieta][iet]->Branch("absrsp",&absrsp,"absrsp/F"); |
427 |
+ |
tEtEmfEtaEt[iemf][ieta][iet]->Branch("relrsp",&relrsp,"relrsp/F"); |
428 |
+ |
tEtEmfEtaEt[iemf][ieta][iet]->Branch("weight",&weight,"weight/F"); |
429 |
+ |
} |
430 |
+ |
} |
431 |
+ |
} |
432 |
+ |
} |
433 |
+ |
} |
434 |
+ |
|
435 |
+ |
cout<<jetalgs[i]<<": trees created."<<endl; |
436 |
+ |
|
437 |
|
// loop over all events |
438 |
|
for (int ievt=0;ievt<nevts;ievt++) { |
439 |
|
|
440 |
< |
t->GetEntry(el->GetEntry(ievt)); |
441 |
< |
|
440 |
> |
//t->GetEntry(el->GetEntry(ievt)); |
441 |
> |
for (int ibrnch=0;ibrnch<(int)branches.size();ibrnch++) |
442 |
> |
branches[ibrnch]->GetEntry(el->GetEntry(ievt)); |
443 |
> |
|
444 |
|
// loop over all jets in the event |
445 |
|
for (int ijt=0;ijt<njt;ijt++) { |
446 |
|
|
447 |
|
if (jtgendr[ijt]<0.0||jtgendr[ijt]>drmax) continue; |
448 |
|
|
449 |
< |
if (ijt>0) continue; |
449 |
> |
//if (ijt>ijtmax) continue; |
450 |
> |
|
451 |
> |
et = jtet[ijt]; |
452 |
> |
genet = jtgenet[ijt]; |
453 |
> |
eta = (abseta) ? std::abs(jteta[ijt]) : jteta[ijt]; |
454 |
> |
phi = jtphi[ijt]; |
455 |
> |
emf = jtemf[ijt]; |
456 |
|
|
457 |
< |
float pt = jtpt[ijt]; |
255 |
< |
float genpt = jtgenpt[ijt]; |
256 |
< |
float eta = (abseta) ? std::abs(jteta[ijt]) : jteta[ijt]; |
257 |
< |
float phi = jtphi[ijt]; |
457 |
> |
etcorr = (applyjes==0) ? et : et*jtjes[ijt][applyjes-1]; |
458 |
|
|
459 |
< |
float ptcorr = (applyjes==0) ? pt : pt*jtjes[ijt][applyjes-1]; |
460 |
< |
float deltapt = genpt - ptcorr; |
459 |
> |
absrsp = genet - etcorr; |
460 |
> |
relrsp = etcorr/genet; |
461 |
|
|
462 |
< |
// 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); |
462 |
> |
if (TMath::IsNaN(emf)) emf = 0.0; |
463 |
|
|
464 |
< |
// pt |
465 |
< |
int ipt = hist::get_ibin(pt,ptbins); |
466 |
< |
hPt[i][ipt] ->Fill(pt,weight); |
467 |
< |
hPtCorr[i][ipt] ->Fill(ptcorr,weight); |
468 |
< |
hRspVsPt[i][ipt]->Fill(deltapt,weight); |
464 |
> |
int igenet = (netbins>0) ? hist::get_ibin(genet,etbins) : -1; |
465 |
> |
int iet = (netbins>0) ? hist::get_ibin(et,etbins) : -1; |
466 |
> |
int ieta = (netabins>0) ? hist::get_ibin(eta,etabins) : -1; |
467 |
> |
int iphi = (nphibins>0) ? hist::get_ibin(phi,phibins) : -1; |
468 |
> |
int iemf = (nemfbins>0) ? hist::get_ibin(emf,emfbins) : -1; |
469 |
|
|
470 |
< |
// eta |
471 |
< |
int ieta = hist::get_ibin(eta,etabins); |
472 |
< |
hEta[i][ieta] ->Fill(eta, weight); |
473 |
< |
hEtaPtCorr[i][ieta]->Fill(ptcorr, weight); |
474 |
< |
hRspVsEta[i][ieta] ->Fill(deltapt,weight); |
475 |
< |
|
476 |
< |
// phi |
477 |
< |
int iphi = hist::get_ibin(phi,phibins); |
478 |
< |
hPhi[i][iphi] ->Fill(phi, weight); |
479 |
< |
hPhiPtCorr[i][iphi]->Fill(ptcorr, weight); |
480 |
< |
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 |
< |
|
470 |
> |
if (netbins>0) tGenEt[igenet]->Fill(); |
471 |
> |
if (netbins>0) tEt[iet]->Fill(); |
472 |
> |
if (netabins>0) tEta[ieta]->Fill(); |
473 |
> |
if (nphibins>0) tPhi[iphi]->Fill(); |
474 |
> |
if (nemfbins>0) tEmf[iemf]->Fill(); |
475 |
> |
|
476 |
> |
if (netbins>0&&netabins>0) tEtEtaEt[ieta][iet]->Fill(); |
477 |
> |
if (netbins>0&&nemfbins>0) tEtEmfEt[iemf][iet]->Fill(); |
478 |
> |
|
479 |
> |
if (netbins>0&&netabins>0&&nemfbins>0) |
480 |
> |
tEtEmfEtaEt[iemf][ieta][iet]->Fill(); |
481 |
|
|
482 |
|
} // jets |
483 |
|
} // evts |
484 |
|
|
485 |
|
cout<<jetalgs[i]<<": all events processed / histograms filled."<<endl; |
297 |
– |
|
486 |
|
|
487 |
< |
TDirectory* d = plotfile->mkdir(jetalgs[i].c_str()); |
487 |
> |
|
488 |
> |
// replace histograms |
489 |
> |
|
490 |
> |
replaceHistos(netbins, |
491 |
> |
hGenEt,hEtCorrGenEt,hAbsRspGenEt,hRelRspGenEt,tGenEt); |
492 |
> |
replaceHistos(netbins,hEt,hEtCorrEt,hAbsRspEt,hRelRspEt,tEt); |
493 |
> |
replaceHistos(netabins,hEta,hEtCorrEta,hAbsRspEta,hRelRspEta,tEta); |
494 |
> |
replaceHistos(nphibins,hPhi,hEtCorrPhi,hAbsRspPhi,hRelRspPhi,tPhi); |
495 |
> |
replaceHistos(nemfbins,hEmf,hEtCorrEmf,hAbsRspEmf,hRelRspEmf,tEmf); |
496 |
> |
replaceHistos(netabins,netbins, |
497 |
> |
hEtEtaEt,hEtCorrEtaEt,hAbsRspEtaEt,hRelRspEtaEt, |
498 |
> |
tEtEtaEt); |
499 |
> |
replaceHistos(nemfbins,netbins, |
500 |
> |
hEtEmfEt,hEtCorrEmfEt,hAbsRspEmfEt,hRelRspEmfEt, |
501 |
> |
tEtEmfEt); |
502 |
> |
replaceHistos(nemfbins,netabins,netbins, |
503 |
> |
hEtEmfEtaEt,hEtCorrEmfEtaEt, |
504 |
> |
hAbsRspEmfEtaEt,hRelRspEmfEtaEt, |
505 |
> |
tEtEmfEtaEt); |
506 |
> |
|
507 |
> |
cout<<jetalgs[i]<<": all histograms replaced, and tree deleted."<<endl; |
508 |
> |
|
509 |
> |
|
510 |
> |
plotfile->cd(); |
511 |
> |
TDirectory* d = plotfile->GetDirectory(jetalgs[i].c_str()); |
512 |
> |
if (0!=d) plotfile->rmdir(jetalgs[i].c_str()); |
513 |
> |
d = plotfile->mkdir(jetalgs[i].c_str()); |
514 |
|
d->cd(); |
515 |
|
|
516 |
+ |
// fit response histos |
517 |
+ |
fitHistos(hAbsRspGenEt,netbins,nsigma); |
518 |
+ |
fitHistos(hRelRspGenEt,netbins,nsigma); |
519 |
|
|
520 |
< |
// response / resolution vs *genpt* |
521 |
< |
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 |
< |
|
520 |
> |
fitHistos(hAbsRspEt, netbins,nsigma); |
521 |
> |
fitHistos(hRelRspEt, netbins,nsigma); |
522 |
|
|
523 |
< |
// |
524 |
< |
// make plots |
525 |
< |
// |
526 |
< |
for (unsigned int ialg=0;ialg<jetalgs.size();ialg++) { |
527 |
< |
|
528 |
< |
TDirectory* d = plotfile->GetDirectory(jetalgs[ialg].c_str()); |
529 |
< |
d->cd(); |
523 |
> |
fitHistos(hAbsRspEta, netabins,nsigma); |
524 |
> |
fitHistos(hRelRspEta, netabins,nsigma); |
525 |
> |
|
526 |
> |
fitHistos(hAbsRspPhi, nphibins,nsigma); |
527 |
> |
fitHistos(hRelRspPhi, nphibins,nsigma); |
528 |
> |
|
529 |
> |
fitHistos(hAbsRspEmf, nemfbins,nsigma); |
530 |
> |
fitHistos(hRelRspEmf, nemfbins,nsigma); |
531 |
> |
|
532 |
> |
fitHistos(hAbsRspEtaEt,netabins,netbins,nsigma); |
533 |
> |
fitHistos(hRelRspEtaEt,netabins,netbins,nsigma); |
534 |
> |
|
535 |
> |
fitHistos(hAbsRspEmfEt,nemfbins,netbins,nsigma); |
536 |
> |
fitHistos(hRelRspEmfEt,nemfbins,netbins,nsigma); |
537 |
> |
|
538 |
> |
fitHistos(hAbsRspEmfEtaEt,nemfbins,netabins,netbins,nsigma); |
539 |
> |
fitHistos(hRelRspEmfEtaEt,nemfbins,netabins,netbins,nsigma); |
540 |
> |
|
541 |
> |
cout<<jetalgs[i]<<": fits done."<<endl; |
542 |
|
|
368 |
– |
// genpt |
369 |
– |
TCanvas* cGenPt = plotVarHistos(ptbins,hGenPt[ialg],true); |
370 |
– |
TCanvas* cGenPtPtCorr = plotVarHistos(ptbins,hGenPtPtCorr[ialg],true); |
371 |
– |
TCanvas* cRspGenPt = plotRspHistos(ptbins,hRspVsGenPt[ialg]); |
372 |
– |
|
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; |
403 |
– |
|
404 |
– |
// eta-pt |
405 |
– |
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 |
– |
} |
543 |
|
|
544 |
< |
} |
545 |
< |
|
546 |
< |
plotfile->cd(); |
544 |
> |
// save histograms |
545 |
> |
saveHistos(hGenEt, netbins); |
546 |
> |
saveHistos(hEtCorrGenEt,netbins); |
547 |
> |
saveHistos(hAbsRspGenEt,netbins); |
548 |
> |
saveHistos(hRelRspGenEt,netbins); |
549 |
> |
|
550 |
> |
saveHistos(hEt, netbins); |
551 |
> |
saveHistos(hEtCorrEt, netbins); |
552 |
> |
saveHistos(hAbsRspEt, netbins); |
553 |
> |
saveHistos(hRelRspEt, netbins); |
554 |
|
|
555 |
< |
// response & resolution vs genpT |
556 |
< |
TCanvas* cRspVsGenPt=plotResponse(jetalgs,applyjes,gRspVsGenPt, |
557 |
< |
"GenPt","jet p_{T}^{gen} [GeV]"); |
558 |
< |
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(); |
555 |
> |
saveHistos(hEta, netabins); |
556 |
> |
saveHistos(hEtCorrEta, netabins); |
557 |
> |
saveHistos(hAbsRspEta, netabins); |
558 |
> |
saveHistos(hRelRspEta, netabins); |
559 |
|
|
560 |
< |
// response & resolution vs eta/pT |
561 |
< |
for (unsigned int ialg=0;ialg<jetalgs.size();ialg++) { |
560 |
> |
saveHistos(hPhi, nphibins); |
561 |
> |
saveHistos(hEtCorrPhi, nphibins); |
562 |
> |
saveHistos(hAbsRspPhi, nphibins); |
563 |
> |
saveHistos(hRelRspPhi, nphibins); |
564 |
|
|
565 |
< |
TDirectory* d = plotfile->GetDirectory(jetalgs[ialg].c_str()); |
566 |
< |
d->cd(); |
565 |
> |
saveHistos(hEmf, nemfbins); |
566 |
> |
saveHistos(hEtCorrEmf, nemfbins); |
567 |
> |
saveHistos(hAbsRspEmf, nemfbins); |
568 |
> |
saveHistos(hRelRspEmf, nemfbins); |
569 |
|
|
570 |
< |
TCanvas* cRspVsEtaPt = plotResponse(jetalgs[ialg],applyjes,fitcorr, |
571 |
< |
etabins,gRspVsEtaPt[ialg], |
572 |
< |
"Eta","#eta","Pt","jet p_{T} [GeV]"); |
573 |
< |
TCanvas* cSgmVsEtaPt = plotResolution(jetalgs[ialg],etabins,gSgmVsEtaPt[ialg], |
574 |
< |
"Eta","#eta","Pt","jet p_{T} [GeV]"); |
570 |
> |
saveHistos(hEtEtaEt, netabins,netbins); |
571 |
> |
saveHistos(hEtCorrEtaEt,netabins,netbins); |
572 |
> |
saveHistos(hAbsRspEtaEt,netabins,netbins); |
573 |
> |
saveHistos(hRelRspEtaEt,netabins,netbins); |
574 |
> |
|
575 |
> |
saveHistos(hEtEmfEt, nemfbins,netbins); |
576 |
> |
saveHistos(hEtCorrEmfEt,nemfbins,netbins); |
577 |
> |
saveHistos(hAbsRspEmfEt,nemfbins,netbins); |
578 |
> |
saveHistos(hRelRspEmfEt,nemfbins,netbins); |
579 |
> |
|
580 |
> |
saveHistos(hEtEmfEtaEt, nemfbins,netabins,netbins); |
581 |
> |
saveHistos(hEtCorrEmfEtaEt,nemfbins,netabins,netbins); |
582 |
> |
saveHistos(hAbsRspEmfEtaEt,nemfbins,netabins,netbins); |
583 |
> |
saveHistos(hRelRspEmfEtaEt,nemfbins,netabins,netbins); |
584 |
|
|
585 |
< |
cRspVsEtaPt->cd(); |
458 |
< |
cSgmVsEtaPt->cd(); |
585 |
> |
cout<<jetalgs[i]<<": histograms saved."<<endl; |
586 |
|
} |
587 |
|
|
588 |
|
plotfile->Write(); |
589 |
|
plotfile->Close(); |
590 |
|
delete plotfile; |
591 |
< |
|
592 |
< |
app->Run(); |
591 |
> |
|
592 |
> |
//app->Run(); |
593 |
|
|
594 |
|
return 0; |
595 |
|
} |
601 |
|
//////////////////////////////////////////////////////////////////////////////// |
602 |
|
|
603 |
|
//______________________________________________________________________________ |
604 |
< |
void fillRspAndSgm(TGraphErrors* gRsp,TGraphErrors* gSgm,float nsigma, |
605 |
< |
const vector<float>& bins,TH1F** hVar,TH1F** hPt,TH1F** hRsp) |
604 |
> |
void replaceHistos(int nbins, |
605 |
> |
TH1F**& hVar,TH1F**& hEtCorr,TH1F**& hAbsRsp,TH1F**& hRelRsp, |
606 |
> |
TTree**& tVar) |
607 |
|
{ |
608 |
< |
int nbins = bins.size()+1; |
608 |
> |
if (nbins==0) return; |
609 |
> |
|
610 |
> |
TH1::SetDefaultSumw2(); |
611 |
|
for (int ibin=0;ibin<nbins;ibin++) { |
612 |
|
|
613 |
< |
float var = hVar[ibin]->GetMean(); |
614 |
< |
float pt = hPt[ibin]->GetMean(); |
615 |
< |
TH1F* h = hRsp[ibin]; |
616 |
< |
|
617 |
< |
float de = h->GetMean(); |
618 |
< |
float errde = h->GetMeanError(); |
619 |
< |
float sgmde = h->GetRMS(); |
620 |
< |
float errsgmde = h->GetRMSError(); |
621 |
< |
|
622 |
< |
if (h->GetEntries()>50) { |
623 |
< |
string fname = "fitfnc" + (string)h->GetName(); |
624 |
< |
TF1* f = new TF1(fname.c_str(),"gaus",de-nsigma*sgmde,de+nsigma*sgmde); |
625 |
< |
f->SetLineColor(kRed); |
626 |
< |
f->SetLineWidth(1); |
627 |
< |
f->SetParameter(1,de); |
628 |
< |
f->SetParameter(2,sgmde); |
629 |
< |
h->Fit(f,"QR"); |
630 |
< |
|
631 |
< |
de = f->GetParameter(1); |
632 |
< |
errde = f->GetParError(1); |
633 |
< |
sgmde = f->GetParameter(2); |
634 |
< |
errsgmde = f->GetParError(2); |
635 |
< |
} |
636 |
< |
|
637 |
< |
float rsp = pt/(pt+de); |
638 |
< |
float errrsp = pt/(pt+de)/(pt+de)*errde; |
639 |
< |
float sgm = pt/(pt+de)/(pt+de)*sgmde; |
640 |
< |
float errsgm = 2*pt*sgmde/(pt+de)/(pt+de)/(pt+de)*errde; |
641 |
< |
|
642 |
< |
sgm /= rsp; |
643 |
< |
errsgm /= rsp; |
644 |
< |
|
645 |
< |
gRsp->SetPoint(ibin,var,rsp); |
646 |
< |
gRsp->SetPointError(ibin,0.0,errrsp); |
647 |
< |
gSgm->SetPoint(ibin,var,sgm); |
648 |
< |
gSgm->SetPointError(ibin,0.0,errsgm); |
613 |
> |
TH1F* h(0); |
614 |
> |
|
615 |
> |
if (tVar[ibin]->GetEntries()==0) continue; |
616 |
> |
|
617 |
> |
tVar[ibin]->Project("h(50)","val","weight*(1)"); |
618 |
> |
h = (TH1F*)gROOT->FindObject("h"); |
619 |
> |
h->SetNameTitle(hVar[ibin]->GetName(),hVar[ibin]->GetTitle()); |
620 |
> |
h->SetXTitle(hVar[ibin]->GetXaxis()->GetTitle()); |
621 |
> |
h->SetYTitle(hVar[ibin]->GetYaxis()->GetTitle()); |
622 |
> |
delete hVar[ibin]; |
623 |
> |
hVar[ibin] = h; |
624 |
> |
|
625 |
> |
tVar[ibin]->Project("h(50)","etcorr","weight*(1)"); |
626 |
> |
h = (TH1F*)gROOT->FindObject("h"); |
627 |
> |
h->SetNameTitle(hEtCorr[ibin]->GetName(),hEtCorr[ibin]->GetTitle()); |
628 |
> |
h->SetXTitle(hEtCorr[ibin]->GetXaxis()->GetTitle()); |
629 |
> |
h->SetYTitle(hEtCorr[ibin]->GetYaxis()->GetTitle()); |
630 |
> |
delete hEtCorr[ibin]; |
631 |
> |
hEtCorr[ibin] = h; |
632 |
> |
|
633 |
> |
tVar[ibin]->Project("h(50)","absrsp","weight*(1)"); |
634 |
> |
h = (TH1F*)gROOT->FindObject("h"); |
635 |
> |
h->SetNameTitle(hAbsRsp[ibin]->GetName(),hAbsRsp[ibin]->GetTitle()); |
636 |
> |
h->SetXTitle(hAbsRsp[ibin]->GetXaxis()->GetTitle()); |
637 |
> |
h->SetYTitle(hAbsRsp[ibin]->GetYaxis()->GetTitle()); |
638 |
> |
delete hAbsRsp[ibin]; |
639 |
> |
hAbsRsp[ibin] = h; |
640 |
> |
|
641 |
> |
tVar[ibin]->Project("h(50)","relrsp","weight*(1)"); |
642 |
> |
h = (TH1F*)gROOT->FindObject("h"); |
643 |
> |
h->SetNameTitle(hRelRsp[ibin]->GetName(),hRelRsp[ibin]->GetTitle()); |
644 |
> |
h->SetXTitle(hRelRsp[ibin]->GetXaxis()->GetTitle()); |
645 |
> |
h->SetYTitle(hRelRsp[ibin]->GetYaxis()->GetTitle()); |
646 |
> |
delete hRelRsp[ibin]; |
647 |
> |
hRelRsp[ibin] = h; |
648 |
> |
|
649 |
> |
delete tVar[ibin]; |
650 |
|
} |
651 |
|
|
652 |
< |
gRsp->Write(); |
522 |
< |
gSgm->Write(); |
652 |
> |
delete [] tVar; |
653 |
|
} |
654 |
|
|
655 |
|
|
656 |
|
//______________________________________________________________________________ |
657 |
< |
TCanvas* plotVarHistos(const vector<float>& bins,TH1F** hVar,bool logy) |
657 |
> |
void replaceHistos(int nbins1,int nbins2, |
658 |
> |
TH1F***& hVar,TH1F***& hEtCorr,TH1F***& hAbsRsp,TH1F***& hRelRsp, |
659 |
> |
TTree***& tVar) |
660 |
|
{ |
661 |
< |
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(); |
661 |
> |
if (nbins1==0) return; |
662 |
|
|
663 |
< |
return c; |
663 |
> |
for (int i1=0;i1<nbins1;i1++) |
664 |
> |
replaceHistos(nbins2,hVar[i1],hEtCorr[i1],hAbsRsp[i1],hRelRsp[i1],tVar[i1]); |
665 |
> |
delete [] tVar; |
666 |
|
} |
667 |
|
|
668 |
|
|
669 |
|
//______________________________________________________________________________ |
670 |
< |
TCanvas* plotRspHistos(const vector<float>& bins,TH1F** hRsp) |
670 |
> |
void replaceHistos(int nbins1,int nbins2,int nbins3, |
671 |
> |
TH1F****& hVar,TH1F****& hEtCorr,TH1F****&hAbsRsp,TH1F****&hRelRsp, |
672 |
> |
TTree****& tVar) |
673 |
|
{ |
674 |
< |
string::size_type pos; |
675 |
< |
string cname,jetalg; |
676 |
< |
|
677 |
< |
|
678 |
< |
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; |
674 |
> |
if (nbins1==0) return; |
675 |
> |
for (int i1=0;i1<nbins1;i1++) |
676 |
> |
replaceHistos(nbins2,nbins3, |
677 |
> |
hVar[i1],hEtCorr[i1],hAbsRsp[i1],hRelRsp[i1],tVar[i1]); |
678 |
> |
delete [] tVar; |
679 |
|
} |
680 |
|
|
681 |
|
|
682 |
|
//______________________________________________________________________________ |
683 |
< |
TCanvas* plotResponse(const vector<string>& jetalgs,short applyjes, |
606 |
< |
const vector<TGraphErrors*> gRsp, |
607 |
< |
const string& varname,const string& xtitle) |
683 |
> |
void fitHistos(TH1F** histos,int n,double nsigma) |
684 |
|
{ |
685 |
< |
assert(jetalgs.size()==gRsp.size()); |
686 |
< |
|
687 |
< |
string cname = "RspVs"+varname; |
688 |
< |
TCanvas* c = new TCanvas(cname.c_str(),cname.c_str()); |
689 |
< |
TLegend* l = new TLegend(0.7,0.25+jetalgs.size()*0.045,0.85,0.25); |
690 |
< |
|
691 |
< |
l->SetFillColor(10); |
692 |
< |
Color_t color=kBlack; |
693 |
< |
TMultiGraph* mg = new TMultiGraph(); |
694 |
< |
|
695 |
< |
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++; |
685 |
> |
for (int i=0;i<n;i++) { |
686 |
> |
float nrm = histos[i]->Integral(); |
687 |
> |
float mean = histos[i]->GetMean(); |
688 |
> |
float rms = histos[i]->GetRMS(); |
689 |
> |
TF1* f = new TF1("gaus","gaus",mean-nsigma*rms,mean+nsigma*rms); |
690 |
> |
f->SetParameter(0,nrm); |
691 |
> |
f->SetParameter(1,mean); |
692 |
> |
f->SetParameter(2,rms); |
693 |
> |
f->SetLineWidth(1); |
694 |
> |
f->SetLineColor(kRed); |
695 |
> |
histos[i]->Fit(f,"QR"); // Avoid TCanvas? |
696 |
|
} |
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; |
697 |
|
} |
698 |
|
|
699 |
|
|
700 |
|
//______________________________________________________________________________ |
701 |
< |
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 |
< |
|
701 |
> |
void fitHistos(TH1F*** histos,int n1,int n2,double nsigma) |
702 |
|
{ |
703 |
< |
int nbins = bins1.size()+1; |
704 |
< |
|
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; |
703 |
> |
if (n1==0||n2==0) return; |
704 |
> |
for (int i1=0;i1<n1;i1++) fitHistos(histos[i1],n2,nsigma); |
705 |
|
} |
706 |
|
|
707 |
|
|
708 |
|
//______________________________________________________________________________ |
709 |
< |
TCanvas* plotResolution(const vector<string>& jetalgs, |
772 |
< |
const vector<TGraphErrors*> gSgm, |
773 |
< |
const string& varname,const string& xtitle) |
709 |
> |
void fitHistos(TH1F**** histos,int n1,int n2,int n3,double nsigma) |
710 |
|
{ |
711 |
< |
string cname = "SgmVs"+varname; |
712 |
< |
TCanvas* c = new TCanvas(cname.c_str(),cname.c_str()); |
713 |
< |
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 |
< |
} |
711 |
> |
if (n1==0||n2==0||n3==0) return; |
712 |
> |
for (int i1=0;i1<n1;i1++) fitHistos(histos[i1],n2,n3,nsigma); |
713 |
> |
} |
714 |
|
|
715 |
< |
g->Write(); |
716 |
< |
|
717 |
< |
color++; |
715 |
> |
|
716 |
> |
//______________________________________________________________________________ |
717 |
> |
void saveHistos(TH1F**& histos,int n) |
718 |
> |
{ |
719 |
> |
if (n==0) return; |
720 |
> |
for (int i=0;i<n;i++) { |
721 |
> |
histos[i]->Write(); |
722 |
> |
delete histos[i]; |
723 |
|
} |
724 |
< |
|
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; |
724 |
> |
delete [] histos; |
725 |
|
} |
726 |
|
|
727 |
|
|
728 |
|
//______________________________________________________________________________ |
729 |
< |
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) |
729 |
> |
void saveHistos(TH1F***& histos,int n1,int n2) |
730 |
|
{ |
731 |
< |
int nbins = bins1.size()+1; |
732 |
< |
|
733 |
< |
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; |
731 |
> |
if (n1==0) return; |
732 |
> |
for (int i1=0;i1<n1;i1++) saveHistos(histos[i1],n2); |
733 |
> |
delete [] histos; |
734 |
|
} |
735 |
|
|
736 |
|
|
737 |
|
//______________________________________________________________________________ |
738 |
< |
float get_xmax(TGraph* g) |
738 |
> |
void saveHistos(TH1F****& histos,int n1,int n2,int n3) |
739 |
|
{ |
740 |
< |
float result(-1e200); |
741 |
< |
for (int i=0;i<g->GetN();++i) if (g->GetX()[i]>result) result=g->GetX()[i]; |
742 |
< |
return result; |
740 |
> |
if (n1==0) return; |
741 |
> |
for (int i1=0;i1<n1;i1++) saveHistos(histos[i1],n2,n3); |
742 |
> |
delete [] histos; |
743 |
|
} |