ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/SchieferD/jetcalib/mcresponse_x.cpp
Revision: 1.4
Committed: Thu Aug 30 09:01:17 2007 UTC (17 years, 8 months ago) by schiefer
Branch: MAIN
Changes since 1.3: +331 -136 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     #include <TROOT.h>
14     #include <TRint.h>
15     #include <TFile.h>
16     #include <TTree.h>
17     #include <TEventList.h>
18     #include <TH1F.h>
19     #include <TF1.h>
20    
21     #include <iostream>
22     #include <iomanip>
23     #include <cmath>
24     #include <string>
25     #include <vector>
26    
27    
28     using namespace std;
29    
30    
31     ////////////////////////////////////////////////////////////////////////////////
32     // defaults
33     ////////////////////////////////////////////////////////////////////////////////
34 schiefer 1.4 string dflt_etbins ="8,12,15,20,25,30,45,65,90,120,180";
35 schiefer 1.3 string dflt_etabins=".2,.4,.6,.8,1.,1.2,1.4,1.6,1.8,2.,2.3,2.5,2.7,3.0,3.3";
36 schiefer 1.1 string dflt_phibins="-2.75,-2.25,-1.75.,-1.25,-.75,-.25,.25,.75,1.25,1.75,2.25,2.75";
37 schiefer 1.4 string dflt_emfbins="0.2,0.4,0.6,0.8";
38 schiefer 1.1
39    
40     ////////////////////////////////////////////////////////////////////////////////
41     // declare global functions
42     ////////////////////////////////////////////////////////////////////////////////
43 schiefer 1.4 void replaceHistos(int nbins,
44     TH1F**& hVar,TH1F**& hEtCorr,TH1F**& hAbsRsp,TH1F**& hRelRsp,
45     TTree**& tVar);
46     void replaceHistos(int nbins1,int nbins2,
47     TH1F***&hVar,TH1F***&hEtCorr,TH1F***&hAbsRsp,TH1F***&hRelRsp,
48     TTree***& tVar);
49     void replaceHistos(int nbins1,int nbins2,int nbins3,
50     TH1F****&hVar,TH1F****&hEtCorr,
51     TH1F****&hAbsRsp,TH1F****&hRelRsp,
52     TTree****& tVar);
53     void fitHistos (TH1F** histos,int n,double nsigma);
54     void fitHistos (TH1F*** histos,int n1,int n2,double nsigma);
55     void fitHistos (TH1F**** histos,int n1,int n2,int n3,double nsigma);
56     void saveHistos (TH1F** histos,int n);
57     void saveHistos (TH1F*** histos,int n1,int n2);
58     void saveHistos (TH1F**** histos,int n1,int n2,int n3);
59 schiefer 1.1
60    
61    
62     ////////////////////////////////////////////////////////////////////////////////
63     // main
64     ////////////////////////////////////////////////////////////////////////////////
65    
66     //______________________________________________________________________________
67     int main(int argc,char**argv)
68     {
69     cmdline cl;
70     cl.parse(argc,argv);
71 schiefer 1.2
72 schiefer 1.1 vector<string> input = cl.get_vector<string>("input");
73 schiefer 1.4 string output = cl.get_value<string> ("output", "mcrsp.root");
74     string selection = cl.get_value<string> ("selection", "njt>0");
75     float drmax = cl.get_value<float> ("drmax", 0.3);
76     short applyjes = cl.get_value<short> ("applyjes", 0);
77     string treename = cl.get_value<string> ("treename", "t");
78     float nsigma = cl.get_value<float> ("nsigma", 3.0);
79     vector<float> etbins = cl.get_vector<float> ("etbins", dflt_etbins);
80     vector<float> etabins = cl.get_vector<float> ("etabins",dflt_etabins);
81     vector<float> phibins = cl.get_vector<float> ("phibins",dflt_phibins);
82     vector<float> emfbins = cl.get_vector<float> ("emfbins",dflt_emfbins);
83     bool abseta = cl.get_value<bool> ("abseta", true);
84 schiefer 1.1
85     if (!cl.check()) return 0;
86     cl.print();
87    
88    
89     // etabins
90     if (!abseta) {
91     int neta=(int)etabins.size();
92     std::reverse(etabins.begin(),etabins.end());
93     for (int ieta=neta-1;ieta>=0;ieta--) etabins.push_back(etabins[ieta]);
94     for (int ieta=0;ieta<neta;ieta++) etabins[ieta]*=-1;
95     }
96    
97    
98     // jetalgs
99     vector<string> jetalgs;
100     for (vector<string>::const_iterator it=input.begin();it!=input.end();++it) {
101     string jetalg = *it;
102     string::size_type pos;
103     while ((pos=jetalg.find('/'))!=string::npos) jetalg = jetalg.substr(pos+1);
104     jetalg = jetalg.substr(0,jetalg.find(".root"));
105     jetalgs.push_back(jetalg);
106     }
107    
108     // prepare branch variables and histograms / graphs
109     float weight;
110     char njt;
111 schiefer 1.3 float jtet[100];
112 schiefer 1.1 float jtpt[100];
113     float jteta[100];
114     float jtphi[100];
115 schiefer 1.4 float jtemf[100];
116 schiefer 1.3 float jtgenet[100];
117 schiefer 1.1 float jtgenpt[100];
118     float jtgendr[100];
119     float jtjes[100][3];
120    
121 schiefer 1.3 int netbins = etbins.size()+1;
122 schiefer 1.1 int netabins = etabins.size()+1;
123     int nphibins = phibins.size()+1;
124 schiefer 1.4 int nemfbins = emfbins.size()+1;
125 schiefer 1.1
126 schiefer 1.4 vector<TH1F**> hGenEt;
127     vector<TH1F**> hEtCorrGenEt;
128     vector<TH1F**> hAbsRspGenEt;
129     vector<TH1F**> hRelRspGenEt;
130    
131     vector<TH1F**> hEt;
132     vector<TH1F**> hEtCorrEt;
133     vector<TH1F**> hAbsRspEt;
134     vector<TH1F**> hRelRspEt;
135    
136     vector<TH1F**> hEta;
137     vector<TH1F**> hEtCorrEta;
138     vector<TH1F**> hAbsRspEta;
139     vector<TH1F**> hRelRspEta;
140    
141     vector<TH1F**> hPhi;
142     vector<TH1F**> hEtCorrPhi;
143     vector<TH1F**> hAbsRspPhi;
144     vector<TH1F**> hRelRspPhi;
145    
146     vector<TH1F**> hEmf;
147     vector<TH1F**> hEtCorrEmf;
148     vector<TH1F**> hAbsRspEmf;
149     vector<TH1F**> hRelRspEmf;
150    
151     vector<TH1F***> hEtEtaEt;
152     vector<TH1F***> hEtCorrEtaEt;
153     vector<TH1F***> hAbsRspEtaEt;
154     vector<TH1F***> hRelRspEtaEt;
155    
156     vector<TH1F***> hEtEmfEt;
157     vector<TH1F***> hEtCorrEmfEt;
158     vector<TH1F***> hAbsRspEmfEt;
159     vector<TH1F***> hRelRspEmfEt;
160    
161     vector<TH1F****> hEtEmfEtaEt;
162     vector<TH1F****> hEtCorrEmfEtaEt;
163     vector<TH1F****> hAbsRspEmfEtaEt;
164     vector<TH1F****> hRelRspEmfEtaEt;
165 schiefer 1.1
166    
167 schiefer 1.3 argc=3; argv[1]="-l"; argv[2]="-b";
168 schiefer 1.4 TRint* app = new TRint(argv[0],&argc,argv); app->Argc();
169 schiefer 1.1
170    
171     //
172     // loop over all input files
173     //
174 schiefer 1.3 TFile* plotfile = new TFile(output.c_str(),"RECREATE");
175 schiefer 1.2
176 schiefer 1.1 for (unsigned int i=0;i<input.size();++i) {
177 schiefer 1.3 TFile* f = new TFile(input[i].c_str(),"UPDATE"); if (!f->IsOpen()) return 0;
178 schiefer 1.1 TTree* t = (TTree*)f->Get(treename.c_str());
179    
180     TEventList* el = new TEventList("sel","sel");
181     t->Draw(">>sel",selection.c_str());
182    
183     t->SetBranchAddress("weight", &weight);
184     t->SetBranchAddress("njt", &njt);
185 schiefer 1.3 t->SetBranchAddress("jtet", jtet);
186 schiefer 1.1 t->SetBranchAddress("jtpt", jtpt);
187     t->SetBranchAddress("jteta", jteta);
188     t->SetBranchAddress("jtphi", jtphi);
189 schiefer 1.4 t->SetBranchAddress("jtemf", jtemf);
190 schiefer 1.3 t->SetBranchAddress("jtgenet",jtgenet);
191 schiefer 1.1 t->SetBranchAddress("jtgenpt",jtgenpt);
192     t->SetBranchAddress("jtgendr",jtgendr);
193    
194     if (t->FindBranch("jtjes")) t->SetBranchAddress("jtjes", jtjes);
195     else
196     for (int i1=0;i1<100;i1++)
197     for (int i2=0;i2<3;i2++) jtjes[i1][i2]=1.0;
198    
199    
200 schiefer 1.3 string absrsp_xtitle = "E_{T}^{gen}-E_{T} [GeV]";
201     string relrsp_xtitle = "E_{T}^{rec}/E_{T}";
202 schiefer 1.2
203 schiefer 1.4 hGenEt .push_back(hist::initHistos("GenEt",jetalgs[i],100,
204     "jet E_{T}^{gen} [GeV]","GenEt",etbins));
205     hEtCorrGenEt.push_back(hist::initHistos("EtCorrGenEt",jetalgs[i],100,
206     "jet E_{T}^{gen} [GeV]","GenEt",etbins));
207     hAbsRspGenEt.push_back(hist::initHistos("AbsRspGenEt",jetalgs[i],100,-50,150,
208     absrsp_xtitle,"GenEt",etbins));
209     hRelRspGenEt.push_back(hist::initHistos("RelRspGenEt",jetalgs[i],100,0,2,
210     relrsp_xtitle,"GenEt",etbins));
211    
212     hEt .push_back(hist::initHistos("Et",jetalgs[i],100,
213     "jet E_{T} [GeV]","Et",etbins));
214     hEtCorrEt .push_back(hist::initHistos("EtCorrEt",jetalgs[i],100,
215     "jet E_{T} [GeV]","Et",etbins));
216     hAbsRspEt .push_back(hist::initHistos("AbsRspEt",jetalgs[i],100,-50,150,
217     absrsp_xtitle,"Et",etbins));
218     hRelRspEt .push_back(hist::initHistos("RelRspEt",jetalgs[i],100,0,2,
219     relrsp_xtitle,"Et",etbins));
220    
221     hEta .push_back(hist::initHistos("Eta",jetalgs[i],100,
222     "jet #eta","Eta",etabins));
223     hEtCorrEta .push_back(hist::initHistos("EtCorrEta",jetalgs[i],100,
224     "jet E_{T} [GeV]","Eta",etabins));
225     hAbsRspEta .push_back(hist::initHistos("AbsRspEta",jetalgs[i],100,-50,150,
226     absrsp_xtitle,"Eta",etabins));
227     hRelRspEta .push_back(hist::initHistos("RelRspEta",jetalgs[i],100,0,2,
228     relrsp_xtitle,"Eta",etabins));
229    
230     hPhi .push_back(hist::initHistos("Phi",jetalgs[i],100,
231     "jet #phi","Phi",phibins));
232     hEtCorrPhi .push_back(hist::initHistos("EtCorrPhi",jetalgs[i],100,
233     "jet E_{T} [Gev]","Phi",phibins));
234     hAbsRspPhi .push_back(hist::initHistos("AbsRspPhi",jetalgs[i],100,-50,150,
235     absrsp_xtitle,"Phi",phibins));
236     hRelRspPhi .push_back(hist::initHistos("RelRspPhi",jetalgs[i],100,0,2,
237     relrsp_xtitle,"Phi",phibins));
238    
239     hEmf .push_back(hist::initHistos("Emf",jetalgs[i],100,
240     "jet emf","Emf",emfbins));
241     hEtCorrEmf .push_back(hist::initHistos("EtCorrEmf",jetalgs[i],100,
242     "jet E_{T} [Gev]","Emf",emfbins));
243     hAbsRspEmf .push_back(hist::initHistos("AbsRspEmf",jetalgs[i],100,-50,150,
244     absrsp_xtitle,"Emf",emfbins));
245     hRelRspEmf .push_back(hist::initHistos("RelRspEmf",jetalgs[i],100,0,2,
246     relrsp_xtitle,"Emf",emfbins));
247 schiefer 1.2
248 schiefer 1.4 hEtEtaEt .push_back(hist::initHistos("Et",jetalgs[i],100,"jet E_{T} [GeV]",
249 schiefer 1.3 "Eta",etabins,"Et",etbins));
250 schiefer 1.4 hEtCorrEtaEt.push_back(hist::initHistos("EtCorrEtaEt",jetalgs[i],100,
251     "jet E_{T} [GeV]",
252 schiefer 1.3 "Eta",etabins,"Et",etbins));
253 schiefer 1.4 hAbsRspEtaEt.push_back(hist::initHistos("AbsRspEtaEt",jetalgs[i],100,-50,150,
254     absrsp_xtitle,
255 schiefer 1.3 "Eta",etabins,"Et",etbins));
256 schiefer 1.4 hRelRspEtaEt.push_back(hist::initHistos("RelRspEtaEt",jetalgs[i],100,0,2,
257     relrsp_xtitle,
258 schiefer 1.3 "Eta",etabins,"Et",etbins));
259 schiefer 1.1
260 schiefer 1.4 hEtEmfEt .push_back(hist::initHistos("Et",jetalgs[i],100,"jet E_{T} [GeV]",
261     "Emf",emfbins,"Et",etbins));
262     hEtCorrEmfEt.push_back(hist::initHistos("EtCorrEmfEt",jetalgs[i],100,
263     "jet E_{T} [GeV]",
264     "Emf",emfbins,"Et",etbins));
265     hAbsRspEmfEt.push_back(hist::initHistos("AbsRspEmfEt",jetalgs[i],100,-50,150,
266     absrsp_xtitle,
267     "Emf",emfbins,"Et",etbins));
268     hRelRspEmfEt.push_back(hist::initHistos("RelRspEmfEt",jetalgs[i],100,0,2,
269     relrsp_xtitle,
270     "Emf",emfbins,"Et",etbins));
271    
272     hEtEmfEtaEt .push_back(hist::initHistos("Et",jetalgs[i],100,"jet E_{T} [GeV]",
273     "Emf",emfbins,
274     "Eta",etabins,
275     "Et", etbins));
276     hEtCorrEmfEtaEt.push_back(hist::initHistos("EtCorrEmfEtaEt",jetalgs[i],100,
277     "jet E_{T} [GeV]",
278     "Emf",emfbins,
279     "Eta",etabins,
280     "Et", etbins));
281     hAbsRspEmfEtaEt.push_back(hist::initHistos("AbsRspEmfEtaEt",jetalgs[i],
282     100,-50,150,absrsp_xtitle,
283     "Emf",emfbins,
284     "Eta",etabins,
285     "Et", etbins));
286     hRelRspEmfEtaEt.push_back(hist::initHistos("RelRspEmfEtaEt",jetalgs[i],
287     100,0,2,relrsp_xtitle,
288     "Emf",emfbins,
289     "Eta",etabins,
290     "Et", etbins));
291    
292 schiefer 1.1
293     int nevts = el->GetN();
294     cout<<jetalgs[i]<<": "<<nevts<<" events selected."<<endl;
295 schiefer 1.4
296 schiefer 1.3 // prepare intermediate trees
297 schiefer 1.4 float genet,et,eta,phi,emf,etcorr,absrsp,relrsp;
298     TTree** tGenEt = new TTree*[netbins];
299     TTree** tEt = new TTree*[netbins];
300     TTree** tEta = new TTree*[netabins];
301     TTree** tPhi = new TTree*[nphibins];
302     TTree** tEmf = new TTree*[nemfbins];
303     TTree*** tEtEtaEt = new TTree**[netabins];
304     TTree*** tEtEmfEt = new TTree**[nemfbins];
305     TTree**** tEtEmfEtaEt = new TTree***[nemfbins];
306    
307 schiefer 1.3 for (int iet=0;iet<netbins;iet++) {
308 schiefer 1.4
309 schiefer 1.3 tGenEt[iet]=new TTree("tGenEt","tGenEt");
310     tGenEt[iet]->Branch("val", &genet, "val/F");
311     tGenEt[iet]->Branch("etcorr",&etcorr,"etcorr/F");
312 schiefer 1.4 tGenEt[iet]->Branch("absrsp",&absrsp,"absrsp/F");
313     tGenEt[iet]->Branch("relrsp",&relrsp,"relrsp/F");
314 schiefer 1.3 tGenEt[iet]->Branch("weight",&weight,"weight/F");
315 schiefer 1.4
316 schiefer 1.3 tEt[iet]=new TTree("tEt","tEt");
317     tEt[iet]->Branch("val", &et, "val/F");
318     tEt[iet]->Branch("etcorr",&etcorr,"etcorr/F");
319 schiefer 1.4 tEt[iet]->Branch("absrsp",&absrsp,"absrsp/F");
320     tEt[iet]->Branch("relrsp",&relrsp,"relrsp/F");
321 schiefer 1.3 tEt[iet]->Branch("weight",&weight,"weight/F");
322     }
323     for (int ieta=0;ieta<netabins;ieta++) {
324     tEta[ieta]=new TTree("tEta","tEta");
325     tEta[ieta]->Branch("val", &eta, "val/F");
326     tEta[ieta]->Branch("etcorr",&etcorr,"etcorr/F");
327 schiefer 1.4 tEta[ieta]->Branch("absrsp",&absrsp,"absrsp/F");
328     tEta[ieta]->Branch("relrsp",&relrsp,"relrsp/F");
329 schiefer 1.3 tEta[ieta]->Branch("weight",&weight,"weight/F");
330     }
331     for (int iphi=0;iphi<nphibins;iphi++) {
332     tPhi[iphi]=new TTree("tPhi","tPhi");
333     tPhi[iphi]->Branch("val", &phi, "val/F");
334     tPhi[iphi]->Branch("etcorr",&etcorr,"etcorr/F");
335 schiefer 1.4 tPhi[iphi]->Branch("absrsp",&absrsp,"absrsp/F");
336     tPhi[iphi]->Branch("relrsp",&relrsp,"relrsp/F");
337 schiefer 1.3 tPhi[iphi]->Branch("weight",&weight,"weight/F");
338     }
339 schiefer 1.4 for (int iemf=0;iemf<nemfbins;iemf++) {
340     tEmf[iemf]=new TTree("tEmf","tEmf");
341     tEmf[iemf]->Branch("val", &emf, "val/F");
342     tEmf[iemf]->Branch("etcorr",&etcorr,"etcorr/F");
343     tEmf[iemf]->Branch("absrsp",&absrsp,"absrsp/F");
344     tEmf[iemf]->Branch("relrsp",&relrsp,"relrsp/F");
345     tEmf[iemf]->Branch("weight",&weight,"weight/F");
346     }
347 schiefer 1.3 for (int ieta=0;ieta<netabins;ieta++) {
348     tEtEtaEt[ieta] = new TTree*[netbins];
349     for (int iet=0;iet<netbins;iet++) {
350     tEtEtaEt[ieta][iet] = new TTree("tEtEtaEt","tEtEtaEt");
351 schiefer 1.4 tEtEtaEt[ieta][iet]->Branch("val", &et, "val/F");
352 schiefer 1.3 tEtEtaEt[ieta][iet]->Branch("etcorr",&etcorr,"etcorr/F");
353 schiefer 1.4 tEtEtaEt[ieta][iet]->Branch("absrsp",&absrsp,"absrsp/F");
354     tEtEtaEt[ieta][iet]->Branch("relrsp",&relrsp,"relrsp/F");
355 schiefer 1.3 tEtEtaEt[ieta][iet]->Branch("weight",&weight,"weight/F");
356     }
357     }
358 schiefer 1.4 for (int iemf=0;iemf<nemfbins;iemf++) {
359     tEtEmfEt[iemf] = new TTree*[netbins];
360     for (int iet=0;iet<netbins;iet++) {
361     tEtEmfEt[iemf][iet] = new TTree("tEtEmfEt","tEtEmfEt");
362     tEtEmfEt[iemf][iet]->Branch("val", &et, "val/F");
363     tEtEmfEt[iemf][iet]->Branch("etcorr",&etcorr,"etcorr/F");
364     tEtEmfEt[iemf][iet]->Branch("absrsp",&absrsp,"absrsp/F");
365     tEtEmfEt[iemf][iet]->Branch("relrsp",&relrsp,"relrsp/F");
366     tEtEmfEt[iemf][iet]->Branch("weight",&weight,"weight/F");
367     }
368     }
369     for (int iemf=0;iemf<nemfbins;iemf++) {
370     tEtEmfEtaEt[iemf] = new TTree**[netabins];
371     for (int ieta=0;ieta<netabins;ieta++) {
372     tEtEmfEtaEt[iemf][ieta] = new TTree*[netbins];
373     for (int iet=0;iet<netbins;iet++) {
374     tEtEmfEtaEt[iemf][ieta][iet] = new TTree("tEtEmfEtaEt","tEtEmfEtaEt");
375     tEtEmfEtaEt[iemf][ieta][iet]->Branch("val", &et, "val/F");
376     tEtEmfEtaEt[iemf][ieta][iet]->Branch("etcorr",&etcorr,"etcorr/F");
377     tEtEmfEtaEt[iemf][ieta][iet]->Branch("absrsp",&absrsp,"absrsp/F");
378     tEtEmfEtaEt[iemf][ieta][iet]->Branch("relrsp",&relrsp,"relrsp/F");
379     tEtEmfEtaEt[iemf][ieta][iet]->Branch("weight",&weight,"weight/F");
380     }
381     }
382     }
383 schiefer 1.1
384    
385     // loop over all events
386     for (int ievt=0;ievt<nevts;ievt++) {
387    
388     t->GetEntry(el->GetEntry(ievt));
389    
390     // loop over all jets in the event
391     for (int ijt=0;ijt<njt;ijt++) {
392    
393     if (jtgendr[ijt]<0.0||jtgendr[ijt]>drmax) continue;
394    
395     if (ijt>0) continue;
396    
397 schiefer 1.3 et = jtet[ijt];
398     genet = jtgenet[ijt];
399     eta = (abseta) ? std::abs(jteta[ijt]) : jteta[ijt];
400     phi = jtphi[ijt];
401 schiefer 1.4 emf = jtemf[ijt];
402 schiefer 1.3 etcorr = (applyjes==0) ? et : et*jtjes[ijt][applyjes-1];
403 schiefer 1.1
404 schiefer 1.4 absrsp = genet - etcorr;
405     relrsp = etcorr/genet;
406 schiefer 1.2
407 schiefer 1.3 // genet
408     int igenet = hist::get_ibin(genet,etbins);
409     hAbsRspGenEt[i][igenet] ->Fill(absrsp,weight);
410     hRelRspGenEt[i][igenet] ->Fill(relrsp,weight);
411     tGenEt[igenet]->Fill();
412 schiefer 1.1
413 schiefer 1.3 // et
414     int iet = hist::get_ibin(et,etbins);
415     hAbsRspEt[i][iet]->Fill(absrsp,weight);
416     hRelRspEt[i][iet]->Fill(relrsp,weight);
417     tEt[iet]->Fill();
418 schiefer 1.1
419     // eta
420 schiefer 1.2 int ieta = hist::get_ibin(eta,etabins);
421 schiefer 1.3 hAbsRspEta[i][ieta]->Fill(absrsp,weight);
422     hRelRspEta[i][ieta]->Fill(relrsp,weight);
423     tEta[ieta]->Fill();
424    
425 schiefer 1.1 // phi
426 schiefer 1.2 int iphi = hist::get_ibin(phi,phibins);
427 schiefer 1.3 hAbsRspPhi[i][iphi]->Fill(absrsp,weight);
428     hRelRspPhi[i][iphi]->Fill(relrsp,weight);
429     tPhi[iphi]->Fill();
430    
431 schiefer 1.4 // emf
432     int iemf = hist::get_ibin(emf,emfbins);
433     hAbsRspEmf[i][iemf]->Fill(absrsp,weight);
434     hRelRspEmf[i][iemf]->Fill(relrsp,weight);
435     tEmf[iemf]->Fill();
436    
437 schiefer 1.3 // et-eta
438     hAbsRspEtaEt[i][ieta][iet]->Fill(absrsp,weight);
439     hRelRspEtaEt[i][ieta][iet]->Fill(relrsp,weight);
440     tEtEtaEt[ieta][iet]->Fill();
441 schiefer 1.1
442 schiefer 1.4 // et-emf
443     hAbsRspEmfEt[i][iemf][iet]->Fill(absrsp,weight);
444     hRelRspEmfEt[i][iemf][iet]->Fill(relrsp,weight);
445     tEtEmfEt[iemf][iet]->Fill();
446    
447     // et-eta-emf
448     hAbsRspEmfEtaEt[i][iemf][ieta][iet]->Fill(absrsp,weight);
449     hRelRspEmfEtaEt[i][iemf][ieta][iet]->Fill(relrsp,weight);
450     tEtEmfEtaEt[iemf][ieta][iet]->Fill();
451    
452 schiefer 1.1 } // jets
453     } // evts
454    
455     cout<<jetalgs[i]<<": all events processed / histograms filled."<<endl;
456 schiefer 1.3
457 schiefer 1.2
458 schiefer 1.3 // replace histograms
459 schiefer 1.4
460     replaceHistos(netbins,
461     hGenEt[i],hEtCorrGenEt[i],hAbsRspGenEt[i],hRelRspGenEt[i],tGenEt);
462     replaceHistos(netbins,hEt[i],hEtCorrEt[i],hAbsRspEt[i],hRelRspEt[i],tEt);
463     replaceHistos(netabins,hEta[i],hEtCorrEta[i],hAbsRspEta[i],hRelRspEta[i],tEta);
464     replaceHistos(nphibins,hPhi[i],hEtCorrPhi[i],hAbsRspPhi[i],hRelRspPhi[i],tPhi);
465     replaceHistos(nemfbins,hEmf[i],hEtCorrEmf[i],hAbsRspEmf[i],hRelRspEmf[i],tEmf);
466     replaceHistos(netabins,netbins,
467     hEtEtaEt[i],hEtCorrEtaEt[i],hAbsRspEtaEt[i],hRelRspEtaEt[i],
468     tEtEtaEt);
469     replaceHistos(nemfbins,netbins,
470     hEtEmfEt[i],hEtCorrEmfEt[i],hAbsRspEmfEt[i],hRelRspEmfEt[i],
471     tEtEmfEt);
472     replaceHistos(nemfbins,netabins,netbins,
473     hEtEmfEtaEt[i],hEtCorrEmfEtaEt[i],
474     hAbsRspEmfEtaEt[i],hRelRspEmfEtaEt[i],
475     tEtEmfEtaEt);
476 schiefer 1.2
477     TDirectory* d = plotfile->mkdir(jetalgs[i].c_str());
478     d->cd();
479 schiefer 1.1
480 schiefer 1.3 // fit response histos
481     fitHistos(hAbsRspGenEt[i],netbins,nsigma);
482     fitHistos(hRelRspGenEt[i],netbins,nsigma);
483    
484     fitHistos(hAbsRspEt[i], netbins,nsigma);
485     fitHistos(hRelRspEt[i], netbins,nsigma);
486    
487     fitHistos(hAbsRspEta[i], netabins,nsigma);
488     fitHistos(hRelRspEta[i], netabins,nsigma);
489    
490     fitHistos(hAbsRspPhi[i], nphibins,nsigma);
491     fitHistos(hRelRspPhi[i], nphibins,nsigma);
492    
493 schiefer 1.4 fitHistos(hAbsRspEmf[i], nemfbins,nsigma);
494     fitHistos(hRelRspEmf[i], nemfbins,nsigma);
495    
496 schiefer 1.3 fitHistos(hAbsRspEtaEt[i],netabins,netbins,nsigma);
497     fitHistos(hRelRspEtaEt[i],netabins,netbins,nsigma);
498    
499 schiefer 1.4 fitHistos(hAbsRspEmfEt[i],nemfbins,netbins,nsigma);
500     fitHistos(hRelRspEmfEt[i],nemfbins,netbins,nsigma);
501    
502     fitHistos(hAbsRspEmfEtaEt[i],nemfbins,netabins,netbins,nsigma);
503     fitHistos(hRelRspEmfEtaEt[i],nemfbins,netabins,netbins,nsigma);
504    
505 schiefer 1.3
506     // save histograms
507     saveHistos(hGenEt[i], netbins);
508     saveHistos(hEtCorrGenEt[i],netbins);
509     saveHistos(hAbsRspGenEt[i],netbins);
510     saveHistos(hRelRspGenEt[i],netbins);
511    
512     saveHistos(hEt[i], netbins);
513     saveHistos(hEtCorrEt[i], netbins);
514     saveHistos(hAbsRspEt[i], netbins);
515     saveHistos(hRelRspEt[i], netbins);
516    
517     saveHistos(hEta[i], netabins);
518     saveHistos(hEtCorrEta[i], netabins);
519     saveHistos(hAbsRspEta[i], netabins);
520     saveHistos(hRelRspEta[i], netabins);
521    
522     saveHistos(hPhi[i], nphibins);
523     saveHistos(hEtCorrPhi[i], nphibins);
524     saveHistos(hAbsRspPhi[i], nphibins);
525     saveHistos(hRelRspPhi[i], nphibins);
526    
527 schiefer 1.4 saveHistos(hEmf[i], nemfbins);
528     saveHistos(hEtCorrEmf[i], nemfbins);
529     saveHistos(hAbsRspEmf[i], nemfbins);
530     saveHistos(hRelRspEmf[i], nemfbins);
531    
532 schiefer 1.3 saveHistos(hEtEtaEt[i], netabins,netbins);
533     saveHistos(hEtCorrEtaEt[i],netabins,netbins);
534     saveHistos(hAbsRspEtaEt[i],netabins,netbins);
535     saveHistos(hRelRspEtaEt[i],netabins,netbins);
536 schiefer 1.4
537     saveHistos(hEtEmfEt[i], nemfbins,netbins);
538     saveHistos(hEtCorrEmfEt[i],nemfbins,netbins);
539     saveHistos(hAbsRspEmfEt[i],nemfbins,netbins);
540     saveHistos(hRelRspEmfEt[i],nemfbins,netbins);
541    
542     saveHistos(hEtEmfEtaEt[i], nemfbins,netabins,netbins);
543     saveHistos(hEtCorrEmfEtaEt[i],nemfbins,netabins,netbins);
544     saveHistos(hAbsRspEmfEtaEt[i],nemfbins,netabins,netbins);
545     saveHistos(hRelRspEmfEtaEt[i],nemfbins,netabins,netbins);
546 schiefer 1.1 }
547    
548 schiefer 1.2 plotfile->Write();
549     plotfile->Close();
550     delete plotfile;
551 schiefer 1.1
552     return 0;
553     }
554    
555    
556    
557     ////////////////////////////////////////////////////////////////////////////////
558     // implementation of global functions
559     ////////////////////////////////////////////////////////////////////////////////
560    
561     //______________________________________________________________________________
562 schiefer 1.4 void replaceHistos(int nbins,
563     TH1F**& hVar,TH1F**& hEtCorr,TH1F**& hAbsRsp,TH1F**& hRelRsp,
564     TTree**& tVar)
565 schiefer 1.1 {
566 schiefer 1.3 TH1::SetDefaultSumw2();
567 schiefer 1.1 for (int ibin=0;ibin<nbins;ibin++) {
568 schiefer 1.4
569 schiefer 1.3 TH1F* h(0);
570 schiefer 1.4
571     if (tVar[ibin]->GetEntries()==0) continue;
572 schiefer 1.1
573 schiefer 1.3 tVar[ibin]->Project("h(50)","val","weight*(1)");
574     h = (TH1F*)gROOT->FindObject("h");
575     h->SetNameTitle(hVar[ibin]->GetName(),hVar[ibin]->GetTitle());
576     h->SetXTitle(hVar[ibin]->GetXaxis()->GetTitle());
577     h->SetYTitle(hVar[ibin]->GetYaxis()->GetTitle());
578     delete hVar[ibin];
579     hVar[ibin] = h;
580    
581     tVar[ibin]->Project("h(50)","etcorr","weight*(1)");
582     h = (TH1F*)gROOT->FindObject("h");
583     h->SetNameTitle(hEtCorr[ibin]->GetName(),hEtCorr[ibin]->GetTitle());
584     h->SetXTitle(hEtCorr[ibin]->GetXaxis()->GetTitle());
585     h->SetYTitle(hEtCorr[ibin]->GetYaxis()->GetTitle());
586     delete hEtCorr[ibin];
587     hEtCorr[ibin] = h;
588 schiefer 1.4
589     tVar[ibin]->Project("h(50)","absrsp","weight*(1)");
590     h = (TH1F*)gROOT->FindObject("h");
591     h->SetNameTitle(hAbsRsp[ibin]->GetName(),hAbsRsp[ibin]->GetTitle());
592     h->SetXTitle(hAbsRsp[ibin]->GetXaxis()->GetTitle());
593     h->SetYTitle(hAbsRsp[ibin]->GetYaxis()->GetTitle());
594     delete hAbsRsp[ibin];
595     hAbsRsp[ibin] = h;
596    
597     tVar[ibin]->Project("h(50)","relrsp","weight*(1)");
598     h = (TH1F*)gROOT->FindObject("h");
599     h->SetNameTitle(hRelRsp[ibin]->GetName(),hRelRsp[ibin]->GetTitle());
600     h->SetXTitle(hRelRsp[ibin]->GetXaxis()->GetTitle());
601     h->SetYTitle(hRelRsp[ibin]->GetYaxis()->GetTitle());
602     delete hRelRsp[ibin];
603     hRelRsp[ibin] = h;
604    
605 schiefer 1.3 delete tVar[ibin];
606 schiefer 1.1 }
607    
608 schiefer 1.3 delete [] tVar;
609 schiefer 1.1 }
610    
611 schiefer 1.4
612     //______________________________________________________________________________
613     void replaceHistos(int nbins1,int nbins2,
614     TH1F***& hVar,TH1F***& hEtCorr,TH1F***& hAbsRsp,TH1F***& hRelRsp,
615     TTree***& tVar)
616     {
617     for (int i1=0;i1<nbins1;i1++)
618     replaceHistos(nbins2,hVar[i1],hEtCorr[i1],hAbsRsp[i1],hRelRsp[i1],tVar[i1]);
619     }
620    
621    
622     //______________________________________________________________________________
623     void replaceHistos(int nbins1,int nbins2,int nbins3,
624     TH1F****& hVar,TH1F****& hEtCorr,TH1F****&hAbsRsp,TH1F****&hRelRsp,
625     TTree****& tVar)
626     {
627     for (int i1=0;i1<nbins1;i1++)
628     replaceHistos(nbins2,nbins3,
629     hVar[i1],hEtCorr[i1],hAbsRsp[i1],hRelRsp[i1],tVar[i1]);
630     }
631    
632    
633 schiefer 1.1 //______________________________________________________________________________
634 schiefer 1.3 void fitHistos(TH1F** histos,int n,double nsigma)
635 schiefer 1.1 {
636 schiefer 1.3 for (int i=0;i<n;i++) {
637     float nrm = histos[i]->Integral();
638     float mean = histos[i]->GetMean();
639     float rms = histos[i]->GetRMS();
640     TF1* f = new TF1("gaus","gaus",mean-nsigma*rms,mean+nsigma*rms);
641     f->SetParameter(0,nrm);
642     f->SetParameter(1,mean);
643     f->SetParameter(2,rms);
644     f->SetLineWidth(1);
645     f->SetLineColor(kRed);
646     histos[i]->Fit(f,"QR");
647 schiefer 1.1 }
648     }
649    
650    
651     //______________________________________________________________________________
652 schiefer 1.3 void fitHistos(TH1F*** histos,int n1,int n2,double nsigma)
653 schiefer 1.1 {
654 schiefer 1.3 for (int i1=0;i1<n1;i1++) fitHistos(histos[i1],n2,nsigma);
655 schiefer 1.1 }
656    
657    
658     //______________________________________________________________________________
659 schiefer 1.4 void fitHistos(TH1F**** histos,int n1,int n2,int n3,double nsigma)
660     {
661     for (int i1=0;i1<n1;i1++) fitHistos(histos[i1],n2,n3,nsigma);
662     }
663    
664    
665     //______________________________________________________________________________
666 schiefer 1.3 void saveHistos(TH1F** histos,int n)
667 schiefer 1.1 {
668 schiefer 1.3 for (int i=0;i<n;i++) histos[i]->Write();
669 schiefer 1.1 }
670    
671    
672     //______________________________________________________________________________
673 schiefer 1.3 void saveHistos(TH1F*** histos,int n1,int n2)
674 schiefer 1.1 {
675 schiefer 1.3 for (int i1=0;i1<n1;i1++) saveHistos(histos[i1],n2);
676 schiefer 1.1 }
677 schiefer 1.4
678    
679     //______________________________________________________________________________
680     void saveHistos(TH1F**** histos,int n1,int n2,int n3)
681     {
682     for (int i1=0;i1<n1;i1++) saveHistos(histos[i1],n2,n3);
683     }