ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/SchieferD/jetcalib/mcresponse_x.cpp
Revision: 1.8
Committed: Sun Sep 2 14:53:13 2007 UTC (17 years, 7 months ago) by schiefer
Branch: MAIN
Changes since 1.7: +37 -31 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.8
377     if (netbins>0) {
378     for (int ieta=0;ieta<netabins;ieta++) {
379     tEtEtaEt[ieta] = new TTree*[netbins];
380     for (int iet=0;iet<netbins;iet++) {
381     tEtEtaEt[ieta][iet] = new TTree("tEtEtaEt","tEtEtaEt");
382     tEtEtaEt[ieta][iet]->Branch("val", &et, "val/F");
383     tEtEtaEt[ieta][iet]->Branch("etcorr",&etcorr,"etcorr/F");
384     tEtEtaEt[ieta][iet]->Branch("absrsp",&absrsp,"absrsp/F");
385     tEtEtaEt[ieta][iet]->Branch("relrsp",&relrsp,"relrsp/F");
386     tEtEtaEt[ieta][iet]->Branch("weight",&weight,"weight/F");
387     }
388 schiefer 1.3 }
389 schiefer 1.8 for (int iemf=0;iemf<nemfbins;iemf++) {
390     tEtEmfEt[iemf] = new TTree*[netbins];
391     for (int iet=0;iet<netbins;iet++) {
392     tEtEmfEt[iemf][iet] = new TTree("tEtEmfEt","tEtEmfEt");
393     tEtEmfEt[iemf][iet]->Branch("val", &et, "val/F");
394     tEtEmfEt[iemf][iet]->Branch("etcorr",&etcorr,"etcorr/F");
395     tEtEmfEt[iemf][iet]->Branch("absrsp",&absrsp,"absrsp/F");
396     tEtEmfEt[iemf][iet]->Branch("relrsp",&relrsp,"relrsp/F");
397     tEtEmfEt[iemf][iet]->Branch("weight",&weight,"weight/F");
398     }
399 schiefer 1.4 }
400 schiefer 1.8 if (netabins>0) {
401     for (int iemf=0;iemf<nemfbins;iemf++) {
402     tEtEmfEtaEt[iemf] = new TTree**[netabins];
403     for (int ieta=0;ieta<netabins;ieta++) {
404     tEtEmfEtaEt[iemf][ieta] = new TTree*[netbins];
405     for (int iet=0;iet<netbins;iet++) {
406     tEtEmfEtaEt[iemf][ieta][iet] = new TTree("tEtEmfEtaEt",
407     "tEtEmfEtaEt");
408     tEtEmfEtaEt[iemf][ieta][iet]->Branch("val", &et, "val/F");
409     tEtEmfEtaEt[iemf][ieta][iet]->Branch("etcorr",&etcorr,"etcorr/F");
410     tEtEmfEtaEt[iemf][ieta][iet]->Branch("absrsp",&absrsp,"absrsp/F");
411     tEtEmfEtaEt[iemf][ieta][iet]->Branch("relrsp",&relrsp,"relrsp/F");
412     tEtEmfEtaEt[iemf][ieta][iet]->Branch("weight",&weight,"weight/F");
413     }
414     }
415 schiefer 1.4 }
416     }
417     }
418 schiefer 1.1
419 schiefer 1.5 cout<<jetalgs[i]<<": trees created."<<endl;
420    
421 schiefer 1.1
422     // loop over all events
423     for (int ievt=0;ievt<nevts;ievt++) {
424    
425     t->GetEntry(el->GetEntry(ievt));
426    
427     // loop over all jets in the event
428     for (int ijt=0;ijt<njt;ijt++) {
429    
430     if (jtgendr[ijt]<0.0||jtgendr[ijt]>drmax) continue;
431    
432 schiefer 1.7 //if (ijt>ijtmax) continue;
433    
434 schiefer 1.3 et = jtet[ijt];
435     genet = jtgenet[ijt];
436     eta = (abseta) ? std::abs(jteta[ijt]) : jteta[ijt];
437     phi = jtphi[ijt];
438 schiefer 1.4 emf = jtemf[ijt];
439 schiefer 1.6
440 schiefer 1.3 etcorr = (applyjes==0) ? et : et*jtjes[ijt][applyjes-1];
441 schiefer 1.1
442 schiefer 1.4 absrsp = genet - etcorr;
443     relrsp = etcorr/genet;
444 schiefer 1.2
445 schiefer 1.6 if (TMath::IsNaN(emf)) emf = 0.0;
446    
447 schiefer 1.7 int igenet = (netbins>0) ? hist::get_ibin(genet,etbins) : -1;
448     int iet = (netbins>0) ? hist::get_ibin(et,etbins) : -1;
449     int ieta = (netabins>0) ? hist::get_ibin(eta,etabins) : -1;
450     int iphi = (nphibins>0) ? hist::get_ibin(phi,phibins) : -1;
451     int iemf = (nemfbins>0) ? hist::get_ibin(emf,emfbins) : -1;
452    
453     if (netbins>0) tGenEt[igenet]->Fill();
454     if (netbins>0) tEt[iet]->Fill();
455     if (netabins>0) tEta[ieta]->Fill();
456     if (nphibins>0) tPhi[iphi]->Fill();
457     if (nemfbins>0) tEmf[iemf]->Fill();
458    
459     if (netbins>0&&netabins>0) tEtEtaEt[ieta][iet]->Fill();
460     if (netbins>0&&nemfbins>0) tEtEmfEt[iemf][iet]->Fill();
461    
462     if (netbins>0&&netabins>0&&nemfbins>0)
463     tEtEmfEtaEt[iemf][ieta][iet]->Fill();
464 schiefer 1.4
465 schiefer 1.1 } // jets
466     } // evts
467    
468     cout<<jetalgs[i]<<": all events processed / histograms filled."<<endl;
469 schiefer 1.3
470 schiefer 1.2
471 schiefer 1.3 // replace histograms
472 schiefer 1.4
473     replaceHistos(netbins,
474 schiefer 1.5 hGenEt,hEtCorrGenEt,hAbsRspGenEt,hRelRspGenEt,tGenEt);
475     replaceHistos(netbins,hEt,hEtCorrEt,hAbsRspEt,hRelRspEt,tEt);
476     replaceHistos(netabins,hEta,hEtCorrEta,hAbsRspEta,hRelRspEta,tEta);
477     replaceHistos(nphibins,hPhi,hEtCorrPhi,hAbsRspPhi,hRelRspPhi,tPhi);
478     replaceHistos(nemfbins,hEmf,hEtCorrEmf,hAbsRspEmf,hRelRspEmf,tEmf);
479 schiefer 1.4 replaceHistos(netabins,netbins,
480 schiefer 1.5 hEtEtaEt,hEtCorrEtaEt,hAbsRspEtaEt,hRelRspEtaEt,
481 schiefer 1.4 tEtEtaEt);
482     replaceHistos(nemfbins,netbins,
483 schiefer 1.5 hEtEmfEt,hEtCorrEmfEt,hAbsRspEmfEt,hRelRspEmfEt,
484 schiefer 1.4 tEtEmfEt);
485     replaceHistos(nemfbins,netabins,netbins,
486 schiefer 1.5 hEtEmfEtaEt,hEtCorrEmfEtaEt,
487     hAbsRspEmfEtaEt,hRelRspEmfEtaEt,
488 schiefer 1.4 tEtEmfEtaEt);
489 schiefer 1.2
490 schiefer 1.5 cout<<jetalgs[i]<<": all histograms replaced, and tree deleted."<<endl;
491    
492 schiefer 1.6
493     TDirectory* d = plotfile->GetDirectory(jetalgs[i].c_str());
494     if (0!=d) plotfile->rmdir(jetalgs[i].c_str());
495     d = plotfile->mkdir(jetalgs[i].c_str());
496 schiefer 1.2 d->cd();
497 schiefer 1.1
498 schiefer 1.3 // fit response histos
499 schiefer 1.5 fitHistos(hAbsRspGenEt,netbins,nsigma);
500     fitHistos(hRelRspGenEt,netbins,nsigma);
501 schiefer 1.3
502 schiefer 1.5 fitHistos(hAbsRspEt, netbins,nsigma);
503     fitHistos(hRelRspEt, netbins,nsigma);
504 schiefer 1.3
505 schiefer 1.5 fitHistos(hAbsRspEta, netabins,nsigma);
506     fitHistos(hRelRspEta, netabins,nsigma);
507 schiefer 1.3
508 schiefer 1.5 fitHistos(hAbsRspPhi, nphibins,nsigma);
509     fitHistos(hRelRspPhi, nphibins,nsigma);
510 schiefer 1.3
511 schiefer 1.5 fitHistos(hAbsRspEmf, nemfbins,nsigma);
512     fitHistos(hRelRspEmf, nemfbins,nsigma);
513 schiefer 1.4
514 schiefer 1.5 fitHistos(hAbsRspEtaEt,netabins,netbins,nsigma);
515     fitHistos(hRelRspEtaEt,netabins,netbins,nsigma);
516 schiefer 1.3
517 schiefer 1.5 fitHistos(hAbsRspEmfEt,nemfbins,netbins,nsigma);
518     fitHistos(hRelRspEmfEt,nemfbins,netbins,nsigma);
519 schiefer 1.4
520 schiefer 1.5 fitHistos(hAbsRspEmfEtaEt,nemfbins,netabins,netbins,nsigma);
521     fitHistos(hRelRspEmfEtaEt,nemfbins,netabins,netbins,nsigma);
522    
523     cout<<jetalgs[i]<<": fits done."<<endl;
524 schiefer 1.4
525 schiefer 1.3
526     // save histograms
527 schiefer 1.5 saveHistos(hGenEt, netbins);
528     saveHistos(hEtCorrGenEt,netbins);
529     saveHistos(hAbsRspGenEt,netbins);
530     saveHistos(hRelRspGenEt,netbins);
531    
532     saveHistos(hEt, netbins);
533     saveHistos(hEtCorrEt, netbins);
534     saveHistos(hAbsRspEt, netbins);
535     saveHistos(hRelRspEt, netbins);
536    
537     saveHistos(hEta, netabins);
538     saveHistos(hEtCorrEta, netabins);
539     saveHistos(hAbsRspEta, netabins);
540     saveHistos(hRelRspEta, netabins);
541    
542     saveHistos(hPhi, nphibins);
543     saveHistos(hEtCorrPhi, nphibins);
544     saveHistos(hAbsRspPhi, nphibins);
545     saveHistos(hRelRspPhi, nphibins);
546    
547     saveHistos(hEmf, nemfbins);
548     saveHistos(hEtCorrEmf, nemfbins);
549     saveHistos(hAbsRspEmf, nemfbins);
550     saveHistos(hRelRspEmf, nemfbins);
551    
552     saveHistos(hEtEtaEt, netabins,netbins);
553     saveHistos(hEtCorrEtaEt,netabins,netbins);
554     saveHistos(hAbsRspEtaEt,netabins,netbins);
555     saveHistos(hRelRspEtaEt,netabins,netbins);
556    
557     saveHistos(hEtEmfEt, nemfbins,netbins);
558     saveHistos(hEtCorrEmfEt,nemfbins,netbins);
559     saveHistos(hAbsRspEmfEt,nemfbins,netbins);
560     saveHistos(hRelRspEmfEt,nemfbins,netbins);
561    
562     saveHistos(hEtEmfEtaEt, nemfbins,netabins,netbins);
563     saveHistos(hEtCorrEmfEtaEt,nemfbins,netabins,netbins);
564     saveHistos(hAbsRspEmfEtaEt,nemfbins,netabins,netbins);
565     saveHistos(hRelRspEmfEtaEt,nemfbins,netabins,netbins);
566    
567     plotfile->Write();
568    
569     cout<<jetalgs[i]<<": histograms saved."<<endl;
570 schiefer 1.1 }
571    
572 schiefer 1.2 plotfile->Write();
573     plotfile->Close();
574     delete plotfile;
575 schiefer 1.1
576     return 0;
577     }
578    
579    
580    
581     ////////////////////////////////////////////////////////////////////////////////
582     // implementation of global functions
583     ////////////////////////////////////////////////////////////////////////////////
584    
585     //______________________________________________________________________________
586 schiefer 1.4 void replaceHistos(int nbins,
587     TH1F**& hVar,TH1F**& hEtCorr,TH1F**& hAbsRsp,TH1F**& hRelRsp,
588     TTree**& tVar)
589 schiefer 1.1 {
590 schiefer 1.7 if (nbins==0) return;
591    
592 schiefer 1.3 TH1::SetDefaultSumw2();
593 schiefer 1.1 for (int ibin=0;ibin<nbins;ibin++) {
594 schiefer 1.4
595 schiefer 1.3 TH1F* h(0);
596 schiefer 1.4
597     if (tVar[ibin]->GetEntries()==0) continue;
598 schiefer 1.1
599 schiefer 1.3 tVar[ibin]->Project("h(50)","val","weight*(1)");
600     h = (TH1F*)gROOT->FindObject("h");
601     h->SetNameTitle(hVar[ibin]->GetName(),hVar[ibin]->GetTitle());
602     h->SetXTitle(hVar[ibin]->GetXaxis()->GetTitle());
603     h->SetYTitle(hVar[ibin]->GetYaxis()->GetTitle());
604     delete hVar[ibin];
605     hVar[ibin] = h;
606    
607     tVar[ibin]->Project("h(50)","etcorr","weight*(1)");
608     h = (TH1F*)gROOT->FindObject("h");
609     h->SetNameTitle(hEtCorr[ibin]->GetName(),hEtCorr[ibin]->GetTitle());
610     h->SetXTitle(hEtCorr[ibin]->GetXaxis()->GetTitle());
611     h->SetYTitle(hEtCorr[ibin]->GetYaxis()->GetTitle());
612     delete hEtCorr[ibin];
613     hEtCorr[ibin] = h;
614 schiefer 1.4
615     tVar[ibin]->Project("h(50)","absrsp","weight*(1)");
616     h = (TH1F*)gROOT->FindObject("h");
617     h->SetNameTitle(hAbsRsp[ibin]->GetName(),hAbsRsp[ibin]->GetTitle());
618     h->SetXTitle(hAbsRsp[ibin]->GetXaxis()->GetTitle());
619     h->SetYTitle(hAbsRsp[ibin]->GetYaxis()->GetTitle());
620     delete hAbsRsp[ibin];
621     hAbsRsp[ibin] = h;
622    
623     tVar[ibin]->Project("h(50)","relrsp","weight*(1)");
624     h = (TH1F*)gROOT->FindObject("h");
625     h->SetNameTitle(hRelRsp[ibin]->GetName(),hRelRsp[ibin]->GetTitle());
626     h->SetXTitle(hRelRsp[ibin]->GetXaxis()->GetTitle());
627     h->SetYTitle(hRelRsp[ibin]->GetYaxis()->GetTitle());
628     delete hRelRsp[ibin];
629     hRelRsp[ibin] = h;
630    
631 schiefer 1.3 delete tVar[ibin];
632 schiefer 1.1 }
633    
634 schiefer 1.3 delete [] tVar;
635 schiefer 1.1 }
636    
637 schiefer 1.4
638     //______________________________________________________________________________
639     void replaceHistos(int nbins1,int nbins2,
640     TH1F***& hVar,TH1F***& hEtCorr,TH1F***& hAbsRsp,TH1F***& hRelRsp,
641     TTree***& tVar)
642     {
643 schiefer 1.7 if (nbins1==0) return;
644    
645 schiefer 1.4 for (int i1=0;i1<nbins1;i1++)
646     replaceHistos(nbins2,hVar[i1],hEtCorr[i1],hAbsRsp[i1],hRelRsp[i1],tVar[i1]);
647 schiefer 1.5 delete [] tVar;
648 schiefer 1.4 }
649    
650    
651     //______________________________________________________________________________
652     void replaceHistos(int nbins1,int nbins2,int nbins3,
653     TH1F****& hVar,TH1F****& hEtCorr,TH1F****&hAbsRsp,TH1F****&hRelRsp,
654     TTree****& tVar)
655     {
656 schiefer 1.7 if (nbins1==0) return;
657 schiefer 1.4 for (int i1=0;i1<nbins1;i1++)
658     replaceHistos(nbins2,nbins3,
659     hVar[i1],hEtCorr[i1],hAbsRsp[i1],hRelRsp[i1],tVar[i1]);
660 schiefer 1.5 delete [] tVar;
661 schiefer 1.4 }
662    
663    
664 schiefer 1.1 //______________________________________________________________________________
665 schiefer 1.3 void fitHistos(TH1F** histos,int n,double nsigma)
666 schiefer 1.1 {
667 schiefer 1.3 for (int i=0;i<n;i++) {
668     float nrm = histos[i]->Integral();
669     float mean = histos[i]->GetMean();
670     float rms = histos[i]->GetRMS();
671     TF1* f = new TF1("gaus","gaus",mean-nsigma*rms,mean+nsigma*rms);
672     f->SetParameter(0,nrm);
673     f->SetParameter(1,mean);
674     f->SetParameter(2,rms);
675     f->SetLineWidth(1);
676     f->SetLineColor(kRed);
677     histos[i]->Fit(f,"QR");
678 schiefer 1.1 }
679     }
680    
681    
682     //______________________________________________________________________________
683 schiefer 1.3 void fitHistos(TH1F*** histos,int n1,int n2,double nsigma)
684 schiefer 1.1 {
685 schiefer 1.3 for (int i1=0;i1<n1;i1++) fitHistos(histos[i1],n2,nsigma);
686 schiefer 1.1 }
687    
688    
689     //______________________________________________________________________________
690 schiefer 1.4 void fitHistos(TH1F**** histos,int n1,int n2,int n3,double nsigma)
691     {
692     for (int i1=0;i1<n1;i1++) fitHistos(histos[i1],n2,n3,nsigma);
693     }
694    
695    
696     //______________________________________________________________________________
697 schiefer 1.5 void saveHistos(TH1F**& histos,int n)
698 schiefer 1.1 {
699 schiefer 1.7 if (n==0) return;
700 schiefer 1.5 for (int i=0;i<n;i++) {
701     histos[i]->Write();
702     delete histos[i];
703     }
704     delete [] histos;
705 schiefer 1.1 }
706    
707    
708     //______________________________________________________________________________
709 schiefer 1.5 void saveHistos(TH1F***& histos,int n1,int n2)
710 schiefer 1.1 {
711 schiefer 1.7 if (n1==0) return;
712 schiefer 1.3 for (int i1=0;i1<n1;i1++) saveHistos(histos[i1],n2);
713 schiefer 1.5 delete [] histos;
714 schiefer 1.1 }
715 schiefer 1.4
716    
717     //______________________________________________________________________________
718 schiefer 1.5 void saveHistos(TH1F****& histos,int n1,int n2,int n3)
719 schiefer 1.4 {
720 schiefer 1.7 if (n1==0) return;
721 schiefer 1.4 for (int i1=0;i1<n1;i1++) saveHistos(histos[i1],n2,n3);
722 schiefer 1.5 delete [] histos;
723 schiefer 1.4 }