ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/SchieferD/jetcalib/mcresponse_x.cpp
(Generate patch)

Comparing UserCode/SchieferD/jetcalib/mcresponse_x.cpp (file contents):
Revision 1.6 by schiefer, Fri Aug 31 13:56:04 2007 UTC vs.
Revision 1.13 by schiefer, Mon Sep 3 17:20:56 2007 UTC

# Line 10 | Line 10
10   #include "utils/cmdline.h"
11   #include "utils/hist.h"
12  
13 + #include <TSystem.h>
14   #include <TROOT.h>
15   #include <TRint.h>
16   #include <TFile.h>
# Line 30 | Line 31 | using namespace std;
31  
32  
33   ////////////////////////////////////////////////////////////////////////////////
33 – // defaults
34 – ////////////////////////////////////////////////////////////////////////////////
35 – string dflt_etbins ="8,12,15,20,25,30,45,65,90,120,180";
36 – 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";
37 – string dflt_phibins="-2.75,-2.25,-1.75.,-1.25,-.75,-.25,.25,.75,1.25,1.75,2.25,2.75";
38 – string dflt_emfbins="0.2,0.4,0.6,0.8";
39 –
40 –
41 – ////////////////////////////////////////////////////////////////////////////////
34   // declare global functions
35   ////////////////////////////////////////////////////////////////////////////////
36   void     replaceHistos(int nbins,
# Line 74 | Line 66 | int main(int argc,char**argv)
66    string         output    = cl.get_value<string> ("output", "mcrsp.root");
67    string         selection = cl.get_value<string> ("selection",   "njt>0");
68    float          drmax     = cl.get_value<float>  ("drmax",           0.3);
69 +  float          jtetmin   = cl.get_value<float>  ("jtetmin",         0.0);
70 +  float          jtetmax   = cl.get_value<float>  ("jtetmax",        1e06);
71 +  float          jtetamin  = cl.get_value<float>  ("jtetamin",        0.0);
72 +  float          jtetamax  = cl.get_value<float>  ("jtetamax",        9.9);
73    short          applyjes  = cl.get_value<short>  ("applyjes",          0);
74    string         treename  = cl.get_value<string> ("treename",        "t");
75    float          nsigma    = cl.get_value<float>  ("nsigma",          3.0);
76 <  vector<float>  etbins    = cl.get_vector<float> ("etbins",  dflt_etbins);
77 <  vector<float>  etabins   = cl.get_vector<float> ("etabins",dflt_etabins);
78 <  vector<float>  phibins   = cl.get_vector<float> ("phibins",dflt_phibins);
79 <  vector<float>  emfbins   = cl.get_vector<float> ("emfbins",dflt_emfbins);
76 >  vector<float>  etbins    = cl.get_vector<float> ("etbins",           "");
77 >  vector<float>  etabins   = cl.get_vector<float> ("etabins",          "");
78 >  vector<float>  phibins   = cl.get_vector<float> ("phibins",          "");
79 >  vector<float>  emfbins   = cl.get_vector<float> ("emfbins",          "");
80    bool           abseta    = cl.get_value<bool>   ("abseta",         true);
81    
82    if (!cl.check()) return 0;
83    cl.print();
84    
85 <
85 >  int netbins  = (etbins.size() >0) ? etbins.size() +1 : 0;
86 >  int netabins = (etabins.size()>0) ? etabins.size()+1 : 0;
87 >  int nphibins = (phibins.size()>0) ? phibins.size()+1 : 0;
88 >  int nemfbins = (emfbins.size()>0) ? emfbins.size()+1 : 0;
89 >  
90 >  if (netbins+netabins+nphibins+nemfbins==0) {
91 >    cout<<"Must specify bins: etbins, etabins, phibins, AND/OR emfbins!"<<endl;
92 >    return 0;
93 >  }
94 >  
95    // etabins
96 <  if (!abseta) {
96 >  if (netabins>0&&!abseta) {
97      int neta=(int)etabins.size();
98      std::reverse(etabins.begin(),etabins.end());
99      for (int ieta=neta-1;ieta>=0;ieta--) etabins.push_back(etabins[ieta]);
# Line 117 | Line 122 | int main(int argc,char**argv)
122    float jtgendr[100];
123    float jtjes[100][3];
124    
120 –  int netbins  = etbins.size()+1;
121 –  int netabins = etabins.size()+1;
122 –  int nphibins = phibins.size()+1;
123 –  int nemfbins = emfbins.size()+1;
124 –  
125    argc=3; argv[1]="-l"; argv[2]="-b";
126    TRint* app = new TRint(argv[0],&argc,argv); app->Argc();
127    
# Line 132 | Line 132 | int main(int argc,char**argv)
132    TFile* plotfile = new TFile(output.c_str(),"UPDATE");
133    
134    for (unsigned int i=0;i<input.size();++i) {
135 <    TFile* f = new TFile(input[i].c_str(),"UPDATE"); if (!f->IsOpen()) return 0;
135 >    TFile* f = new TFile(input[i].c_str(),"READ"); if (!f->IsOpen()) return 0;
136      TTree* t = (TTree*)f->Get(treename.c_str());
137      
138      TEventList* el = new TEventList("sel","sel");
# Line 146 | Line 146 | int main(int argc,char**argv)
146      t->SetBranchAddress("jtemf",  jtemf);
147      t->SetBranchAddress("jtgenet",jtgenet);
148      t->SetBranchAddress("jtgendr",jtgendr);
149 <
150 <    if (t->FindBranch("jtjes")) t->SetBranchAddress("jtjes", jtjes);
151 <    else
149 >    
150 >    vector<TBranch*> branches;
151 >    branches.push_back(t->GetBranch("weight"));
152 >    branches.push_back(t->GetBranch("njt"));
153 >    branches.push_back(t->GetBranch("jtet"));
154 >    branches.push_back(t->GetBranch("jteta"));
155 >    branches.push_back(t->GetBranch("jtphi"));
156 >    branches.push_back(t->GetBranch("jtemf"));
157 >    branches.push_back(t->GetBranch("jtgenet"));
158 >    branches.push_back(t->GetBranch("jtgendr"));
159 >    
160 >    if (t->FindBranch("jtjes")) {
161 >      t->SetBranchAddress("jtjes", jtjes);
162 >      branches.push_back(t->GetBranch("jtjes"));
163 >    }
164 >    else {
165        for (int i1=0;i1<100;i1++)
166          for (int i2=0;i2<3;i2++) jtjes[i1][i2]=1.0;
167 <    
167 >    }
168      string absrsp_xtitle = "E_{T}^{gen}-E_{T} [GeV]";
169      string relrsp_xtitle = "E_{T}^{rec}/E_{T}";
170  
158 –    TH1F** hGenEt       = hist::initHistos("GenEt",jetalgs[i],100,
159 –                                           "jet E_{T}^{gen} [GeV]","GenEt",etbins);
160 –    TH1F** hEtCorrGenEt = hist::initHistos("EtCorrGenEt",jetalgs[i],100,
161 –                                           "jet E_{T}^{gen} [GeV]","GenEt",etbins);
162 –    TH1F** hAbsRspGenEt = hist::initHistos("AbsRspGenEt",jetalgs[i],100,-50,150,
163 –                                           absrsp_xtitle,"GenEt",etbins);
164 –    TH1F** hRelRspGenEt = hist::initHistos("RelRspGenEt",jetalgs[i],100,0,2,
165 –                                           relrsp_xtitle,"GenEt",etbins);
166 –    
167 –    TH1F** hEt       = hist::initHistos("Et",jetalgs[i],100,
168 –                                        "jet E_{T} [GeV]","Et",etbins);
169 –    TH1F** hEtCorrEt = hist::initHistos("EtCorrEt",jetalgs[i],100,
170 –                                        "jet E_{T} [GeV]","Et",etbins);
171 –    TH1F** hAbsRspEt = hist::initHistos("AbsRspEt",jetalgs[i],100,-50,150,
172 –                                        absrsp_xtitle,"Et",etbins);
173 –    TH1F** hRelRspEt = hist::initHistos("RelRspEt",jetalgs[i],100,0,2,
174 –                                        relrsp_xtitle,"Et",etbins);
175 –    
176 –    TH1F** hEta       = hist::initHistos("Eta",jetalgs[i],100,
177 –                                         "jet #eta","Eta",etabins);
178 –    TH1F** hEtCorrEta = hist::initHistos("EtCorrEta",jetalgs[i],100,
179 –                                         "jet E_{T} [GeV]","Eta",etabins);
180 –    TH1F** hAbsRspEta = hist::initHistos("AbsRspEta",jetalgs[i],100,-50,150,
181 –                                         absrsp_xtitle,"Eta",etabins);
182 –    TH1F** hRelRspEta = hist::initHistos("RelRspEta",jetalgs[i],100,0,2,
183 –                                         relrsp_xtitle,"Eta",etabins);
184 –    
185 –    TH1F** hPhi       = hist::initHistos("Phi",jetalgs[i],100,
186 –                                         "jet #phi","Phi",phibins);
187 –    TH1F** hEtCorrPhi = hist::initHistos("EtCorrPhi",jetalgs[i],100,
188 –                                         "jet E_{T} [Gev]","Phi",phibins);
189 –    TH1F** hAbsRspPhi = hist::initHistos("AbsRspPhi",jetalgs[i],100,-50,150,
190 –                                         absrsp_xtitle,"Phi",phibins);
191 –    TH1F** hRelRspPhi = hist::initHistos("RelRspPhi",jetalgs[i],100,0,2,
192 –                                         relrsp_xtitle,"Phi",phibins);
193 –    
194 –    TH1F** hEmf       = hist::initHistos("Emf",jetalgs[i],100,
195 –                                         "jet emf","Emf",emfbins);
196 –    TH1F** hEtCorrEmf = hist::initHistos("EtCorrEmf",jetalgs[i],100,
197 –                                         "jet E_{T} [Gev]","Emf",emfbins);
198 –    TH1F** hAbsRspEmf = hist::initHistos("AbsRspEmf",jetalgs[i],100,-50,150,
199 –                                         absrsp_xtitle,"Emf",emfbins);
200 –    TH1F** hRelRspEmf = hist::initHistos("RelRspEmf",jetalgs[i],100,0,2,
201 –                                         relrsp_xtitle,"Emf",emfbins);
202 –    
203 –    TH1F*** hEtEtaEt     = hist::initHistos("Et",jetalgs[i],100,"jet E_{T} [GeV]",
204 –                                            "Eta",etabins,"Et",etbins);
205 –    TH1F*** hEtCorrEtaEt = hist::initHistos("EtCorrEtaEt",jetalgs[i],100,
206 –                                            "jet E_{T} [GeV]",
207 –                                            "Eta",etabins,"Et",etbins);
208 –    TH1F*** hAbsRspEtaEt = hist::initHistos("AbsRspEtaEt",jetalgs[i],100,-50,150,
209 –                                            absrsp_xtitle,
210 –                                            "Eta",etabins,"Et",etbins);
211 –    TH1F*** hRelRspEtaEt = hist::initHistos("RelRspEtaEt",jetalgs[i],100,0,2,
212 –                                            relrsp_xtitle,
213 –                                            "Eta",etabins,"Et",etbins);
214 –    
215 –    TH1F*** hEtEmfEt     = hist::initHistos("Et",jetalgs[i],100,"jet E_{T} [GeV]",
216 –                                            "Emf",emfbins,"Et",etbins);
217 –    TH1F*** hEtCorrEmfEt = hist::initHistos("EtCorrEmfEt",jetalgs[i],100,
218 –                                            "jet E_{T} [GeV]",
219 –                                            "Emf",emfbins,"Et",etbins);
220 –    TH1F*** hAbsRspEmfEt = hist::initHistos("AbsRspEmfEt",jetalgs[i],100,-50,150,
221 –                                            absrsp_xtitle,
222 –                                            "Emf",emfbins,"Et",etbins);
223 –    TH1F*** hRelRspEmfEt = hist::initHistos("RelRspEmfEt",jetalgs[i],100,0,2,
224 –                                            relrsp_xtitle,
225 –                                            "Emf",emfbins,"Et",etbins);
226 –    
227 –    TH1F**** hEtEmfEtaEt     = hist::initHistos("Et",jetalgs[i],100,"jet E_{T} [GeV]",
228 –                                                "Emf",emfbins,
229 –                                                "Eta",etabins,
230 –                                                "Et", etbins);
231 –    TH1F**** hEtCorrEmfEtaEt = hist::initHistos("EtCorrEmfEtaEt",jetalgs[i],100,
232 –                                                "jet E_{T} [GeV]",
233 –                                                "Emf",emfbins,
234 –                                                "Eta",etabins,
235 –                                                "Et", etbins);
236 –    TH1F**** hAbsRspEmfEtaEt = hist::initHistos("AbsRspEmfEtaEt",jetalgs[i],
237 –                                                100,-50,150,absrsp_xtitle,
238 –                                                "Emf",emfbins,
239 –                                                "Eta",etabins,
240 –                                                "Et", etbins);
241 –    TH1F**** hRelRspEmfEtaEt = hist::initHistos("RelRspEmfEtaEt",jetalgs[i],
242 –                                                100,0,2,relrsp_xtitle,
243 –                                                "Emf",emfbins,
244 –                                                "Eta",etabins,
245 –                                                "Et", etbins);
171      
172 +    TH1F** hGenEt(0);
173 +    TH1F** hEtCorrGenEt(0);
174 +    TH1F** hAbsRspGenEt(0);
175 +    TH1F** hRelRspGenEt(0);
176 +    
177 +    if (netbins>0) {
178 +      hGenEt       = hist::initHistos("GenEt",jetalgs[i],100,
179 +                                      "jet E_{T}^{gen} [GeV]",
180 +                                      "GenEt",etbins);
181 +      hEtCorrGenEt = hist::initHistos("EtCorrGenEt",jetalgs[i],100,
182 +                                      "jet E_{T}^{gen} [GeV]",
183 +                                      "GenEt",etbins);
184 +      hAbsRspGenEt = hist::initHistos("AbsRspGenEt",jetalgs[i],
185 +                                      100,-50,150,
186 +                                      absrsp_xtitle,
187 +                                      "GenEt",etbins);
188 +      hRelRspGenEt = hist::initHistos("RelRspGenEt",jetalgs[i],
189 +                                      100,0,2,
190 +                                      relrsp_xtitle,
191 +                                      "GenEt",etbins);
192 +    }
193 +    
194 +
195 +    TH1F** hEt(0);
196 +    TH1F** hEtCorrEt(0);
197 +    TH1F** hAbsRspEt(0);
198 +    TH1F** hRelRspEt(0);
199 +    
200 +    if (netbins>0) {
201 +      hEt       = hist::initHistos("Et",jetalgs[i],100,
202 +                                   "jet E_{T} [GeV]","Et",etbins);
203 +      hEtCorrEt = hist::initHistos("EtCorrEt",jetalgs[i],100,
204 +                                   "jet E_{T} [GeV]","Et",etbins);
205 +      hAbsRspEt = hist::initHistos("AbsRspEt",jetalgs[i],100,-50,150,
206 +                                   absrsp_xtitle,"Et",etbins);
207 +      hRelRspEt = hist::initHistos("RelRspEt",jetalgs[i],100,0,2,
208 +                                   relrsp_xtitle,"Et",etbins);
209 +    }
210 +    
211 +
212 +    TH1F** hEta(0);
213 +    TH1F** hEtCorrEta(0);
214 +    TH1F** hAbsRspEta(0);
215 +    TH1F** hRelRspEta(0);
216 +    
217 +    if (netabins>0) {
218 +      hEta       = hist::initHistos("Eta",jetalgs[i],100,
219 +                                    "jet #eta","Eta",etabins);
220 +      hEtCorrEta = hist::initHistos("EtCorrEta",jetalgs[i],100,
221 +                                    "jet E_{T} [GeV]","Eta",etabins);
222 +      hAbsRspEta = hist::initHistos("AbsRspEta",jetalgs[i],100,-50,150,
223 +                                    absrsp_xtitle,"Eta",etabins);
224 +      hRelRspEta = hist::initHistos("RelRspEta",jetalgs[i],100,0,2,
225 +                                    relrsp_xtitle,"Eta",etabins);
226 +    }
227 +    
228 +
229 +    TH1F** hPhi(0);
230 +    TH1F** hEtCorrPhi(0);
231 +    TH1F** hAbsRspPhi(0);
232 +    TH1F** hRelRspPhi(0);
233 +    
234 +    if (nphibins>0) {
235 +      hPhi       = hist::initHistos("Phi",jetalgs[i],100,
236 +                                    "jet #phi","Phi",phibins);
237 +      hEtCorrPhi = hist::initHistos("EtCorrPhi",jetalgs[i],100,
238 +                                    "jet E_{T} [Gev]","Phi",phibins);
239 +      hAbsRspPhi = hist::initHistos("AbsRspPhi",jetalgs[i],100,-50,150,
240 +                                    absrsp_xtitle,"Phi",phibins);
241 +      hRelRspPhi = hist::initHistos("RelRspPhi",jetalgs[i],100,0,2,
242 +                                    relrsp_xtitle,"Phi",phibins);
243 +    }
244 +    
245 +
246 +    TH1F** hEmf(0);
247 +    TH1F** hEtCorrEmf(0);
248 +    TH1F** hAbsRspEmf(0);
249 +    TH1F** hRelRspEmf(0);
250 +    
251 +    if (nemfbins>0) {
252 +      hEmf       = hist::initHistos("Emf",jetalgs[i],100,
253 +                                    "jet emf","Emf",emfbins);
254 +      hEtCorrEmf = hist::initHistos("EtCorrEmf",jetalgs[i],100,
255 +                                    "jet E_{T} [Gev]","Emf",emfbins);
256 +      hAbsRspEmf = hist::initHistos("AbsRspEmf",jetalgs[i],100,-50,150,
257 +                                    absrsp_xtitle,"Emf",emfbins);
258 +      hRelRspEmf = hist::initHistos("RelRspEmf",jetalgs[i],100,0,2,
259 +                                    relrsp_xtitle,"Emf",emfbins);
260 +    }
261 +    
262 +    
263 +    TH1F*** hEtEtaEt(0);
264 +    TH1F*** hEtCorrEtaEt(0);
265 +    TH1F*** hAbsRspEtaEt(0);
266 +    TH1F*** hRelRspEtaEt(0);
267 +    
268 +    if (netbins>0&&netabins>0) {
269 +      hEtEtaEt     = hist::initHistos("Et",jetalgs[i],100,"jet E_{T} [GeV]",
270 +                                      "Eta",etabins,"Et",etbins);
271 +      hEtCorrEtaEt = hist::initHistos("EtCorrEtaEt",jetalgs[i],100,
272 +                                      "jet E_{T} [GeV]",
273 +                                      "Eta",etabins,"Et",etbins);
274 +      hAbsRspEtaEt = hist::initHistos("AbsRspEtaEt",jetalgs[i],100,-50,150,
275 +                                      absrsp_xtitle,
276 +                                      "Eta",etabins,"Et",etbins);
277 +      hRelRspEtaEt = hist::initHistos("RelRspEtaEt",jetalgs[i],100,0,2,
278 +                                      relrsp_xtitle,
279 +                                      "Eta",etabins,"Et",etbins);
280 +    }
281 +
282 +
283 +    TH1F*** hEtEmfEt(0);
284 +    TH1F*** hEtCorrEmfEt(0);
285 +    TH1F*** hAbsRspEmfEt(0);
286 +    TH1F*** hRelRspEmfEt(0);
287 +    
288 +    if (netbins>0&&nemfbins>0) {
289 +      hEtEmfEt     = hist::initHistos("Et",jetalgs[i],100,"jet E_{T} [GeV]",
290 +                                      "Emf",emfbins,"Et",etbins);
291 +      hEtCorrEmfEt = hist::initHistos("EtCorrEmfEt",jetalgs[i],100,
292 +                                      "jet E_{T} [GeV]",
293 +                                      "Emf",emfbins,"Et",etbins);
294 +      hAbsRspEmfEt = hist::initHistos("AbsRspEmfEt",jetalgs[i],100,-50,150,
295 +                                      absrsp_xtitle,
296 +                                      "Emf",emfbins,"Et",etbins);
297 +      hRelRspEmfEt = hist::initHistos("RelRspEmfEt",jetalgs[i],100,0,2,
298 +                                      relrsp_xtitle,
299 +                                      "Emf",emfbins,"Et",etbins);
300 +    }
301 +    
302 +
303 +    TH1F**** hEtEmfEtaEt(0);
304 +    TH1F**** hEtCorrEmfEtaEt(0);
305 +    TH1F**** hAbsRspEmfEtaEt(0);
306 +    TH1F**** hRelRspEmfEtaEt(0);
307 +    
308 +    if (netbins>0&&netabins>0&&nemfbins>0) {
309 +      hEtEmfEtaEt     = hist::initHistos("Et",jetalgs[i],100,"jet E_{T} [GeV]",
310 +                                         "Emf",emfbins,
311 +                                         "Eta",etabins,
312 +                                         "Et", etbins);
313 +      hEtCorrEmfEtaEt = hist::initHistos("EtCorrEmfEtaEt",jetalgs[i],100,
314 +                                         "jet E_{T} [GeV]",
315 +                                         "Emf",emfbins,
316 +                                         "Eta",etabins,
317 +                                         "Et", etbins);
318 +      hAbsRspEmfEtaEt = hist::initHistos("AbsRspEmfEtaEt",jetalgs[i],
319 +                                         100,-50,150,absrsp_xtitle,
320 +                                         "Emf",emfbins,
321 +                                         "Eta",etabins,
322 +                                         "Et", etbins);
323 +      hRelRspEmfEtaEt = hist::initHistos("RelRspEmfEtaEt",jetalgs[i],
324 +                                         100,0,2,relrsp_xtitle,
325 +                                         "Emf",emfbins,
326 +                                         "Eta",etabins,
327 +                                         "Et", etbins);
328 +    }
329      
330      int nevts = el->GetN();
331      cout<<jetalgs[i]<<": "<<nevts<<" events selected."<<endl;
332      
333      // prepare intermediate trees
334      float genet,et,eta,phi,emf,etcorr,absrsp,relrsp;
335 <    TTree**   tGenEt      = new TTree*[netbins];
336 <    TTree**   tEt         = new TTree*[netbins];
337 <    TTree**   tEta        = new TTree*[netabins];
338 <    TTree**   tPhi        = new TTree*[nphibins];
339 <    TTree**   tEmf        = new TTree*[nemfbins];
340 <    TTree***  tEtEtaEt    = new TTree**[netabins];
341 <    TTree***  tEtEmfEt    = new TTree**[nemfbins];
342 <    TTree**** tEtEmfEtaEt = new TTree***[nemfbins];
335 >
336 >    TTree**   tGenEt(0);
337 >    TTree**   tEt(0);
338 >    TTree**   tEta(0);
339 >    TTree**   tPhi(0);
340 >    TTree**   tEmf(0);
341 >    TTree***  tEtEtaEt(0);
342 >    TTree***  tEtEmfEt(0);
343 >    TTree**** tEtEmfEtaEt(0);
344 >
345 >    gROOT->cd();
346 >    
347 >    if (netbins>0)  tGenEt = new TTree*[netbins];
348 >    if (netbins>0)  tEt    = new TTree*[netbins];
349 >    if (netabins>0) tEta   = new TTree*[netabins];
350 >    if (nphibins>0) tPhi   = new TTree*[nphibins];
351 >    if (nemfbins>0) tEmf   = new TTree*[nemfbins];
352 >    if (netbins>0&&netabins>0) tEtEtaEt = new TTree**[netabins];
353 >    if (netbins>0&&nemfbins>0) tEtEmfEt = new TTree**[nemfbins];
354 >    if (netbins>0&&netabins>0&&nemfbins>0) tEtEmfEtaEt = new TTree***[nemfbins];
355      
356      for (int iet=0;iet<netbins;iet++) {
357        
# Line 299 | Line 393 | int main(int argc,char**argv)
393        tEmf[iemf]->Branch("relrsp",&relrsp,"relrsp/F");
394        tEmf[iemf]->Branch("weight",&weight,"weight/F");
395      }
396 <    for (int ieta=0;ieta<netabins;ieta++) {
397 <      tEtEtaEt[ieta] = new TTree*[netbins];
304 <      for (int iet=0;iet<netbins;iet++) {
305 <        tEtEtaEt[ieta][iet] = new TTree("tEtEtaEt","tEtEtaEt");
306 <        tEtEtaEt[ieta][iet]->Branch("val",   &et,    "val/F");
307 <        tEtEtaEt[ieta][iet]->Branch("etcorr",&etcorr,"etcorr/F");
308 <        tEtEtaEt[ieta][iet]->Branch("absrsp",&absrsp,"absrsp/F");
309 <        tEtEtaEt[ieta][iet]->Branch("relrsp",&relrsp,"relrsp/F");
310 <        tEtEtaEt[ieta][iet]->Branch("weight",&weight,"weight/F");
311 <      }
312 <    }
313 <    for (int iemf=0;iemf<nemfbins;iemf++) {
314 <      tEtEmfEt[iemf] = new TTree*[netbins];
315 <      for (int iet=0;iet<netbins;iet++) {
316 <        tEtEmfEt[iemf][iet] = new TTree("tEtEmfEt","tEtEmfEt");
317 <        tEtEmfEt[iemf][iet]->Branch("val",   &et,    "val/F");
318 <        tEtEmfEt[iemf][iet]->Branch("etcorr",&etcorr,"etcorr/F");
319 <        tEtEmfEt[iemf][iet]->Branch("absrsp",&absrsp,"absrsp/F");
320 <        tEtEmfEt[iemf][iet]->Branch("relrsp",&relrsp,"relrsp/F");
321 <        tEtEmfEt[iemf][iet]->Branch("weight",&weight,"weight/F");
322 <      }
323 <    }
324 <    for (int iemf=0;iemf<nemfbins;iemf++) {
325 <      tEtEmfEtaEt[iemf] = new TTree**[netabins];
396 >
397 >    if (netbins>0) {
398        for (int ieta=0;ieta<netabins;ieta++) {
399 <        tEtEmfEtaEt[iemf][ieta] = new TTree*[netbins];
399 >        tEtEtaEt[ieta] = new TTree*[netbins];
400 >        for (int iet=0;iet<netbins;iet++) {
401 >          tEtEtaEt[ieta][iet] = new TTree("tEtEtaEt","tEtEtaEt");
402 >          tEtEtaEt[ieta][iet]->Branch("val",   &et,    "val/F");
403 >          tEtEtaEt[ieta][iet]->Branch("etcorr",&etcorr,"etcorr/F");
404 >          tEtEtaEt[ieta][iet]->Branch("absrsp",&absrsp,"absrsp/F");
405 >          tEtEtaEt[ieta][iet]->Branch("relrsp",&relrsp,"relrsp/F");
406 >          tEtEtaEt[ieta][iet]->Branch("weight",&weight,"weight/F");
407 >        }
408 >      }
409 >      for (int iemf=0;iemf<nemfbins;iemf++) {
410 >        tEtEmfEt[iemf] = new TTree*[netbins];
411          for (int iet=0;iet<netbins;iet++) {
412 <          tEtEmfEtaEt[iemf][ieta][iet] = new TTree("tEtEmfEtaEt","tEtEmfEtaEt");
413 <          tEtEmfEtaEt[iemf][ieta][iet]->Branch("val",   &et,    "val/F");
414 <          tEtEmfEtaEt[iemf][ieta][iet]->Branch("etcorr",&etcorr,"etcorr/F");
415 <          tEtEmfEtaEt[iemf][ieta][iet]->Branch("absrsp",&absrsp,"absrsp/F");
416 <          tEtEmfEtaEt[iemf][ieta][iet]->Branch("relrsp",&relrsp,"relrsp/F");
417 <          tEtEmfEtaEt[iemf][ieta][iet]->Branch("weight",&weight,"weight/F");
412 >          tEtEmfEt[iemf][iet] = new TTree("tEtEmfEt","tEtEmfEt");
413 >          tEtEmfEt[iemf][iet]->Branch("val",   &et,    "val/F");
414 >          tEtEmfEt[iemf][iet]->Branch("etcorr",&etcorr,"etcorr/F");
415 >          tEtEmfEt[iemf][iet]->Branch("absrsp",&absrsp,"absrsp/F");
416 >          tEtEmfEt[iemf][iet]->Branch("relrsp",&relrsp,"relrsp/F");
417 >          tEtEmfEt[iemf][iet]->Branch("weight",&weight,"weight/F");
418 >        }
419 >      }
420 >      if (netabins>0) {
421 >        for (int iemf=0;iemf<nemfbins;iemf++) {
422 >          tEtEmfEtaEt[iemf] = new TTree**[netabins];
423 >          for (int ieta=0;ieta<netabins;ieta++) {
424 >            tEtEmfEtaEt[iemf][ieta] = new TTree*[netbins];
425 >            for (int iet=0;iet<netbins;iet++) {
426 >              tEtEmfEtaEt[iemf][ieta][iet] = new TTree("tEtEmfEtaEt",
427 >                                                       "tEtEmfEtaEt");
428 >              tEtEmfEtaEt[iemf][ieta][iet]->Branch("val",   &et,    "val/F");
429 >              tEtEmfEtaEt[iemf][ieta][iet]->Branch("etcorr",&etcorr,"etcorr/F");
430 >              tEtEmfEtaEt[iemf][ieta][iet]->Branch("absrsp",&absrsp,"absrsp/F");
431 >              tEtEmfEtaEt[iemf][ieta][iet]->Branch("relrsp",&relrsp,"relrsp/F");
432 >              tEtEmfEtaEt[iemf][ieta][iet]->Branch("weight",&weight,"weight/F");
433 >            }
434 >          }
435          }
436        }
437      }
438      
439      cout<<jetalgs[i]<<": trees created."<<endl;
440  
341 –    
441      // loop over all events
442      for (int ievt=0;ievt<nevts;ievt++) {
443        
444 <      t->GetEntry(el->GetEntry(ievt));
445 <      
444 >      //t->GetEntry(el->GetEntry(ievt));
445 >      for (int ibrnch=0;ibrnch<(int)branches.size();ibrnch++)
446 >        branches[ibrnch]->GetEntry(el->GetEntry(ievt));
447 >
448        // loop over all jets in the event
449        for (int ijt=0;ijt<njt;ijt++) {
450          
451 <        if (jtgendr[ijt]<0.0||jtgendr[ijt]>drmax) continue;
452 <
453 <        // if (ijt>0) continue;
451 >        if (jtgendr[ijt]<0.0   ||jtgendr[ijt]>drmax)  continue;
452 >        if (jtet[ijt]<jtetmin  ||jtet[ijt]>jtetmax)   continue;
453 >        if (jteta[ijt]<jtetamin||jteta[ijt]>jtetamax) continue;
454  
455          et      = jtet[ijt];
456          genet   = jtgenet[ijt];
# Line 364 | Line 465 | int main(int argc,char**argv)
465  
466          if (TMath::IsNaN(emf)) emf = 0.0;
467          
468 <        
469 <        // genet
470 <        int igenet = hist::get_ibin(genet,etbins);
471 <        tGenEt[igenet]->Fill();
472 <        
473 <        // et
474 <        int iet = hist::get_ibin(et,etbins);
475 <        tEt[iet]->Fill();
476 <        
477 <        // eta
478 <        int ieta = hist::get_ibin(eta,etabins);
479 <        tEta[ieta]->Fill();
480 <        
481 <        // phi
482 <        int iphi = hist::get_ibin(phi,phibins);
483 <        tPhi[iphi]->Fill();
484 <        
384 <        // emf
385 <        int iemf = hist::get_ibin(emf,emfbins);
386 <        tEmf[iemf]->Fill();
387 <        
388 <        // et-eta
389 <        tEtEtaEt[ieta][iet]->Fill();
390 <        
391 <        // et-emf
392 <        tEtEmfEt[iemf][iet]->Fill();
393 <        
394 <        // et-eta-emf
395 <        tEtEmfEtaEt[iemf][ieta][iet]->Fill();
468 >        int igenet = (netbins>0)  ? hist::get_ibin(genet,etbins) : -1;
469 >        int iet    = (netbins>0)  ? hist::get_ibin(et,etbins)    : -1;
470 >        int ieta   = (netabins>0) ? hist::get_ibin(eta,etabins)  : -1;
471 >        int iphi   = (nphibins>0) ? hist::get_ibin(phi,phibins)  : -1;
472 >        int iemf   = (nemfbins>0) ? hist::get_ibin(emf,emfbins)  : -1;
473 >        
474 >        if (netbins>0)  tGenEt[igenet]->Fill();
475 >        if (netbins>0)  tEt[iet]->Fill();
476 >        if (netabins>0) tEta[ieta]->Fill();
477 >        if (nphibins>0) tPhi[iphi]->Fill();
478 >        if (nemfbins>0) tEmf[iemf]->Fill();
479 >
480 >        if (netbins>0&&netabins>0) tEtEtaEt[ieta][iet]->Fill();
481 >        if (netbins>0&&nemfbins>0) tEtEmfEt[iemf][iet]->Fill();
482 >
483 >        if (netbins>0&&netabins>0&&nemfbins>0)
484 >          tEtEmfEtaEt[iemf][ieta][iet]->Fill();
485          
486        } // jets
487      } // evts
# Line 422 | Line 511 | int main(int argc,char**argv)
511      cout<<jetalgs[i]<<": all histograms replaced, and tree deleted."<<endl;
512  
513      
514 +    plotfile->cd();
515      TDirectory* d = plotfile->GetDirectory(jetalgs[i].c_str());
516      if (0!=d) plotfile->rmdir(jetalgs[i].c_str());
517      d = plotfile->mkdir(jetalgs[i].c_str());
# Line 430 | Line 520 | int main(int argc,char**argv)
520      // fit response histos
521      fitHistos(hAbsRspGenEt,netbins,nsigma);
522      fitHistos(hRelRspGenEt,netbins,nsigma);
523 <
523 >    
524      fitHistos(hAbsRspEt,   netbins,nsigma);
525      fitHistos(hRelRspEt,   netbins,nsigma);
526  
# Line 439 | Line 529 | int main(int argc,char**argv)
529  
530      fitHistos(hAbsRspPhi,  nphibins,nsigma);
531      fitHistos(hRelRspPhi,  nphibins,nsigma);
532 <    
532 >
533      fitHistos(hAbsRspEmf,  nemfbins,nsigma);
534      fitHistos(hRelRspEmf,  nemfbins,nsigma);
535 <    
535 >
536      fitHistos(hAbsRspEtaEt,netabins,netbins,nsigma);
537      fitHistos(hRelRspEtaEt,netabins,netbins,nsigma);
538  
# Line 496 | Line 586 | int main(int argc,char**argv)
586      saveHistos(hAbsRspEmfEtaEt,nemfbins,netabins,netbins);
587      saveHistos(hRelRspEmfEtaEt,nemfbins,netabins,netbins);
588  
499 –    plotfile->Write();
500 –  
589      cout<<jetalgs[i]<<": histograms saved."<<endl;
590    }
591    
592    plotfile->Write();
593    plotfile->Close();
594    delete plotfile;
595 +
596 +  //app->Run();
597    
598    return 0;
599   }
# Line 519 | Line 609 | void replaceHistos(int nbins,
609                     TH1F**& hVar,TH1F**& hEtCorr,TH1F**& hAbsRsp,TH1F**& hRelRsp,
610                     TTree**& tVar)
611   {
612 +  if (nbins==0) return;
613 +  
614    TH1::SetDefaultSumw2();
615    for (int ibin=0;ibin<nbins;ibin++) {
616      
# Line 570 | Line 662 | void replaceHistos(int nbins1,int nbins2
662                     TH1F***& hVar,TH1F***& hEtCorr,TH1F***& hAbsRsp,TH1F***& hRelRsp,
663                     TTree***& tVar)
664   {
665 +  if (nbins1==0) return;
666 +  
667    for (int i1=0;i1<nbins1;i1++)
668      replaceHistos(nbins2,hVar[i1],hEtCorr[i1],hAbsRsp[i1],hRelRsp[i1],tVar[i1]);
669    delete [] tVar;
# Line 581 | Line 675 | void replaceHistos(int nbins1,int nbins2
675                     TH1F****& hVar,TH1F****& hEtCorr,TH1F****&hAbsRsp,TH1F****&hRelRsp,
676                     TTree****& tVar)
677   {
678 +  if (nbins1==0) return;
679    for (int i1=0;i1<nbins1;i1++)
680      replaceHistos(nbins2,nbins3,
681                    hVar[i1],hEtCorr[i1],hAbsRsp[i1],hRelRsp[i1],tVar[i1]);
# Line 601 | Line 696 | void fitHistos(TH1F** histos,int n,doubl
696      f->SetParameter(2,rms);
697      f->SetLineWidth(1);
698      f->SetLineColor(kRed);
699 <    histos[i]->Fit(f,"QR");
699 >    histos[i]->Fit(f,"QR"); // Avoid TCanvas?
700    }
701   }
702  
# Line 609 | Line 704 | void fitHistos(TH1F** histos,int n,doubl
704   //______________________________________________________________________________
705   void fitHistos(TH1F*** histos,int n1,int n2,double nsigma)
706   {
707 +  if (n1==0||n2==0) return;
708    for (int i1=0;i1<n1;i1++) fitHistos(histos[i1],n2,nsigma);
709   }
710  
# Line 616 | Line 712 | void fitHistos(TH1F*** histos,int n1,int
712   //______________________________________________________________________________
713   void fitHistos(TH1F**** histos,int n1,int n2,int n3,double nsigma)
714   {
715 +  if (n1==0||n2==0||n3==0) return;
716    for (int i1=0;i1<n1;i1++) fitHistos(histos[i1],n2,n3,nsigma);
717   }
718  
# Line 623 | Line 720 | void fitHistos(TH1F**** histos,int n1,in
720   //______________________________________________________________________________
721   void saveHistos(TH1F**& histos,int n)
722   {
723 +  if (n==0) return;
724    for (int i=0;i<n;i++) {
725      histos[i]->Write();
726      delete histos[i];
# Line 634 | Line 732 | void saveHistos(TH1F**& histos,int n)
732   //______________________________________________________________________________
733   void saveHistos(TH1F***& histos,int n1,int n2)
734   {
735 +  if (n1==0) return;
736    for (int i1=0;i1<n1;i1++) saveHistos(histos[i1],n2);
737    delete [] histos;
738   }
# Line 642 | Line 741 | void saveHistos(TH1F***& histos,int n1,i
741   //______________________________________________________________________________
742   void saveHistos(TH1F****& histos,int n1,int n2,int n3)
743   {
744 +  if (n1==0) return;
745    for (int i1=0;i1<n1;i1++) saveHistos(histos[i1],n2,n3);
746    delete [] histos;
747   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines