ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/SchieferD/jetcalib/mcresponse_x.cpp
Revision: 1.5
Committed: Fri Aug 31 12:03:49 2007 UTC (17 years, 8 months ago) by schiefer
Branch: MAIN
Changes since 1.4: +173 -213 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 schiefer 1.5 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.3 argc=3; argv[1]="-l"; argv[2]="-b";
127 schiefer 1.4 TRint* app = new TRint(argv[0],&argc,argv); app->Argc();
128 schiefer 1.1
129    
130     //
131     // loop over all input files
132     //
133 schiefer 1.3 TFile* plotfile = new TFile(output.c_str(),"RECREATE");
134 schiefer 1.2
135 schiefer 1.1 for (unsigned int i=0;i<input.size();++i) {
136 schiefer 1.3 TFile* f = new TFile(input[i].c_str(),"UPDATE"); if (!f->IsOpen()) return 0;
137 schiefer 1.1 TTree* t = (TTree*)f->Get(treename.c_str());
138    
139     TEventList* el = new TEventList("sel","sel");
140     t->Draw(">>sel",selection.c_str());
141    
142     t->SetBranchAddress("weight", &weight);
143     t->SetBranchAddress("njt", &njt);
144 schiefer 1.3 t->SetBranchAddress("jtet", jtet);
145 schiefer 1.1 t->SetBranchAddress("jtpt", jtpt);
146     t->SetBranchAddress("jteta", jteta);
147     t->SetBranchAddress("jtphi", jtphi);
148 schiefer 1.4 t->SetBranchAddress("jtemf", jtemf);
149 schiefer 1.3 t->SetBranchAddress("jtgenet",jtgenet);
150 schiefer 1.1 t->SetBranchAddress("jtgenpt",jtgenpt);
151     t->SetBranchAddress("jtgendr",jtgendr);
152    
153     if (t->FindBranch("jtjes")) t->SetBranchAddress("jtjes", jtjes);
154     else
155     for (int i1=0;i1<100;i1++)
156     for (int i2=0;i2<3;i2++) jtjes[i1][i2]=1.0;
157    
158 schiefer 1.3 string absrsp_xtitle = "E_{T}^{gen}-E_{T} [GeV]";
159     string relrsp_xtitle = "E_{T}^{rec}/E_{T}";
160 schiefer 1.5
161     TH1F** hGenEt = hist::initHistos("GenEt",jetalgs[i],100,
162     "jet E_{T}^{gen} [GeV]","GenEt",etbins);
163     TH1F** hEtCorrGenEt = hist::initHistos("EtCorrGenEt",jetalgs[i],100,
164     "jet E_{T}^{gen} [GeV]","GenEt",etbins);
165     TH1F** hAbsRspGenEt = hist::initHistos("AbsRspGenEt",jetalgs[i],100,-50,150,
166     absrsp_xtitle,"GenEt",etbins);
167     TH1F** hRelRspGenEt = hist::initHistos("RelRspGenEt",jetalgs[i],100,0,2,
168     relrsp_xtitle,"GenEt",etbins);
169    
170     TH1F** hEt = hist::initHistos("Et",jetalgs[i],100,
171     "jet E_{T} [GeV]","Et",etbins);
172     TH1F** hEtCorrEt = hist::initHistos("EtCorrEt",jetalgs[i],100,
173     "jet E_{T} [GeV]","Et",etbins);
174     TH1F** hAbsRspEt = hist::initHistos("AbsRspEt",jetalgs[i],100,-50,150,
175     absrsp_xtitle,"Et",etbins);
176     TH1F** hRelRspEt = hist::initHistos("RelRspEt",jetalgs[i],100,0,2,
177     relrsp_xtitle,"Et",etbins);
178    
179     TH1F** hEta = hist::initHistos("Eta",jetalgs[i],100,
180     "jet #eta","Eta",etabins);
181     TH1F** hEtCorrEta = hist::initHistos("EtCorrEta",jetalgs[i],100,
182     "jet E_{T} [GeV]","Eta",etabins);
183     TH1F** hAbsRspEta = hist::initHistos("AbsRspEta",jetalgs[i],100,-50,150,
184     absrsp_xtitle,"Eta",etabins);
185     TH1F** hRelRspEta = hist::initHistos("RelRspEta",jetalgs[i],100,0,2,
186     relrsp_xtitle,"Eta",etabins);
187    
188     TH1F** hPhi = hist::initHistos("Phi",jetalgs[i],100,
189     "jet #phi","Phi",phibins);
190     TH1F** hEtCorrPhi = hist::initHistos("EtCorrPhi",jetalgs[i],100,
191     "jet E_{T} [Gev]","Phi",phibins);
192     TH1F** hAbsRspPhi = hist::initHistos("AbsRspPhi",jetalgs[i],100,-50,150,
193     absrsp_xtitle,"Phi",phibins);
194     TH1F** hRelRspPhi = hist::initHistos("RelRspPhi",jetalgs[i],100,0,2,
195     relrsp_xtitle,"Phi",phibins);
196    
197     TH1F** hEmf = hist::initHistos("Emf",jetalgs[i],100,
198     "jet emf","Emf",emfbins);
199     TH1F** hEtCorrEmf = hist::initHistos("EtCorrEmf",jetalgs[i],100,
200     "jet E_{T} [Gev]","Emf",emfbins);
201     TH1F** hAbsRspEmf = hist::initHistos("AbsRspEmf",jetalgs[i],100,-50,150,
202     absrsp_xtitle,"Emf",emfbins);
203     TH1F** hRelRspEmf = hist::initHistos("RelRspEmf",jetalgs[i],100,0,2,
204     relrsp_xtitle,"Emf",emfbins);
205    
206     TH1F*** hEtEtaEt = hist::initHistos("Et",jetalgs[i],100,"jet E_{T} [GeV]",
207     "Eta",etabins,"Et",etbins);
208     TH1F*** hEtCorrEtaEt = hist::initHistos("EtCorrEtaEt",jetalgs[i],100,
209 schiefer 1.4 "jet E_{T} [GeV]",
210 schiefer 1.5 "Eta",etabins,"Et",etbins);
211     TH1F*** hAbsRspEtaEt = hist::initHistos("AbsRspEtaEt",jetalgs[i],100,-50,150,
212 schiefer 1.4 absrsp_xtitle,
213 schiefer 1.5 "Eta",etabins,"Et",etbins);
214     TH1F*** hRelRspEtaEt = hist::initHistos("RelRspEtaEt",jetalgs[i],100,0,2,
215 schiefer 1.4 relrsp_xtitle,
216 schiefer 1.5 "Eta",etabins,"Et",etbins);
217 schiefer 1.1
218 schiefer 1.5 TH1F*** hEtEmfEt = hist::initHistos("Et",jetalgs[i],100,"jet E_{T} [GeV]",
219     "Emf",emfbins,"Et",etbins);
220     TH1F*** hEtCorrEmfEt = hist::initHistos("EtCorrEmfEt",jetalgs[i],100,
221 schiefer 1.4 "jet E_{T} [GeV]",
222 schiefer 1.5 "Emf",emfbins,"Et",etbins);
223     TH1F*** hAbsRspEmfEt = hist::initHistos("AbsRspEmfEt",jetalgs[i],100,-50,150,
224 schiefer 1.4 absrsp_xtitle,
225 schiefer 1.5 "Emf",emfbins,"Et",etbins);
226     TH1F*** hRelRspEmfEt = hist::initHistos("RelRspEmfEt",jetalgs[i],100,0,2,
227 schiefer 1.4 relrsp_xtitle,
228 schiefer 1.5 "Emf",emfbins,"Et",etbins);
229 schiefer 1.4
230 schiefer 1.5 TH1F**** hEtEmfEtaEt = hist::initHistos("Et",jetalgs[i],100,"jet E_{T} [GeV]",
231     "Emf",emfbins,
232     "Eta",etabins,
233     "Et", etbins);
234     TH1F**** hEtCorrEmfEtaEt = hist::initHistos("EtCorrEmfEtaEt",jetalgs[i],100,
235     "jet E_{T} [GeV]",
236     "Emf",emfbins,
237     "Eta",etabins,
238     "Et", etbins);
239     TH1F**** hAbsRspEmfEtaEt = hist::initHistos("AbsRspEmfEtaEt",jetalgs[i],
240     100,-50,150,absrsp_xtitle,
241     "Emf",emfbins,
242     "Eta",etabins,
243     "Et", etbins);
244     TH1F**** hRelRspEmfEtaEt = hist::initHistos("RelRspEmfEtaEt",jetalgs[i],
245     100,0,2,relrsp_xtitle,
246     "Emf",emfbins,
247     "Eta",etabins,
248     "Et", etbins);
249 schiefer 1.4
250 schiefer 1.1
251     int nevts = el->GetN();
252     cout<<jetalgs[i]<<": "<<nevts<<" events selected."<<endl;
253 schiefer 1.4
254 schiefer 1.3 // prepare intermediate trees
255 schiefer 1.4 float genet,et,eta,phi,emf,etcorr,absrsp,relrsp;
256     TTree** tGenEt = new TTree*[netbins];
257     TTree** tEt = new TTree*[netbins];
258     TTree** tEta = new TTree*[netabins];
259     TTree** tPhi = new TTree*[nphibins];
260     TTree** tEmf = new TTree*[nemfbins];
261     TTree*** tEtEtaEt = new TTree**[netabins];
262     TTree*** tEtEmfEt = new TTree**[nemfbins];
263     TTree**** tEtEmfEtaEt = new TTree***[nemfbins];
264    
265 schiefer 1.3 for (int iet=0;iet<netbins;iet++) {
266 schiefer 1.4
267 schiefer 1.3 tGenEt[iet]=new TTree("tGenEt","tGenEt");
268     tGenEt[iet]->Branch("val", &genet, "val/F");
269     tGenEt[iet]->Branch("etcorr",&etcorr,"etcorr/F");
270 schiefer 1.4 tGenEt[iet]->Branch("absrsp",&absrsp,"absrsp/F");
271     tGenEt[iet]->Branch("relrsp",&relrsp,"relrsp/F");
272 schiefer 1.3 tGenEt[iet]->Branch("weight",&weight,"weight/F");
273 schiefer 1.4
274 schiefer 1.3 tEt[iet]=new TTree("tEt","tEt");
275     tEt[iet]->Branch("val", &et, "val/F");
276     tEt[iet]->Branch("etcorr",&etcorr,"etcorr/F");
277 schiefer 1.4 tEt[iet]->Branch("absrsp",&absrsp,"absrsp/F");
278     tEt[iet]->Branch("relrsp",&relrsp,"relrsp/F");
279 schiefer 1.3 tEt[iet]->Branch("weight",&weight,"weight/F");
280     }
281     for (int ieta=0;ieta<netabins;ieta++) {
282     tEta[ieta]=new TTree("tEta","tEta");
283     tEta[ieta]->Branch("val", &eta, "val/F");
284     tEta[ieta]->Branch("etcorr",&etcorr,"etcorr/F");
285 schiefer 1.4 tEta[ieta]->Branch("absrsp",&absrsp,"absrsp/F");
286     tEta[ieta]->Branch("relrsp",&relrsp,"relrsp/F");
287 schiefer 1.3 tEta[ieta]->Branch("weight",&weight,"weight/F");
288     }
289     for (int iphi=0;iphi<nphibins;iphi++) {
290     tPhi[iphi]=new TTree("tPhi","tPhi");
291     tPhi[iphi]->Branch("val", &phi, "val/F");
292     tPhi[iphi]->Branch("etcorr",&etcorr,"etcorr/F");
293 schiefer 1.4 tPhi[iphi]->Branch("absrsp",&absrsp,"absrsp/F");
294     tPhi[iphi]->Branch("relrsp",&relrsp,"relrsp/F");
295 schiefer 1.3 tPhi[iphi]->Branch("weight",&weight,"weight/F");
296     }
297 schiefer 1.4 for (int iemf=0;iemf<nemfbins;iemf++) {
298     tEmf[iemf]=new TTree("tEmf","tEmf");
299     tEmf[iemf]->Branch("val", &emf, "val/F");
300     tEmf[iemf]->Branch("etcorr",&etcorr,"etcorr/F");
301     tEmf[iemf]->Branch("absrsp",&absrsp,"absrsp/F");
302     tEmf[iemf]->Branch("relrsp",&relrsp,"relrsp/F");
303     tEmf[iemf]->Branch("weight",&weight,"weight/F");
304     }
305 schiefer 1.3 for (int ieta=0;ieta<netabins;ieta++) {
306     tEtEtaEt[ieta] = new TTree*[netbins];
307     for (int iet=0;iet<netbins;iet++) {
308     tEtEtaEt[ieta][iet] = new TTree("tEtEtaEt","tEtEtaEt");
309 schiefer 1.4 tEtEtaEt[ieta][iet]->Branch("val", &et, "val/F");
310 schiefer 1.3 tEtEtaEt[ieta][iet]->Branch("etcorr",&etcorr,"etcorr/F");
311 schiefer 1.4 tEtEtaEt[ieta][iet]->Branch("absrsp",&absrsp,"absrsp/F");
312     tEtEtaEt[ieta][iet]->Branch("relrsp",&relrsp,"relrsp/F");
313 schiefer 1.3 tEtEtaEt[ieta][iet]->Branch("weight",&weight,"weight/F");
314     }
315     }
316 schiefer 1.4 for (int iemf=0;iemf<nemfbins;iemf++) {
317     tEtEmfEt[iemf] = new TTree*[netbins];
318     for (int iet=0;iet<netbins;iet++) {
319     tEtEmfEt[iemf][iet] = new TTree("tEtEmfEt","tEtEmfEt");
320     tEtEmfEt[iemf][iet]->Branch("val", &et, "val/F");
321     tEtEmfEt[iemf][iet]->Branch("etcorr",&etcorr,"etcorr/F");
322     tEtEmfEt[iemf][iet]->Branch("absrsp",&absrsp,"absrsp/F");
323     tEtEmfEt[iemf][iet]->Branch("relrsp",&relrsp,"relrsp/F");
324     tEtEmfEt[iemf][iet]->Branch("weight",&weight,"weight/F");
325     }
326     }
327     for (int iemf=0;iemf<nemfbins;iemf++) {
328     tEtEmfEtaEt[iemf] = new TTree**[netabins];
329     for (int ieta=0;ieta<netabins;ieta++) {
330     tEtEmfEtaEt[iemf][ieta] = new TTree*[netbins];
331     for (int iet=0;iet<netbins;iet++) {
332     tEtEmfEtaEt[iemf][ieta][iet] = new TTree("tEtEmfEtaEt","tEtEmfEtaEt");
333     tEtEmfEtaEt[iemf][ieta][iet]->Branch("val", &et, "val/F");
334     tEtEmfEtaEt[iemf][ieta][iet]->Branch("etcorr",&etcorr,"etcorr/F");
335     tEtEmfEtaEt[iemf][ieta][iet]->Branch("absrsp",&absrsp,"absrsp/F");
336     tEtEmfEtaEt[iemf][ieta][iet]->Branch("relrsp",&relrsp,"relrsp/F");
337     tEtEmfEtaEt[iemf][ieta][iet]->Branch("weight",&weight,"weight/F");
338     }
339     }
340     }
341 schiefer 1.1
342 schiefer 1.5 cout<<jetalgs[i]<<": trees created."<<endl;
343    
344 schiefer 1.1
345     // loop over all events
346     for (int ievt=0;ievt<nevts;ievt++) {
347    
348     t->GetEntry(el->GetEntry(ievt));
349    
350     // loop over all jets in the event
351     for (int ijt=0;ijt<njt;ijt++) {
352    
353     if (jtgendr[ijt]<0.0||jtgendr[ijt]>drmax) continue;
354    
355     if (ijt>0) continue;
356    
357 schiefer 1.3 et = jtet[ijt];
358     genet = jtgenet[ijt];
359     eta = (abseta) ? std::abs(jteta[ijt]) : jteta[ijt];
360     phi = jtphi[ijt];
361 schiefer 1.4 emf = jtemf[ijt];
362 schiefer 1.3 etcorr = (applyjes==0) ? et : et*jtjes[ijt][applyjes-1];
363 schiefer 1.1
364 schiefer 1.4 absrsp = genet - etcorr;
365     relrsp = etcorr/genet;
366 schiefer 1.2
367 schiefer 1.3 // genet
368     int igenet = hist::get_ibin(genet,etbins);
369     tGenEt[igenet]->Fill();
370 schiefer 1.1
371 schiefer 1.3 // et
372     int iet = hist::get_ibin(et,etbins);
373     tEt[iet]->Fill();
374 schiefer 1.1
375     // eta
376 schiefer 1.2 int ieta = hist::get_ibin(eta,etabins);
377 schiefer 1.3 tEta[ieta]->Fill();
378    
379 schiefer 1.1 // phi
380 schiefer 1.2 int iphi = hist::get_ibin(phi,phibins);
381 schiefer 1.3 tPhi[iphi]->Fill();
382    
383 schiefer 1.4 // emf
384     int iemf = hist::get_ibin(emf,emfbins);
385     tEmf[iemf]->Fill();
386    
387 schiefer 1.3 // et-eta
388     tEtEtaEt[ieta][iet]->Fill();
389 schiefer 1.1
390 schiefer 1.4 // et-emf
391     tEtEmfEt[iemf][iet]->Fill();
392    
393     // et-eta-emf
394     tEtEmfEtaEt[iemf][ieta][iet]->Fill();
395    
396 schiefer 1.1 } // jets
397     } // evts
398    
399     cout<<jetalgs[i]<<": all events processed / histograms filled."<<endl;
400 schiefer 1.3
401 schiefer 1.2
402 schiefer 1.3 // replace histograms
403 schiefer 1.4
404     replaceHistos(netbins,
405 schiefer 1.5 hGenEt,hEtCorrGenEt,hAbsRspGenEt,hRelRspGenEt,tGenEt);
406     replaceHistos(netbins,hEt,hEtCorrEt,hAbsRspEt,hRelRspEt,tEt);
407     replaceHistos(netabins,hEta,hEtCorrEta,hAbsRspEta,hRelRspEta,tEta);
408     replaceHistos(nphibins,hPhi,hEtCorrPhi,hAbsRspPhi,hRelRspPhi,tPhi);
409     replaceHistos(nemfbins,hEmf,hEtCorrEmf,hAbsRspEmf,hRelRspEmf,tEmf);
410 schiefer 1.4 replaceHistos(netabins,netbins,
411 schiefer 1.5 hEtEtaEt,hEtCorrEtaEt,hAbsRspEtaEt,hRelRspEtaEt,
412 schiefer 1.4 tEtEtaEt);
413     replaceHistos(nemfbins,netbins,
414 schiefer 1.5 hEtEmfEt,hEtCorrEmfEt,hAbsRspEmfEt,hRelRspEmfEt,
415 schiefer 1.4 tEtEmfEt);
416     replaceHistos(nemfbins,netabins,netbins,
417 schiefer 1.5 hEtEmfEtaEt,hEtCorrEmfEtaEt,
418     hAbsRspEmfEtaEt,hRelRspEmfEtaEt,
419 schiefer 1.4 tEtEmfEtaEt);
420 schiefer 1.2
421 schiefer 1.5 cout<<jetalgs[i]<<": all histograms replaced, and tree deleted."<<endl;
422    
423 schiefer 1.2 TDirectory* d = plotfile->mkdir(jetalgs[i].c_str());
424     d->cd();
425 schiefer 1.1
426 schiefer 1.3 // fit response histos
427 schiefer 1.5 fitHistos(hAbsRspGenEt,netbins,nsigma);
428     fitHistos(hRelRspGenEt,netbins,nsigma);
429 schiefer 1.3
430 schiefer 1.5 fitHistos(hAbsRspEt, netbins,nsigma);
431     fitHistos(hRelRspEt, netbins,nsigma);
432 schiefer 1.3
433 schiefer 1.5 fitHistos(hAbsRspEta, netabins,nsigma);
434     fitHistos(hRelRspEta, netabins,nsigma);
435 schiefer 1.3
436 schiefer 1.5 fitHistos(hAbsRspPhi, nphibins,nsigma);
437     fitHistos(hRelRspPhi, nphibins,nsigma);
438 schiefer 1.3
439 schiefer 1.5 fitHistos(hAbsRspEmf, nemfbins,nsigma);
440     fitHistos(hRelRspEmf, nemfbins,nsigma);
441 schiefer 1.4
442 schiefer 1.5 fitHistos(hAbsRspEtaEt,netabins,netbins,nsigma);
443     fitHistos(hRelRspEtaEt,netabins,netbins,nsigma);
444 schiefer 1.3
445 schiefer 1.5 fitHistos(hAbsRspEmfEt,nemfbins,netbins,nsigma);
446     fitHistos(hRelRspEmfEt,nemfbins,netbins,nsigma);
447 schiefer 1.4
448 schiefer 1.5 fitHistos(hAbsRspEmfEtaEt,nemfbins,netabins,netbins,nsigma);
449     fitHistos(hRelRspEmfEtaEt,nemfbins,netabins,netbins,nsigma);
450    
451     cout<<jetalgs[i]<<": fits done."<<endl;
452 schiefer 1.4
453 schiefer 1.3
454     // save histograms
455 schiefer 1.5 saveHistos(hGenEt, netbins);
456     saveHistos(hEtCorrGenEt,netbins);
457     saveHistos(hAbsRspGenEt,netbins);
458     saveHistos(hRelRspGenEt,netbins);
459    
460     saveHistos(hEt, netbins);
461     saveHistos(hEtCorrEt, netbins);
462     saveHistos(hAbsRspEt, netbins);
463     saveHistos(hRelRspEt, netbins);
464    
465     saveHistos(hEta, netabins);
466     saveHistos(hEtCorrEta, netabins);
467     saveHistos(hAbsRspEta, netabins);
468     saveHistos(hRelRspEta, netabins);
469    
470     saveHistos(hPhi, nphibins);
471     saveHistos(hEtCorrPhi, nphibins);
472     saveHistos(hAbsRspPhi, nphibins);
473     saveHistos(hRelRspPhi, nphibins);
474    
475     saveHistos(hEmf, nemfbins);
476     saveHistos(hEtCorrEmf, nemfbins);
477     saveHistos(hAbsRspEmf, nemfbins);
478     saveHistos(hRelRspEmf, nemfbins);
479    
480     saveHistos(hEtEtaEt, netabins,netbins);
481     saveHistos(hEtCorrEtaEt,netabins,netbins);
482     saveHistos(hAbsRspEtaEt,netabins,netbins);
483     saveHistos(hRelRspEtaEt,netabins,netbins);
484    
485     saveHistos(hEtEmfEt, nemfbins,netbins);
486     saveHistos(hEtCorrEmfEt,nemfbins,netbins);
487     saveHistos(hAbsRspEmfEt,nemfbins,netbins);
488     saveHistos(hRelRspEmfEt,nemfbins,netbins);
489    
490     saveHistos(hEtEmfEtaEt, nemfbins,netabins,netbins);
491     saveHistos(hEtCorrEmfEtaEt,nemfbins,netabins,netbins);
492     saveHistos(hAbsRspEmfEtaEt,nemfbins,netabins,netbins);
493     saveHistos(hRelRspEmfEtaEt,nemfbins,netabins,netbins);
494    
495     plotfile->Write();
496    
497     cout<<jetalgs[i]<<": histograms saved."<<endl;
498 schiefer 1.1 }
499    
500 schiefer 1.2 plotfile->Write();
501     plotfile->Close();
502     delete plotfile;
503 schiefer 1.1
504     return 0;
505     }
506    
507    
508    
509     ////////////////////////////////////////////////////////////////////////////////
510     // implementation of global functions
511     ////////////////////////////////////////////////////////////////////////////////
512    
513     //______________________________________________________________________________
514 schiefer 1.4 void replaceHistos(int nbins,
515     TH1F**& hVar,TH1F**& hEtCorr,TH1F**& hAbsRsp,TH1F**& hRelRsp,
516     TTree**& tVar)
517 schiefer 1.1 {
518 schiefer 1.3 TH1::SetDefaultSumw2();
519 schiefer 1.1 for (int ibin=0;ibin<nbins;ibin++) {
520 schiefer 1.4
521 schiefer 1.3 TH1F* h(0);
522 schiefer 1.4
523     if (tVar[ibin]->GetEntries()==0) continue;
524 schiefer 1.1
525 schiefer 1.3 tVar[ibin]->Project("h(50)","val","weight*(1)");
526     h = (TH1F*)gROOT->FindObject("h");
527     h->SetNameTitle(hVar[ibin]->GetName(),hVar[ibin]->GetTitle());
528     h->SetXTitle(hVar[ibin]->GetXaxis()->GetTitle());
529     h->SetYTitle(hVar[ibin]->GetYaxis()->GetTitle());
530     delete hVar[ibin];
531     hVar[ibin] = h;
532    
533     tVar[ibin]->Project("h(50)","etcorr","weight*(1)");
534     h = (TH1F*)gROOT->FindObject("h");
535     h->SetNameTitle(hEtCorr[ibin]->GetName(),hEtCorr[ibin]->GetTitle());
536     h->SetXTitle(hEtCorr[ibin]->GetXaxis()->GetTitle());
537     h->SetYTitle(hEtCorr[ibin]->GetYaxis()->GetTitle());
538     delete hEtCorr[ibin];
539     hEtCorr[ibin] = h;
540 schiefer 1.4
541     tVar[ibin]->Project("h(50)","absrsp","weight*(1)");
542     h = (TH1F*)gROOT->FindObject("h");
543     h->SetNameTitle(hAbsRsp[ibin]->GetName(),hAbsRsp[ibin]->GetTitle());
544     h->SetXTitle(hAbsRsp[ibin]->GetXaxis()->GetTitle());
545     h->SetYTitle(hAbsRsp[ibin]->GetYaxis()->GetTitle());
546     delete hAbsRsp[ibin];
547     hAbsRsp[ibin] = h;
548    
549     tVar[ibin]->Project("h(50)","relrsp","weight*(1)");
550     h = (TH1F*)gROOT->FindObject("h");
551     h->SetNameTitle(hRelRsp[ibin]->GetName(),hRelRsp[ibin]->GetTitle());
552     h->SetXTitle(hRelRsp[ibin]->GetXaxis()->GetTitle());
553     h->SetYTitle(hRelRsp[ibin]->GetYaxis()->GetTitle());
554     delete hRelRsp[ibin];
555     hRelRsp[ibin] = h;
556    
557 schiefer 1.3 delete tVar[ibin];
558 schiefer 1.1 }
559    
560 schiefer 1.3 delete [] tVar;
561 schiefer 1.1 }
562    
563 schiefer 1.4
564     //______________________________________________________________________________
565     void replaceHistos(int nbins1,int nbins2,
566     TH1F***& hVar,TH1F***& hEtCorr,TH1F***& hAbsRsp,TH1F***& hRelRsp,
567     TTree***& tVar)
568     {
569     for (int i1=0;i1<nbins1;i1++)
570     replaceHistos(nbins2,hVar[i1],hEtCorr[i1],hAbsRsp[i1],hRelRsp[i1],tVar[i1]);
571 schiefer 1.5 delete [] tVar;
572 schiefer 1.4 }
573    
574    
575     //______________________________________________________________________________
576     void replaceHistos(int nbins1,int nbins2,int nbins3,
577     TH1F****& hVar,TH1F****& hEtCorr,TH1F****&hAbsRsp,TH1F****&hRelRsp,
578     TTree****& tVar)
579     {
580     for (int i1=0;i1<nbins1;i1++)
581     replaceHistos(nbins2,nbins3,
582     hVar[i1],hEtCorr[i1],hAbsRsp[i1],hRelRsp[i1],tVar[i1]);
583 schiefer 1.5 delete [] tVar;
584 schiefer 1.4 }
585    
586    
587 schiefer 1.1 //______________________________________________________________________________
588 schiefer 1.3 void fitHistos(TH1F** histos,int n,double nsigma)
589 schiefer 1.1 {
590 schiefer 1.3 for (int i=0;i<n;i++) {
591     float nrm = histos[i]->Integral();
592     float mean = histos[i]->GetMean();
593     float rms = histos[i]->GetRMS();
594     TF1* f = new TF1("gaus","gaus",mean-nsigma*rms,mean+nsigma*rms);
595     f->SetParameter(0,nrm);
596     f->SetParameter(1,mean);
597     f->SetParameter(2,rms);
598     f->SetLineWidth(1);
599     f->SetLineColor(kRed);
600     histos[i]->Fit(f,"QR");
601 schiefer 1.1 }
602     }
603    
604    
605     //______________________________________________________________________________
606 schiefer 1.3 void fitHistos(TH1F*** histos,int n1,int n2,double nsigma)
607 schiefer 1.1 {
608 schiefer 1.3 for (int i1=0;i1<n1;i1++) fitHistos(histos[i1],n2,nsigma);
609 schiefer 1.1 }
610    
611    
612     //______________________________________________________________________________
613 schiefer 1.4 void fitHistos(TH1F**** histos,int n1,int n2,int n3,double nsigma)
614     {
615     for (int i1=0;i1<n1;i1++) fitHistos(histos[i1],n2,n3,nsigma);
616     }
617    
618    
619     //______________________________________________________________________________
620 schiefer 1.5 void saveHistos(TH1F**& histos,int n)
621 schiefer 1.1 {
622 schiefer 1.5 for (int i=0;i<n;i++) {
623     histos[i]->Write();
624     delete histos[i];
625     }
626     delete [] histos;
627 schiefer 1.1 }
628    
629    
630     //______________________________________________________________________________
631 schiefer 1.5 void saveHistos(TH1F***& histos,int n1,int n2)
632 schiefer 1.1 {
633 schiefer 1.3 for (int i1=0;i1<n1;i1++) saveHistos(histos[i1],n2);
634 schiefer 1.5 delete [] histos;
635 schiefer 1.1 }
636 schiefer 1.4
637    
638     //______________________________________________________________________________
639 schiefer 1.5 void saveHistos(TH1F****& histos,int n1,int n2,int n3)
640 schiefer 1.4 {
641     for (int i1=0;i1<n1;i1++) saveHistos(histos[i1],n2,n3);
642 schiefer 1.5 delete [] histos;
643 schiefer 1.4 }