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