ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/SchieferD/jetcalib/mcresponse_x.cpp
Revision: 1.14
Committed: Tue Sep 4 12:10:40 2007 UTC (17 years, 7 months ago) by schiefer
Branch: MAIN
Changes since 1.13: +18 -18 lines
Log Message:
*** empty log message ***

File Contents

# User Rev Content
1 schiefer 1.1 ////////////////////////////////////////////////////////////////////////////////
2     //
3     // mcresponse_x
4     // ------------
5     //
6     // 06/16/2007 Philipp Schieferdecker <philipp.schieferdecker@cern.ch>
7     ////////////////////////////////////////////////////////////////////////////////
8    
9    
10     #include "utils/cmdline.h"
11 schiefer 1.2 #include "utils/hist.h"
12 schiefer 1.1
13 schiefer 1.9 #include <TSystem.h>
14 schiefer 1.1 #include <TROOT.h>
15     #include <TRint.h>
16     #include <TFile.h>
17     #include <TTree.h>
18     #include <TEventList.h>
19     #include <TH1F.h>
20     #include <TF1.h>
21 schiefer 1.6 #include <TMath.h>
22 schiefer 1.1
23     #include <iostream>
24     #include <iomanip>
25     #include <cmath>
26     #include <string>
27     #include <vector>
28    
29    
30     using namespace std;
31    
32    
33     ////////////////////////////////////////////////////////////////////////////////
34     // declare global functions
35     ////////////////////////////////////////////////////////////////////////////////
36 schiefer 1.4 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 schiefer 1.5 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 schiefer 1.1
53    
54    
55     ////////////////////////////////////////////////////////////////////////////////
56     // main
57     ////////////////////////////////////////////////////////////////////////////////
58    
59     //______________________________________________________________________________
60     int main(int argc,char**argv)
61     {
62     cmdline cl;
63     cl.parse(argc,argv);
64 schiefer 1.2
65 schiefer 1.1 vector<string> input = cl.get_vector<string>("input");
66 schiefer 1.4 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 schiefer 1.13 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 schiefer 1.4 short applyjes = cl.get_value<short> ("applyjes", 0);
74     string treename = cl.get_value<string> ("treename", "t");
75     float nsigma = cl.get_value<float> ("nsigma", 3.0);
76 schiefer 1.7 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 schiefer 1.4 bool abseta = cl.get_value<bool> ("abseta", true);
81 schiefer 1.1
82     if (!cl.check()) return 0;
83     cl.print();
84    
85 schiefer 1.7 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 schiefer 1.1 // etabins
96 schiefer 1.7 if (netabins>0&&!abseta) {
97 schiefer 1.1 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]);
100     for (int ieta=0;ieta<neta;ieta++) etabins[ieta]*=-1;
101     }
102    
103    
104     // jetalgs
105     vector<string> jetalgs;
106     for (vector<string>::const_iterator it=input.begin();it!=input.end();++it) {
107     string jetalg = *it;
108     string::size_type pos;
109     while ((pos=jetalg.find('/'))!=string::npos) jetalg = jetalg.substr(pos+1);
110     jetalg = jetalg.substr(0,jetalg.find(".root"));
111     jetalgs.push_back(jetalg);
112     }
113    
114     // prepare branch variables and histograms / graphs
115     float weight;
116     char njt;
117 schiefer 1.3 float jtet[100];
118 schiefer 1.1 float jteta[100];
119     float jtphi[100];
120 schiefer 1.4 float jtemf[100];
121 schiefer 1.3 float jtgenet[100];
122 schiefer 1.1 float jtgendr[100];
123     float jtjes[100][3];
124    
125 schiefer 1.3 argc=3; argv[1]="-l"; argv[2]="-b";
126 schiefer 1.4 TRint* app = new TRint(argv[0],&argc,argv); app->Argc();
127 schiefer 1.1
128    
129     //
130     // loop over all input files
131     //
132 schiefer 1.6 TFile* plotfile = new TFile(output.c_str(),"UPDATE");
133 schiefer 1.2
134 schiefer 1.1 for (unsigned int i=0;i<input.size();++i) {
135 schiefer 1.9 TFile* f = new TFile(input[i].c_str(),"READ"); if (!f->IsOpen()) return 0;
136 schiefer 1.1 TTree* t = (TTree*)f->Get(treename.c_str());
137    
138     TEventList* el = new TEventList("sel","sel");
139     t->Draw(">>sel",selection.c_str());
140    
141     t->SetBranchAddress("weight", &weight);
142     t->SetBranchAddress("njt", &njt);
143 schiefer 1.3 t->SetBranchAddress("jtet", jtet);
144 schiefer 1.1 t->SetBranchAddress("jteta", jteta);
145     t->SetBranchAddress("jtphi", jtphi);
146 schiefer 1.4 t->SetBranchAddress("jtemf", jtemf);
147 schiefer 1.3 t->SetBranchAddress("jtgenet",jtgenet);
148 schiefer 1.1 t->SetBranchAddress("jtgendr",jtgendr);
149 schiefer 1.10
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 schiefer 1.1 for (int i1=0;i1<100;i1++)
166     for (int i2=0;i2<3;i2++) jtjes[i1][i2]=1.0;
167 schiefer 1.10 }
168 schiefer 1.3 string absrsp_xtitle = "E_{T}^{gen}-E_{T} [GeV]";
169 schiefer 1.14 string relrsp_xtitle = "E_{T}/E_{T}^{gen}";
170 schiefer 1.5
171 schiefer 1.4
172 schiefer 1.7 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 schiefer 1.14 "E_{T}^{gen} [GeV]",
180 schiefer 1.7 "GenEt",etbins);
181     hEtCorrGenEt = hist::initHistos("EtCorrGenEt",jetalgs[i],100,
182 schiefer 1.14 "E_{T}^{gen} [GeV]",
183 schiefer 1.7 "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 schiefer 1.14 "E_{T} [GeV]","Et",etbins);
203 schiefer 1.7 hEtCorrEt = hist::initHistos("EtCorrEt",jetalgs[i],100,
204 schiefer 1.14 "E_{T} [GeV]","Et",etbins);
205 schiefer 1.7 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 schiefer 1.14 "#eta","Eta",etabins);
220 schiefer 1.7 hEtCorrEta = hist::initHistos("EtCorrEta",jetalgs[i],100,
221 schiefer 1.14 "E_{T} [GeV]","Eta",etabins);
222 schiefer 1.7 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 schiefer 1.14 "#phi","Phi",phibins);
237 schiefer 1.7 hEtCorrPhi = hist::initHistos("EtCorrPhi",jetalgs[i],100,
238 schiefer 1.14 "E_{T} [Gev]","Phi",phibins);
239 schiefer 1.7 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 schiefer 1.14 "emf","Emf",emfbins);
254 schiefer 1.7 hEtCorrEmf = hist::initHistos("EtCorrEmf",jetalgs[i],100,
255 schiefer 1.14 "E_{T} [Gev]","Emf",emfbins);
256 schiefer 1.7 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 schiefer 1.14 hEtEtaEt = hist::initHistos("Et",jetalgs[i],100,"E_{T} [GeV]",
270 schiefer 1.7 "Eta",etabins,"Et",etbins);
271     hEtCorrEtaEt = hist::initHistos("EtCorrEtaEt",jetalgs[i],100,
272 schiefer 1.14 "E_{T} [GeV]",
273 schiefer 1.7 "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 schiefer 1.14 hEtEmfEt = hist::initHistos("Et",jetalgs[i],100,"E_{T} [GeV]",
290 schiefer 1.7 "Emf",emfbins,"Et",etbins);
291     hEtCorrEmfEt = hist::initHistos("EtCorrEmfEt",jetalgs[i],100,
292 schiefer 1.14 "E_{T} [GeV]",
293 schiefer 1.7 "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 schiefer 1.14 hEtEmfEtaEt = hist::initHistos("Et",jetalgs[i],100,"E_{T} [GeV]",
310 schiefer 1.7 "Emf",emfbins,
311     "Eta",etabins,
312     "Et", etbins);
313     hEtCorrEmfEtaEt = hist::initHistos("EtCorrEmfEtaEt",jetalgs[i],100,
314 schiefer 1.14 "E_{T} [GeV]",
315 schiefer 1.7 "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 schiefer 1.1
330     int nevts = el->GetN();
331     cout<<jetalgs[i]<<": "<<nevts<<" events selected."<<endl;
332 schiefer 1.4
333 schiefer 1.3 // prepare intermediate trees
334 schiefer 1.4 float genet,et,eta,phi,emf,etcorr,absrsp,relrsp;
335 schiefer 1.7
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 schiefer 1.10 gROOT->cd();
346    
347 schiefer 1.7 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 schiefer 1.4
356 schiefer 1.3 for (int iet=0;iet<netbins;iet++) {
357 schiefer 1.4
358 schiefer 1.3 tGenEt[iet]=new TTree("tGenEt","tGenEt");
359     tGenEt[iet]->Branch("val", &genet, "val/F");
360 schiefer 1.14 tGenEt[iet]->Branch("etcorr",&genet, "etcorr/F");
361 schiefer 1.4 tGenEt[iet]->Branch("absrsp",&absrsp,"absrsp/F");
362     tGenEt[iet]->Branch("relrsp",&relrsp,"relrsp/F");
363 schiefer 1.3 tGenEt[iet]->Branch("weight",&weight,"weight/F");
364 schiefer 1.4
365 schiefer 1.3 tEt[iet]=new TTree("tEt","tEt");
366     tEt[iet]->Branch("val", &et, "val/F");
367     tEt[iet]->Branch("etcorr",&etcorr,"etcorr/F");
368 schiefer 1.4 tEt[iet]->Branch("absrsp",&absrsp,"absrsp/F");
369     tEt[iet]->Branch("relrsp",&relrsp,"relrsp/F");
370 schiefer 1.3 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 schiefer 1.4 tEta[ieta]->Branch("absrsp",&absrsp,"absrsp/F");
377     tEta[ieta]->Branch("relrsp",&relrsp,"relrsp/F");
378 schiefer 1.3 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 schiefer 1.4 tPhi[iphi]->Branch("absrsp",&absrsp,"absrsp/F");
385     tPhi[iphi]->Branch("relrsp",&relrsp,"relrsp/F");
386 schiefer 1.3 tPhi[iphi]->Branch("weight",&weight,"weight/F");
387     }
388 schiefer 1.4 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 schiefer 1.8
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 schiefer 1.3 }
409 schiefer 1.8 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 schiefer 1.4 }
420 schiefer 1.8 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 schiefer 1.4 }
436     }
437     }
438 schiefer 1.1
439 schiefer 1.5 cout<<jetalgs[i]<<": trees created."<<endl;
440    
441 schiefer 1.1 // loop over all events
442     for (int ievt=0;ievt<nevts;ievt++) {
443    
444 schiefer 1.10 //t->GetEntry(el->GetEntry(ievt));
445     for (int ibrnch=0;ibrnch<(int)branches.size();ibrnch++)
446     branches[ibrnch]->GetEntry(el->GetEntry(ievt));
447    
448 schiefer 1.1 // loop over all jets in the event
449     for (int ijt=0;ijt<njt;ijt++) {
450    
451 schiefer 1.13 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 schiefer 1.1
455 schiefer 1.3 et = jtet[ijt];
456     genet = jtgenet[ijt];
457     eta = (abseta) ? std::abs(jteta[ijt]) : jteta[ijt];
458     phi = jtphi[ijt];
459 schiefer 1.4 emf = jtemf[ijt];
460 schiefer 1.6
461 schiefer 1.3 etcorr = (applyjes==0) ? et : et*jtjes[ijt][applyjes-1];
462 schiefer 1.1
463 schiefer 1.4 absrsp = genet - etcorr;
464     relrsp = etcorr/genet;
465 schiefer 1.2
466 schiefer 1.6 if (TMath::IsNaN(emf)) emf = 0.0;
467    
468 schiefer 1.7 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     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 schiefer 1.4
486 schiefer 1.1 } // jets
487     } // evts
488    
489     cout<<jetalgs[i]<<": all events processed / histograms filled."<<endl;
490 schiefer 1.3
491 schiefer 1.2
492 schiefer 1.3 // replace histograms
493 schiefer 1.4
494     replaceHistos(netbins,
495 schiefer 1.5 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 schiefer 1.4 replaceHistos(netabins,netbins,
501 schiefer 1.5 hEtEtaEt,hEtCorrEtaEt,hAbsRspEtaEt,hRelRspEtaEt,
502 schiefer 1.4 tEtEtaEt);
503     replaceHistos(nemfbins,netbins,
504 schiefer 1.5 hEtEmfEt,hEtCorrEmfEt,hAbsRspEmfEt,hRelRspEmfEt,
505 schiefer 1.4 tEtEmfEt);
506     replaceHistos(nemfbins,netabins,netbins,
507 schiefer 1.5 hEtEmfEtaEt,hEtCorrEmfEtaEt,
508     hAbsRspEmfEtaEt,hRelRspEmfEtaEt,
509 schiefer 1.4 tEtEmfEtaEt);
510 schiefer 1.2
511 schiefer 1.5 cout<<jetalgs[i]<<": all histograms replaced, and tree deleted."<<endl;
512    
513 schiefer 1.6
514 schiefer 1.9 plotfile->cd();
515 schiefer 1.6 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 schiefer 1.2 d->cd();
519 schiefer 1.1
520 schiefer 1.3 // fit response histos
521 schiefer 1.5 fitHistos(hAbsRspGenEt,netbins,nsigma);
522     fitHistos(hRelRspGenEt,netbins,nsigma);
523 schiefer 1.9
524 schiefer 1.5 fitHistos(hAbsRspEt, netbins,nsigma);
525     fitHistos(hRelRspEt, netbins,nsigma);
526 schiefer 1.3
527 schiefer 1.5 fitHistos(hAbsRspEta, netabins,nsigma);
528     fitHistos(hRelRspEta, netabins,nsigma);
529 schiefer 1.3
530 schiefer 1.5 fitHistos(hAbsRspPhi, nphibins,nsigma);
531     fitHistos(hRelRspPhi, nphibins,nsigma);
532 schiefer 1.9
533 schiefer 1.5 fitHistos(hAbsRspEmf, nemfbins,nsigma);
534     fitHistos(hRelRspEmf, nemfbins,nsigma);
535 schiefer 1.9
536 schiefer 1.5 fitHistos(hAbsRspEtaEt,netabins,netbins,nsigma);
537     fitHistos(hRelRspEtaEt,netabins,netbins,nsigma);
538 schiefer 1.3
539 schiefer 1.5 fitHistos(hAbsRspEmfEt,nemfbins,netbins,nsigma);
540     fitHistos(hRelRspEmfEt,nemfbins,netbins,nsigma);
541 schiefer 1.4
542 schiefer 1.5 fitHistos(hAbsRspEmfEtaEt,nemfbins,netabins,netbins,nsigma);
543     fitHistos(hRelRspEmfEtaEt,nemfbins,netabins,netbins,nsigma);
544    
545     cout<<jetalgs[i]<<": fits done."<<endl;
546 schiefer 1.4
547 schiefer 1.3
548     // save histograms
549 schiefer 1.5 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 schiefer 1.9 cout<<jetalgs[i]<<": histograms saved."<<endl;
590 schiefer 1.1 }
591    
592 schiefer 1.2 plotfile->Write();
593     plotfile->Close();
594     delete plotfile;
595 schiefer 1.9
596 schiefer 1.12 //app->Run();
597 schiefer 1.1
598     return 0;
599     }
600    
601    
602    
603     ////////////////////////////////////////////////////////////////////////////////
604     // implementation of global functions
605     ////////////////////////////////////////////////////////////////////////////////
606    
607     //______________________________________________________________________________
608 schiefer 1.4 void replaceHistos(int nbins,
609     TH1F**& hVar,TH1F**& hEtCorr,TH1F**& hAbsRsp,TH1F**& hRelRsp,
610     TTree**& tVar)
611 schiefer 1.1 {
612 schiefer 1.7 if (nbins==0) return;
613    
614 schiefer 1.3 TH1::SetDefaultSumw2();
615 schiefer 1.1 for (int ibin=0;ibin<nbins;ibin++) {
616 schiefer 1.4
617 schiefer 1.3 TH1F* h(0);
618 schiefer 1.4
619     if (tVar[ibin]->GetEntries()==0) continue;
620 schiefer 1.1
621 schiefer 1.3 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     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 schiefer 1.4
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     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 schiefer 1.3 delete tVar[ibin];
654 schiefer 1.1 }
655    
656 schiefer 1.3 delete [] tVar;
657 schiefer 1.1 }
658    
659 schiefer 1.4
660     //______________________________________________________________________________
661     void replaceHistos(int nbins1,int nbins2,
662     TH1F***& hVar,TH1F***& hEtCorr,TH1F***& hAbsRsp,TH1F***& hRelRsp,
663     TTree***& tVar)
664     {
665 schiefer 1.7 if (nbins1==0) return;
666    
667 schiefer 1.4 for (int i1=0;i1<nbins1;i1++)
668     replaceHistos(nbins2,hVar[i1],hEtCorr[i1],hAbsRsp[i1],hRelRsp[i1],tVar[i1]);
669 schiefer 1.5 delete [] tVar;
670 schiefer 1.4 }
671    
672    
673     //______________________________________________________________________________
674     void replaceHistos(int nbins1,int nbins2,int nbins3,
675     TH1F****& hVar,TH1F****& hEtCorr,TH1F****&hAbsRsp,TH1F****&hRelRsp,
676     TTree****& tVar)
677     {
678 schiefer 1.7 if (nbins1==0) return;
679 schiefer 1.4 for (int i1=0;i1<nbins1;i1++)
680     replaceHistos(nbins2,nbins3,
681     hVar[i1],hEtCorr[i1],hAbsRsp[i1],hRelRsp[i1],tVar[i1]);
682 schiefer 1.5 delete [] tVar;
683 schiefer 1.4 }
684    
685    
686 schiefer 1.1 //______________________________________________________________________________
687 schiefer 1.3 void fitHistos(TH1F** histos,int n,double nsigma)
688 schiefer 1.1 {
689 schiefer 1.3 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 schiefer 1.10 histos[i]->Fit(f,"QR"); // Avoid TCanvas?
700 schiefer 1.1 }
701     }
702    
703    
704     //______________________________________________________________________________
705 schiefer 1.3 void fitHistos(TH1F*** histos,int n1,int n2,double nsigma)
706 schiefer 1.1 {
707 schiefer 1.9 if (n1==0||n2==0) return;
708 schiefer 1.3 for (int i1=0;i1<n1;i1++) fitHistos(histos[i1],n2,nsigma);
709 schiefer 1.1 }
710    
711    
712     //______________________________________________________________________________
713 schiefer 1.4 void fitHistos(TH1F**** histos,int n1,int n2,int n3,double nsigma)
714     {
715 schiefer 1.9 if (n1==0||n2==0||n3==0) return;
716 schiefer 1.4 for (int i1=0;i1<n1;i1++) fitHistos(histos[i1],n2,n3,nsigma);
717     }
718    
719    
720     //______________________________________________________________________________
721 schiefer 1.5 void saveHistos(TH1F**& histos,int n)
722 schiefer 1.1 {
723 schiefer 1.7 if (n==0) return;
724 schiefer 1.5 for (int i=0;i<n;i++) {
725     histos[i]->Write();
726     delete histos[i];
727     }
728     delete [] histos;
729 schiefer 1.1 }
730    
731    
732     //______________________________________________________________________________
733 schiefer 1.5 void saveHistos(TH1F***& histos,int n1,int n2)
734 schiefer 1.1 {
735 schiefer 1.7 if (n1==0) return;
736 schiefer 1.3 for (int i1=0;i1<n1;i1++) saveHistos(histos[i1],n2);
737 schiefer 1.5 delete [] histos;
738 schiefer 1.1 }
739 schiefer 1.4
740    
741     //______________________________________________________________________________
742 schiefer 1.5 void saveHistos(TH1F****& histos,int n1,int n2,int n3)
743 schiefer 1.4 {
744 schiefer 1.7 if (n1==0) return;
745 schiefer 1.4 for (int i1=0;i1<n1;i1++) saveHistos(histos[i1],n2,n3);
746 schiefer 1.5 delete [] histos;
747 schiefer 1.4 }