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.2 by schiefer, Tue Aug 21 06:42:45 2007 UTC vs.
Revision 1.12 by schiefer, Mon Sep 3 14:36:27 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 17 | Line 18
18   #include <TEventList.h>
19   #include <TH1F.h>
20   #include <TF1.h>
21 < #include <TGraphErrors.h>
21 < #include <TMultiGraph.h>
22 < #include <TCanvas.h>
23 < #include <TLegend.h>
21 > #include <TMath.h>
22  
23   #include <iostream>
24   #include <iomanip>
27 #include <fstream>
25   #include <cmath>
26   #include <string>
27   #include <vector>
# Line 34 | Line 31 | using namespace std;
31  
32  
33   ////////////////////////////////////////////////////////////////////////////////
37 // defaults
38 ////////////////////////////////////////////////////////////////////////////////
39 string dflt_ptbins ="7,9,12,15,20,25,30,45,65,90,120";
40 string dflt_etabins=".2,.4,.7,1.,1.3,1.5,1.75,2.,2.3,2.5,2.7,2.9,3.2";
41 //string dflt_etabins=".3,.6,.9,1.2,1.5,1.8,2.1,2.4,2.7";
42 string dflt_phibins="-2.75,-2.25,-1.75.,-1.25,-.75,-.25,.25,.75,1.25,1.75,2.25,2.75";
43
44
45 ////////////////////////////////////////////////////////////////////////////////
34   // declare global functions
35   ////////////////////////////////////////////////////////////////////////////////
36 + void     replaceHistos(int nbins,
37 +                       TH1F**& hVar,TH1F**& hEtCorr,TH1F**& hAbsRsp,TH1F**& hRelRsp,
38 +                       TTree**& tVar);
39 + void     replaceHistos(int nbins1,int nbins2,
40 +                       TH1F***&hVar,TH1F***&hEtCorr,TH1F***&hAbsRsp,TH1F***&hRelRsp,
41 +                       TTree***& tVar);
42 + void     replaceHistos(int nbins1,int nbins2,int nbins3,
43 +                       TH1F****&hVar,TH1F****&hEtCorr,
44 +                       TH1F****&hAbsRsp,TH1F****&hRelRsp,
45 +                       TTree****& tVar);
46 + void     fitHistos    (TH1F**    histos,int n,double nsigma);
47 + void     fitHistos    (TH1F***   histos,int n1,int n2,double nsigma);
48 + void     fitHistos    (TH1F****  histos,int n1,int n2,int n3,double nsigma);
49 + void     saveHistos   (TH1F**&   histos,int n);
50 + void     saveHistos   (TH1F***&  histos,int n1,int n2);
51 + void     saveHistos   (TH1F****& histos,int n1,int n2,int n3);
52  
49 void     fillRspAndSgm(TGraphErrors* gRsp,TGraphErrors* gSgm,float nsigma,
50                       const vector<float>& bins,TH1F** hVar,TH1F** hPt,TH1F** hRsp);
51 TCanvas* plotVarHistos(const vector<float>& bins,TH1F** hVar,bool logy=false);
52 TCanvas* plotRspHistos(const vector<float>& bins,TH1F** hRsp);
53 TCanvas* plotResponse (const vector<string>& jetalgs,short applyjes,
54                       const vector<TGraphErrors*> gRsp,
55                       const string& varname,const string& xtitle);
56 TCanvas* plotResponse (const string& jetalg,short applyjes,bool fitcorr,
57                       const vector<float>& bins1,TGraphErrors** gRsp,
58                       const string& varname1,const string vartitle1,
59                       const string& varname2,const string& xtitle);
60 TCanvas* plotResolution(const vector<string>& jetalgs,
61                        const vector<TGraphErrors*> gSgm,
62                        const string& varname,const string& xtitle);
63 TCanvas* plotResolution(const string& jetalg,const vector<float>& bins1,
64                        TGraphErrors** gSgm,
65                        const string& varname1,const string vartitle1,
66                        const string& varname2,const string& xtitle);
67
68 float  get_xmax(TGraph* g);
53  
54  
55   ////////////////////////////////////////////////////////////////////////////////
# Line 79 | Line 63 | int main(int argc,char**argv)
63    cl.parse(argc,argv);
64    
65    vector<string> input     = cl.get_vector<string>("input");
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 +  short          applyjes  = cl.get_value<short>  ("applyjes",          0);
70    string         treename  = cl.get_value<string> ("treename",        "t");
71    float          nsigma    = cl.get_value<float>  ("nsigma",          3.0);
72 <  short          applyjes  = cl.get_value<short>  ("applyjes",          0);
73 <  bool           fitcorr   = cl.get_value<bool>   ("fitcorr",       false);
74 <  vector<float>  ptbins    = cl.get_vector<float> ("ptbins",  dflt_ptbins);
75 <  vector<float>  etabins   = cl.get_vector<float> ("etabins",dflt_etabins);
90 <  vector<float>  phibins   = cl.get_vector<float> ("phibins",dflt_phibins);
72 >  vector<float>  etbins    = cl.get_vector<float> ("etbins",           "");
73 >  vector<float>  etabins   = cl.get_vector<float> ("etabins",          "");
74 >  vector<float>  phibins   = cl.get_vector<float> ("phibins",          "");
75 >  vector<float>  emfbins   = cl.get_vector<float> ("emfbins",          "");
76    bool           abseta    = cl.get_value<bool>   ("abseta",         true);
77    
78    if (!cl.check()) return 0;
79    cl.print();
80    
81 <
81 >  int netbins  = (etbins.size() >0) ? etbins.size() +1 : 0;
82 >  int netabins = (etabins.size()>0) ? etabins.size()+1 : 0;
83 >  int nphibins = (phibins.size()>0) ? phibins.size()+1 : 0;
84 >  int nemfbins = (emfbins.size()>0) ? emfbins.size()+1 : 0;
85 >  
86 >  if (netbins+netabins+nphibins+nemfbins==0) {
87 >    cout<<"Must specify bins: etbins, etabins, phibins, AND/OR emfbins!"<<endl;
88 >    return 0;
89 >  }
90 >  
91    // etabins
92 <  if (!abseta) {
92 >  if (netabins>0&&!abseta) {
93      int neta=(int)etabins.size();
94      std::reverse(etabins.begin(),etabins.end());
95      for (int ieta=neta-1;ieta>=0;ieta--) etabins.push_back(etabins[ieta]);
# Line 116 | Line 110 | int main(int argc,char**argv)
110    // prepare branch variables and histograms / graphs
111    float weight;
112    char  njt;
113 <  float jtpt[100];
113 >  float jtet[100];
114    float jteta[100];
115    float jtphi[100];
116 <  float jtgenpt[100];
116 >  float jtemf[100];
117 >  float jtgenet[100];
118    float jtgendr[100];
119    float jtjes[100][3];
120    
121 <  int nptbins  = ptbins.size()+1;
122 <  int netabins = etabins.size()+1;
128 <  int nphibins = phibins.size()+1;
129 <  
130 <  vector<TH1F**>         hGenPt;
131 <  vector<TH1F**>         hGenPtPtCorr;
132 <  vector<TH1F**>         hRspVsGenPt;
133 <  vector<TH1F**>         hPt;
134 <  vector<TH1F**>         hPtCorr;
135 <  vector<TH1F**>         hRspVsPt;
136 <  vector<TH1F**>         hEta;
137 <  vector<TH1F**>         hEtaPtCorr;
138 <  vector<TH1F**>         hRspVsEta;
139 <  vector<TH1F**>         hPhi;
140 <  vector<TH1F**>         hPhiPtCorr;
141 <  vector<TH1F**>         hRspVsPhi;
142 <  
143 <  vector<TH1F***>        hEtaPtPt;
144 <  vector<TH1F***>        hEtaPtPtCorr;
145 <  vector<TH1F***>        hRspVsEtaPt;
146 <  
147 <  vector<TGraphErrors*>  gRspVsGenPt;
148 <  vector<TGraphErrors*>  gSgmVsGenPt;
149 <  vector<TGraphErrors*>  gRspVsPt;
150 <  vector<TGraphErrors*>  gSgmVsPt;
151 <  vector<TGraphErrors*>  gRspVsEta;
152 <  vector<TGraphErrors*>  gSgmVsEta;
153 <  vector<TGraphErrors*>  gRspVsPhi;
154 <  vector<TGraphErrors*>  gSgmVsPhi;
155 <  
156 <  vector<TGraphErrors**> gRspVsEtaPt;
157 <  vector<TGraphErrors**> gSgmVsEtaPt;
158 <  
159 <  argc=1; //argv[1]="-b";
160 <  TRint* app = new TRint(argv[0],&argc,argv);
121 >  argc=3; argv[1]="-l"; argv[2]="-b";
122 >  TRint* app = new TRint(argv[0],&argc,argv); app->Argc();
123    
124    
125    //
126    // loop over all input files
127    //
128 <  TFile* plotfile = new TFile("mcresponse.root","RECREATE");
128 >  TFile* plotfile = new TFile(output.c_str(),"UPDATE");
129    
130    for (unsigned int i=0;i<input.size();++i) {
131      TFile* f = new TFile(input[i].c_str(),"READ"); if (!f->IsOpen()) return 0;
# Line 174 | Line 136 | int main(int argc,char**argv)
136  
137      t->SetBranchAddress("weight", &weight);
138      t->SetBranchAddress("njt",    &njt);
139 <    t->SetBranchAddress("jtpt",   jtpt);
139 >    t->SetBranchAddress("jtet",   jtet);
140      t->SetBranchAddress("jteta",  jteta);
141      t->SetBranchAddress("jtphi",  jtphi);
142 <    t->SetBranchAddress("jtgenpt",jtgenpt);
142 >    t->SetBranchAddress("jtemf",  jtemf);
143 >    t->SetBranchAddress("jtgenet",jtgenet);
144      t->SetBranchAddress("jtgendr",jtgendr);
145 <
146 <    if (t->FindBranch("jtjes")) t->SetBranchAddress("jtjes", jtjes);
147 <    else
145 >    
146 >    vector<TBranch*> branches;
147 >    branches.push_back(t->GetBranch("weight"));
148 >    branches.push_back(t->GetBranch("njt"));
149 >    branches.push_back(t->GetBranch("jtet"));
150 >    branches.push_back(t->GetBranch("jteta"));
151 >    branches.push_back(t->GetBranch("jtphi"));
152 >    branches.push_back(t->GetBranch("jtemf"));
153 >    branches.push_back(t->GetBranch("jtgenet"));
154 >    branches.push_back(t->GetBranch("jtgendr"));
155 >    
156 >    if (t->FindBranch("jtjes")) {
157 >      t->SetBranchAddress("jtjes", jtjes);
158 >      branches.push_back(t->GetBranch("jtjes"));
159 >    }
160 >    else {
161        for (int i1=0;i1<100;i1++)
162          for (int i2=0;i2<3;i2++) jtjes[i1][i2]=1.0;
163 +    }
164 +    string absrsp_xtitle = "E_{T}^{gen}-E_{T} [GeV]";
165 +    string relrsp_xtitle = "E_{T}^{rec}/E_{T}";
166 +
167 +    
168 +    TH1F** hGenEt(0);
169 +    TH1F** hEtCorrGenEt(0);
170 +    TH1F** hAbsRspGenEt(0);
171 +    TH1F** hRelRspGenEt(0);
172 +    
173 +    if (netbins>0) {
174 +      hGenEt       = hist::initHistos("GenEt",jetalgs[i],100,
175 +                                      "jet E_{T}^{gen} [GeV]",
176 +                                      "GenEt",etbins);
177 +      hEtCorrGenEt = hist::initHistos("EtCorrGenEt",jetalgs[i],100,
178 +                                      "jet E_{T}^{gen} [GeV]",
179 +                                      "GenEt",etbins);
180 +      hAbsRspGenEt = hist::initHistos("AbsRspGenEt",jetalgs[i],
181 +                                      100,-50,150,
182 +                                      absrsp_xtitle,
183 +                                      "GenEt",etbins);
184 +      hRelRspGenEt = hist::initHistos("RelRspGenEt",jetalgs[i],
185 +                                      100,0,2,
186 +                                      relrsp_xtitle,
187 +                                      "GenEt",etbins);
188 +    }
189 +    
190 +
191 +    TH1F** hEt(0);
192 +    TH1F** hEtCorrEt(0);
193 +    TH1F** hAbsRspEt(0);
194 +    TH1F** hRelRspEt(0);
195 +    
196 +    if (netbins>0) {
197 +      hEt       = hist::initHistos("Et",jetalgs[i],100,
198 +                                   "jet E_{T} [GeV]","Et",etbins);
199 +      hEtCorrEt = hist::initHistos("EtCorrEt",jetalgs[i],100,
200 +                                   "jet E_{T} [GeV]","Et",etbins);
201 +      hAbsRspEt = hist::initHistos("AbsRspEt",jetalgs[i],100,-50,150,
202 +                                   absrsp_xtitle,"Et",etbins);
203 +      hRelRspEt = hist::initHistos("RelRspEt",jetalgs[i],100,0,2,
204 +                                   relrsp_xtitle,"Et",etbins);
205 +    }
206 +    
207 +
208 +    TH1F** hEta(0);
209 +    TH1F** hEtCorrEta(0);
210 +    TH1F** hAbsRspEta(0);
211 +    TH1F** hRelRspEta(0);
212 +    
213 +    if (netabins>0) {
214 +      hEta       = hist::initHistos("Eta",jetalgs[i],100,
215 +                                    "jet #eta","Eta",etabins);
216 +      hEtCorrEta = hist::initHistos("EtCorrEta",jetalgs[i],100,
217 +                                    "jet E_{T} [GeV]","Eta",etabins);
218 +      hAbsRspEta = hist::initHistos("AbsRspEta",jetalgs[i],100,-50,150,
219 +                                    absrsp_xtitle,"Eta",etabins);
220 +      hRelRspEta = hist::initHistos("RelRspEta",jetalgs[i],100,0,2,
221 +                                    relrsp_xtitle,"Eta",etabins);
222 +    }
223      
224 +
225 +    TH1F** hPhi(0);
226 +    TH1F** hEtCorrPhi(0);
227 +    TH1F** hAbsRspPhi(0);
228 +    TH1F** hRelRspPhi(0);
229 +    
230 +    if (nphibins>0) {
231 +      hPhi       = hist::initHistos("Phi",jetalgs[i],100,
232 +                                    "jet #phi","Phi",phibins);
233 +      hEtCorrPhi = hist::initHistos("EtCorrPhi",jetalgs[i],100,
234 +                                    "jet E_{T} [Gev]","Phi",phibins);
235 +      hAbsRspPhi = hist::initHistos("AbsRspPhi",jetalgs[i],100,-50,150,
236 +                                    absrsp_xtitle,"Phi",phibins);
237 +      hRelRspPhi = hist::initHistos("RelRspPhi",jetalgs[i],100,0,2,
238 +                                    relrsp_xtitle,"Phi",phibins);
239 +    }
240 +    
241 +
242 +    TH1F** hEmf(0);
243 +    TH1F** hEtCorrEmf(0);
244 +    TH1F** hAbsRspEmf(0);
245 +    TH1F** hRelRspEmf(0);
246 +    
247 +    if (nemfbins>0) {
248 +      hEmf       = hist::initHistos("Emf",jetalgs[i],100,
249 +                                    "jet emf","Emf",emfbins);
250 +      hEtCorrEmf = hist::initHistos("EtCorrEmf",jetalgs[i],100,
251 +                                    "jet E_{T} [Gev]","Emf",emfbins);
252 +      hAbsRspEmf = hist::initHistos("AbsRspEmf",jetalgs[i],100,-50,150,
253 +                                    absrsp_xtitle,"Emf",emfbins);
254 +      hRelRspEmf = hist::initHistos("RelRspEmf",jetalgs[i],100,0,2,
255 +                                    relrsp_xtitle,"Emf",emfbins);
256 +    }
257      
189    string rsp_xtitle = "E_{T}^{gen}-E_{T}^{rec} [GeV]";
258      
259 <    hGenPt      .push_back(hist::initHistos("JetGenPt",jetalgs[i],
260 <                                            4000,0,4000,"jet p_{T}^{gen} [GeV]",
261 <                                            "GenPt",ptbins));
262 <    hGenPtPtCorr.push_back(hist::initHistos("JetPtCorr",jetalgs[i],
263 <                                            4000,0,4000,"jet p_{T}^{gen} [GeV]",
264 <                                            "GenPt",ptbins));
265 <    hRspVsGenPt .push_back(hist::initHistos("JetRspGenPt",jetalgs[i],
266 <                                            100,-100,100,rsp_xtitle,
267 <                                            "GenPt",ptbins));
268 <    
269 <    hPt         .push_back(hist::initHistos("JetPt",jetalgs[i],
270 <                                            4000,0,4000,"jet p_{T} [GeV]",
271 <                                            "Pt",ptbins));
272 <    hPtCorr     .push_back(hist::initHistos("JetPtCorr",jetalgs[i],
273 <                                            4000,0,4000,"jet p_{T} [GeV]",
274 <                                            "Pt",ptbins));
275 <    hRspVsPt    .push_back(hist::initHistos("JetRspPt",jetalgs[i],
276 <                                            100,-100,100,rsp_xtitle,
277 <                                            "Pt",ptbins));
278 <    
279 <    hEta        .push_back(hist::initHistos("JetEta",jetalgs[i],
280 <                                            100,0,5,"jet #eta","Eta",etabins));
281 <    hEtaPtCorr  .push_back(hist::initHistos("JetPtCorr",jetalgs[i],
282 <                                            4000,0,4000,"jet p_{T} [GeV]",
283 <                                            "Eta",etabins));
284 <    hRspVsEta   .push_back(hist::initHistos("JetRspEta",jetalgs[i],
285 <                                            100,-100,100,rsp_xtitle,"Eta",etabins));
286 <    
287 <    hPhi        .push_back(hist::initHistos("JetPhi",jetalgs[i],
288 <                                            100,-3.2,3.2,"jet #phi","Phi",phibins));
289 <    hPhiPtCorr  .push_back(hist::initHistos("JetPtCorr",jetalgs[i],
290 <                                            4000,0,4000,"jet p_{T} [Gev]",
291 <                                            "Phi",phibins));
292 <    hRspVsPhi   .push_back(hist::initHistos("JetRspPhi",jetalgs[i],
293 <                                            100,-100,100,rsp_xtitle,"Phi",phibins));
294 <    
295 <    hEtaPtPt    .push_back(hist::initHistos("JetPt",jetalgs[i],
296 <                                            4000,0,4000,"jet p_{T} [GeV]",
229 <                                            "Eta",etabins,"Pt",ptbins));
230 <    hEtaPtPtCorr.push_back(hist::initHistos("JetPtCorr",jetalgs[i],
231 <                                            4000,0,4000,"jet p_{T} [GeV]",
232 <                                            "Eta",etabins,"Pt",ptbins));
233 <    hRspVsEtaPt .push_back(hist::initHistos("JetRspEtaPt",jetalgs[i],
234 <                                            100,-100,100,rsp_xtitle,
235 <                                            "Eta",etabins,"Pt",ptbins));
259 >    TH1F*** hEtEtaEt(0);
260 >    TH1F*** hEtCorrEtaEt(0);
261 >    TH1F*** hAbsRspEtaEt(0);
262 >    TH1F*** hRelRspEtaEt(0);
263 >    
264 >    if (netbins>0&&netabins>0) {
265 >      hEtEtaEt     = hist::initHistos("Et",jetalgs[i],100,"jet E_{T} [GeV]",
266 >                                      "Eta",etabins,"Et",etbins);
267 >      hEtCorrEtaEt = hist::initHistos("EtCorrEtaEt",jetalgs[i],100,
268 >                                      "jet E_{T} [GeV]",
269 >                                      "Eta",etabins,"Et",etbins);
270 >      hAbsRspEtaEt = hist::initHistos("AbsRspEtaEt",jetalgs[i],100,-50,150,
271 >                                      absrsp_xtitle,
272 >                                      "Eta",etabins,"Et",etbins);
273 >      hRelRspEtaEt = hist::initHistos("RelRspEtaEt",jetalgs[i],100,0,2,
274 >                                      relrsp_xtitle,
275 >                                      "Eta",etabins,"Et",etbins);
276 >    }
277 >
278 >
279 >    TH1F*** hEtEmfEt(0);
280 >    TH1F*** hEtCorrEmfEt(0);
281 >    TH1F*** hAbsRspEmfEt(0);
282 >    TH1F*** hRelRspEmfEt(0);
283 >    
284 >    if (netbins>0&&nemfbins>0) {
285 >      hEtEmfEt     = hist::initHistos("Et",jetalgs[i],100,"jet E_{T} [GeV]",
286 >                                      "Emf",emfbins,"Et",etbins);
287 >      hEtCorrEmfEt = hist::initHistos("EtCorrEmfEt",jetalgs[i],100,
288 >                                      "jet E_{T} [GeV]",
289 >                                      "Emf",emfbins,"Et",etbins);
290 >      hAbsRspEmfEt = hist::initHistos("AbsRspEmfEt",jetalgs[i],100,-50,150,
291 >                                      absrsp_xtitle,
292 >                                      "Emf",emfbins,"Et",etbins);
293 >      hRelRspEmfEt = hist::initHistos("RelRspEmfEt",jetalgs[i],100,0,2,
294 >                                      relrsp_xtitle,
295 >                                      "Emf",emfbins,"Et",etbins);
296 >    }
297      
298 +
299 +    TH1F**** hEtEmfEtaEt(0);
300 +    TH1F**** hEtCorrEmfEtaEt(0);
301 +    TH1F**** hAbsRspEmfEtaEt(0);
302 +    TH1F**** hRelRspEmfEtaEt(0);
303 +    
304 +    if (netbins>0&&netabins>0&&nemfbins>0) {
305 +      hEtEmfEtaEt     = hist::initHistos("Et",jetalgs[i],100,"jet E_{T} [GeV]",
306 +                                         "Emf",emfbins,
307 +                                         "Eta",etabins,
308 +                                         "Et", etbins);
309 +      hEtCorrEmfEtaEt = hist::initHistos("EtCorrEmfEtaEt",jetalgs[i],100,
310 +                                         "jet E_{T} [GeV]",
311 +                                         "Emf",emfbins,
312 +                                         "Eta",etabins,
313 +                                         "Et", etbins);
314 +      hAbsRspEmfEtaEt = hist::initHistos("AbsRspEmfEtaEt",jetalgs[i],
315 +                                         100,-50,150,absrsp_xtitle,
316 +                                         "Emf",emfbins,
317 +                                         "Eta",etabins,
318 +                                         "Et", etbins);
319 +      hRelRspEmfEtaEt = hist::initHistos("RelRspEmfEtaEt",jetalgs[i],
320 +                                         100,0,2,relrsp_xtitle,
321 +                                         "Emf",emfbins,
322 +                                         "Eta",etabins,
323 +                                         "Et", etbins);
324 +    }
325      
326      int nevts = el->GetN();
327      cout<<jetalgs[i]<<": "<<nevts<<" events selected."<<endl;
328      
329 +    // prepare intermediate trees
330 +    float genet,et,eta,phi,emf,etcorr,absrsp,relrsp;
331 +
332 +    TTree**   tGenEt(0);
333 +    TTree**   tEt(0);
334 +    TTree**   tEta(0);
335 +    TTree**   tPhi(0);
336 +    TTree**   tEmf(0);
337 +    TTree***  tEtEtaEt(0);
338 +    TTree***  tEtEmfEt(0);
339 +    TTree**** tEtEmfEtaEt(0);
340 +
341 +    gROOT->cd();
342 +    
343 +    if (netbins>0)  tGenEt = new TTree*[netbins];
344 +    if (netbins>0)  tEt    = new TTree*[netbins];
345 +    if (netabins>0) tEta   = new TTree*[netabins];
346 +    if (nphibins>0) tPhi   = new TTree*[nphibins];
347 +    if (nemfbins>0) tEmf   = new TTree*[nemfbins];
348 +    if (netbins>0&&netabins>0) tEtEtaEt = new TTree**[netabins];
349 +    if (netbins>0&&nemfbins>0) tEtEmfEt = new TTree**[nemfbins];
350 +    if (netbins>0&&netabins>0&&nemfbins>0) tEtEmfEtaEt = new TTree***[nemfbins];
351      
352 +    for (int iet=0;iet<netbins;iet++) {
353 +      
354 +      tGenEt[iet]=new TTree("tGenEt","tGenEt");
355 +      tGenEt[iet]->Branch("val",   &genet, "val/F");
356 +      tGenEt[iet]->Branch("etcorr",&etcorr,"etcorr/F");
357 +      tGenEt[iet]->Branch("absrsp",&absrsp,"absrsp/F");
358 +      tGenEt[iet]->Branch("relrsp",&relrsp,"relrsp/F");
359 +      tGenEt[iet]->Branch("weight",&weight,"weight/F");
360 +
361 +      tEt[iet]=new TTree("tEt","tEt");
362 +      tEt[iet]->Branch("val",   &et,    "val/F");
363 +      tEt[iet]->Branch("etcorr",&etcorr,"etcorr/F");
364 +      tEt[iet]->Branch("absrsp",&absrsp,"absrsp/F");
365 +      tEt[iet]->Branch("relrsp",&relrsp,"relrsp/F");
366 +      tEt[iet]->Branch("weight",&weight,"weight/F");
367 +    }
368 +    for (int ieta=0;ieta<netabins;ieta++) {
369 +      tEta[ieta]=new TTree("tEta","tEta");
370 +      tEta[ieta]->Branch("val",   &eta,   "val/F");
371 +      tEta[ieta]->Branch("etcorr",&etcorr,"etcorr/F");
372 +      tEta[ieta]->Branch("absrsp",&absrsp,"absrsp/F");
373 +      tEta[ieta]->Branch("relrsp",&relrsp,"relrsp/F");
374 +      tEta[ieta]->Branch("weight",&weight,"weight/F");
375 +    }
376 +    for (int iphi=0;iphi<nphibins;iphi++) {
377 +      tPhi[iphi]=new TTree("tPhi","tPhi");
378 +      tPhi[iphi]->Branch("val",   &phi,   "val/F");
379 +      tPhi[iphi]->Branch("etcorr",&etcorr,"etcorr/F");
380 +      tPhi[iphi]->Branch("absrsp",&absrsp,"absrsp/F");
381 +      tPhi[iphi]->Branch("relrsp",&relrsp,"relrsp/F");
382 +      tPhi[iphi]->Branch("weight",&weight,"weight/F");
383 +    }
384 +    for (int iemf=0;iemf<nemfbins;iemf++) {
385 +      tEmf[iemf]=new TTree("tEmf","tEmf");
386 +      tEmf[iemf]->Branch("val",   &emf,   "val/F");
387 +      tEmf[iemf]->Branch("etcorr",&etcorr,"etcorr/F");
388 +      tEmf[iemf]->Branch("absrsp",&absrsp,"absrsp/F");
389 +      tEmf[iemf]->Branch("relrsp",&relrsp,"relrsp/F");
390 +      tEmf[iemf]->Branch("weight",&weight,"weight/F");
391 +    }
392 +
393 +    if (netbins>0) {
394 +      for (int ieta=0;ieta<netabins;ieta++) {
395 +        tEtEtaEt[ieta] = new TTree*[netbins];
396 +        for (int iet=0;iet<netbins;iet++) {
397 +          tEtEtaEt[ieta][iet] = new TTree("tEtEtaEt","tEtEtaEt");
398 +          tEtEtaEt[ieta][iet]->Branch("val",   &et,    "val/F");
399 +          tEtEtaEt[ieta][iet]->Branch("etcorr",&etcorr,"etcorr/F");
400 +          tEtEtaEt[ieta][iet]->Branch("absrsp",&absrsp,"absrsp/F");
401 +          tEtEtaEt[ieta][iet]->Branch("relrsp",&relrsp,"relrsp/F");
402 +          tEtEtaEt[ieta][iet]->Branch("weight",&weight,"weight/F");
403 +        }
404 +      }
405 +      for (int iemf=0;iemf<nemfbins;iemf++) {
406 +        tEtEmfEt[iemf] = new TTree*[netbins];
407 +        for (int iet=0;iet<netbins;iet++) {
408 +          tEtEmfEt[iemf][iet] = new TTree("tEtEmfEt","tEtEmfEt");
409 +          tEtEmfEt[iemf][iet]->Branch("val",   &et,    "val/F");
410 +          tEtEmfEt[iemf][iet]->Branch("etcorr",&etcorr,"etcorr/F");
411 +          tEtEmfEt[iemf][iet]->Branch("absrsp",&absrsp,"absrsp/F");
412 +          tEtEmfEt[iemf][iet]->Branch("relrsp",&relrsp,"relrsp/F");
413 +          tEtEmfEt[iemf][iet]->Branch("weight",&weight,"weight/F");
414 +        }
415 +      }
416 +      if (netabins>0) {
417 +        for (int iemf=0;iemf<nemfbins;iemf++) {
418 +          tEtEmfEtaEt[iemf] = new TTree**[netabins];
419 +          for (int ieta=0;ieta<netabins;ieta++) {
420 +            tEtEmfEtaEt[iemf][ieta] = new TTree*[netbins];
421 +            for (int iet=0;iet<netbins;iet++) {
422 +              tEtEmfEtaEt[iemf][ieta][iet] = new TTree("tEtEmfEtaEt",
423 +                                                       "tEtEmfEtaEt");
424 +              tEtEmfEtaEt[iemf][ieta][iet]->Branch("val",   &et,    "val/F");
425 +              tEtEmfEtaEt[iemf][ieta][iet]->Branch("etcorr",&etcorr,"etcorr/F");
426 +              tEtEmfEtaEt[iemf][ieta][iet]->Branch("absrsp",&absrsp,"absrsp/F");
427 +              tEtEmfEtaEt[iemf][ieta][iet]->Branch("relrsp",&relrsp,"relrsp/F");
428 +              tEtEmfEtaEt[iemf][ieta][iet]->Branch("weight",&weight,"weight/F");
429 +            }
430 +          }
431 +        }
432 +      }
433 +    }
434 +    
435 +    cout<<jetalgs[i]<<": trees created."<<endl;
436 +
437      // loop over all events
438      for (int ievt=0;ievt<nevts;ievt++) {
439        
440 <      t->GetEntry(el->GetEntry(ievt));
441 <      
440 >      //t->GetEntry(el->GetEntry(ievt));
441 >      for (int ibrnch=0;ibrnch<(int)branches.size();ibrnch++)
442 >        branches[ibrnch]->GetEntry(el->GetEntry(ievt));
443 >
444        // loop over all jets in the event
445        for (int ijt=0;ijt<njt;ijt++) {
446          
447          if (jtgendr[ijt]<0.0||jtgendr[ijt]>drmax) continue;
448  
449 <        if (ijt>0) continue;
449 >        //if (ijt>ijtmax) continue;
450 >        
451 >        et      = jtet[ijt];
452 >        genet   = jtgenet[ijt];
453 >        eta     = (abseta) ? std::abs(jteta[ijt]) : jteta[ijt];
454 >        phi     = jtphi[ijt];
455 >        emf     = jtemf[ijt];
456  
457 <        float pt      = jtpt[ijt];
255 <        float genpt   = jtgenpt[ijt];
256 <        float eta     = (abseta) ? std::abs(jteta[ijt]) : jteta[ijt];
257 <        float phi     = jtphi[ijt];
457 >        etcorr  = (applyjes==0) ? et : et*jtjes[ijt][applyjes-1];
458          
459 <        float ptcorr  = (applyjes==0) ? pt : pt*jtjes[ijt][applyjes-1];
460 <        float deltapt = genpt - ptcorr;
459 >        absrsp  = genet - etcorr;
460 >        relrsp  = etcorr/genet;
461  
462 <        // genpt
263 <        int igenpt = hist::get_ibin(genpt,ptbins);
264 <        hGenPt[i][igenpt]      ->Fill(genpt,weight);
265 <        //hGenPtPtCorr[i][igenpt]->Fill(ptcorr,weight);
266 <        hGenPtPtCorr[i][igenpt]->Fill(genpt,weight);
267 <        hRspVsGenPt[i][igenpt] ->Fill(deltapt,weight);
462 >        if (TMath::IsNaN(emf)) emf = 0.0;
463          
464 <        // pt
465 <        int ipt = hist::get_ibin(pt,ptbins);
466 <        hPt[i][ipt]     ->Fill(pt,weight);
467 <        hPtCorr[i][ipt] ->Fill(ptcorr,weight);
468 <        hRspVsPt[i][ipt]->Fill(deltapt,weight);
464 >        int igenet = (netbins>0)  ? hist::get_ibin(genet,etbins) : -1;
465 >        int iet    = (netbins>0)  ? hist::get_ibin(et,etbins)    : -1;
466 >        int ieta   = (netabins>0) ? hist::get_ibin(eta,etabins)  : -1;
467 >        int iphi   = (nphibins>0) ? hist::get_ibin(phi,phibins)  : -1;
468 >        int iemf   = (nemfbins>0) ? hist::get_ibin(emf,emfbins)  : -1;
469          
470 <        // eta
471 <        int ieta = hist::get_ibin(eta,etabins);
472 <        hEta[i][ieta]      ->Fill(eta,    weight);
473 <        hEtaPtCorr[i][ieta]->Fill(ptcorr, weight);
474 <        hRspVsEta[i][ieta] ->Fill(deltapt,weight);
475 <                
476 <        // phi
477 <        int iphi = hist::get_ibin(phi,phibins);
478 <        hPhi[i][iphi]      ->Fill(phi,    weight);
479 <        hPhiPtCorr[i][iphi]->Fill(ptcorr, weight);
480 <        hRspVsPhi[i][iphi] ->Fill(deltapt,weight);
286 <                
287 <        // pt-eta
288 <        hEtaPtPt[i][ieta][ipt]    ->Fill(pt,     weight);
289 <        hEtaPtPtCorr[i][ieta][ipt]->Fill(ptcorr, weight);
290 <        hRspVsEtaPt[i][ieta][ipt] ->Fill(deltapt,weight);
291 <                
470 >        if (netbins>0)  tGenEt[igenet]->Fill();
471 >        if (netbins>0)  tEt[iet]->Fill();
472 >        if (netabins>0) tEta[ieta]->Fill();
473 >        if (nphibins>0) tPhi[iphi]->Fill();
474 >        if (nemfbins>0) tEmf[iemf]->Fill();
475 >
476 >        if (netbins>0&&netabins>0) tEtEtaEt[ieta][iet]->Fill();
477 >        if (netbins>0&&nemfbins>0) tEtEmfEt[iemf][iet]->Fill();
478 >
479 >        if (netbins>0&&netabins>0&&nemfbins>0)
480 >          tEtEmfEtaEt[iemf][ieta][iet]->Fill();
481          
482        } // jets
483      } // evts
484      
485      cout<<jetalgs[i]<<": all events processed / histograms filled."<<endl;
297
486      
487 <    TDirectory* d = plotfile->mkdir(jetalgs[i].c_str());
487 >
488 >    // replace histograms
489 >
490 >    replaceHistos(netbins,
491 >                  hGenEt,hEtCorrGenEt,hAbsRspGenEt,hRelRspGenEt,tGenEt);
492 >    replaceHistos(netbins,hEt,hEtCorrEt,hAbsRspEt,hRelRspEt,tEt);
493 >    replaceHistos(netabins,hEta,hEtCorrEta,hAbsRspEta,hRelRspEta,tEta);
494 >    replaceHistos(nphibins,hPhi,hEtCorrPhi,hAbsRspPhi,hRelRspPhi,tPhi);
495 >    replaceHistos(nemfbins,hEmf,hEtCorrEmf,hAbsRspEmf,hRelRspEmf,tEmf);
496 >    replaceHistos(netabins,netbins,
497 >                  hEtEtaEt,hEtCorrEtaEt,hAbsRspEtaEt,hRelRspEtaEt,
498 >                  tEtEtaEt);
499 >    replaceHistos(nemfbins,netbins,
500 >                  hEtEmfEt,hEtCorrEmfEt,hAbsRspEmfEt,hRelRspEmfEt,
501 >                  tEtEmfEt);
502 >    replaceHistos(nemfbins,netabins,netbins,
503 >                  hEtEmfEtaEt,hEtCorrEmfEtaEt,
504 >                  hAbsRspEmfEtaEt,hRelRspEmfEtaEt,
505 >                  tEtEmfEtaEt);
506 >    
507 >    cout<<jetalgs[i]<<": all histograms replaced, and tree deleted."<<endl;
508 >
509 >    
510 >    plotfile->cd();
511 >    TDirectory* d = plotfile->GetDirectory(jetalgs[i].c_str());
512 >    if (0!=d) plotfile->rmdir(jetalgs[i].c_str());
513 >    d = plotfile->mkdir(jetalgs[i].c_str());
514      d->cd();
515      
516 +    // fit response histos
517 +    fitHistos(hAbsRspGenEt,netbins,nsigma);
518 +    fitHistos(hRelRspGenEt,netbins,nsigma);
519      
520 <    // response / resolution vs *genpt*
521 <    gRspVsGenPt.push_back(new TGraphErrors(nptbins));
305 <    gSgmVsGenPt.push_back(new TGraphErrors(nptbins));
306 <    gRspVsGenPt.back()->SetName(("RspVsGenPt_"+jetalgs[i]).c_str());
307 <    gSgmVsGenPt.back()->SetName(("SgmVsGenPt_"+jetalgs[i]).c_str());
308 <    fillRspAndSgm(gRspVsGenPt.back(),gSgmVsGenPt.back(),nsigma,
309 <                  ptbins,hGenPt[i],hGenPtPtCorr[i],hRspVsGenPt[i]);
310 <    
311 <    // response / resolution vs *pt*
312 <    gRspVsPt.push_back(new TGraphErrors(nptbins));
313 <    gSgmVsPt.push_back(new TGraphErrors(nptbins));
314 <    gRspVsPt.back()->SetName(("RspVsPt_"+jetalgs[i]).c_str());
315 <    gSgmVsPt.back()->SetName(("SgmVsPt_"+jetalgs[i]).c_str());
316 <    fillRspAndSgm(gRspVsPt.back(),gSgmVsPt.back(),nsigma,
317 <                  ptbins,hPt[i],hPtCorr[i],hRspVsPt[i]);
318 <    
319 <    // response / resolution vs *eta*
320 <    gRspVsEta.push_back(new TGraphErrors(netabins));
321 <    gSgmVsEta.push_back(new TGraphErrors(netabins));
322 <    gRspVsEta.back()->SetName(("RspVsEta_"+jetalgs[i]).c_str());
323 <    gSgmVsEta.back()->SetName(("SgmVsEta_"+jetalgs[i]).c_str());
324 <    fillRspAndSgm(gRspVsEta.back(),gSgmVsEta.back(),nsigma,
325 <                  etabins,hEta[i],hEtaPtCorr[i],hRspVsEta[i]);
326 <    
327 <    // response / resolution vs *phi*
328 <    gRspVsPhi.push_back(new TGraphErrors(nphibins));
329 <    gSgmVsPhi.push_back(new TGraphErrors(nphibins));
330 <    gRspVsPhi.back()->SetName(("RspVsPhi_"+jetalgs[i]).c_str());
331 <    gSgmVsPhi.back()->SetName(("SgmVsPhi_"+jetalgs[i]).c_str());
332 <    fillRspAndSgm(gRspVsPhi.back(),gSgmVsPhi.back(),nsigma,
333 <                  phibins,hPhi[i],hPhiPtCorr[i],hRspVsPhi[i]);
334 <    
335 <    // response / resolution vs *pt-eta*
336 <    TGraphErrors** g;
337 <    g=new TGraphErrors*[netabins];
338 <    for (int ieta=0;ieta<netabins;ieta++) g[ieta] = new TGraphErrors(nptbins);
339 <    gRspVsEtaPt.push_back(g);
340 <    g=new TGraphErrors*[netabins];
341 <    for (int ieta=0;ieta<netabins;ieta++) g[ieta] = new TGraphErrors(nptbins);
342 <    gSgmVsEtaPt.push_back(g);
343 <    for (int ieta=0;ieta<netabins;ieta++) {
344 <      gRspVsEtaPt[i][ieta]->SetName(("RspVsPhi_"+
345 <                                     hist::get_bin_name("Eta",ieta,etabins)+
346 <                                     "_"+jetalgs[i]).c_str());
347 <      gSgmVsEtaPt[i][ieta]->SetName(("SgmVsPhi_"+
348 <                                     hist::get_bin_name("Eta",ieta,etabins)+
349 <                                     "_"+jetalgs[i]).c_str());
350 <      fillRspAndSgm(gRspVsEtaPt[i][ieta],gSgmVsEtaPt[i][ieta],nsigma,ptbins,
351 <                    hEtaPtPt[i][ieta],hEtaPtPtCorr[i][ieta],hRspVsEtaPt[i][ieta]);
352 <    }
353 <  }
354 <  
355 <  
356 <  TCanvas* c = (TCanvas*)gROOT->FindObject("c1");
357 <  if (0!=c) delete c;
358 <  
520 >    fitHistos(hAbsRspEt,   netbins,nsigma);
521 >    fitHistos(hRelRspEt,   netbins,nsigma);
522  
523 <  //
524 <  // make plots
525 <  //
526 <  for (unsigned int ialg=0;ialg<jetalgs.size();ialg++) {
527 <    
528 <    TDirectory* d = plotfile->GetDirectory(jetalgs[ialg].c_str());
529 <    d->cd();
523 >    fitHistos(hAbsRspEta,  netabins,nsigma);
524 >    fitHistos(hRelRspEta,  netabins,nsigma);
525 >
526 >    fitHistos(hAbsRspPhi,  nphibins,nsigma);
527 >    fitHistos(hRelRspPhi,  nphibins,nsigma);
528 >
529 >    fitHistos(hAbsRspEmf,  nemfbins,nsigma);
530 >    fitHistos(hRelRspEmf,  nemfbins,nsigma);
531 >
532 >    fitHistos(hAbsRspEtaEt,netabins,netbins,nsigma);
533 >    fitHistos(hRelRspEtaEt,netabins,netbins,nsigma);
534 >
535 >    fitHistos(hAbsRspEmfEt,nemfbins,netbins,nsigma);
536 >    fitHistos(hRelRspEmfEt,nemfbins,netbins,nsigma);
537 >
538 >    fitHistos(hAbsRspEmfEtaEt,nemfbins,netabins,netbins,nsigma);
539 >    fitHistos(hRelRspEmfEtaEt,nemfbins,netabins,netbins,nsigma);
540 >
541 >    cout<<jetalgs[i]<<": fits done."<<endl;
542  
368    // genpt
369    TCanvas* cGenPt       = plotVarHistos(ptbins,hGenPt[ialg],true);
370    TCanvas* cGenPtPtCorr = plotVarHistos(ptbins,hGenPtPtCorr[ialg],true);
371    TCanvas* cRspGenPt    = plotRspHistos(ptbins,hRspVsGenPt[ialg]);
372    
373    delete cGenPt;
374    delete cGenPtPtCorr;
375    delete cRspGenPt;
376    
377    // pt
378    TCanvas* cPt     = plotVarHistos(ptbins,hPt[ialg],true);
379    TCanvas* cPtCorr = plotVarHistos(ptbins,hPtCorr[ialg],true);
380    TCanvas* cRspPt  = plotRspHistos(ptbins,hRspVsPt[ialg]);
381    
382    delete cPt;
383    delete cPtCorr;
384    delete cRspPt;
385    
386    // eta
387    TCanvas* cEta       = plotVarHistos(etabins,hEta[ialg]);
388    TCanvas* cEtaPtCorr = plotVarHistos(etabins,hEtaPtCorr[ialg],true);
389    TCanvas* cRspEta    = plotRspHistos(etabins,hRspVsEta[ialg]);
390    
391    delete cEta;
392    delete cEtaPtCorr;
393    delete cRspEta;
394    
395    // phi
396    TCanvas* cPhi       = plotVarHistos(phibins,hPhi[ialg]);
397    TCanvas* cPhiPtCorr = plotVarHistos(phibins,hPhiPtCorr[ialg],true);
398    TCanvas* cRspPhi    = plotRspHistos(phibins,hRspVsPhi[ialg]);
399
400    delete cPhi;
401    delete cPhiPtCorr;
402    delete cRspPhi;
403    
404    // eta-pt
405    for (int ieta=0;ieta<netabins;ieta++) {
406      TCanvas* cEtaPtPt     = plotVarHistos(ptbins,hEtaPtPt[ialg][ieta],true);
407      TCanvas* cEtaPtPtCorr = plotVarHistos(ptbins,hEtaPtPtCorr[ialg][ieta],true);
408      TCanvas* cRspEtaPt    = plotRspHistos(ptbins,hRspVsEtaPt[ialg][ieta]);
409      
410      delete cEtaPtPt;
411      delete cEtaPtPtCorr;
412      delete cRspEtaPt;
413    }
543      
544 <  }
545 <  
546 <  plotfile->cd();
544 >    // save histograms
545 >    saveHistos(hGenEt,      netbins);
546 >    saveHistos(hEtCorrGenEt,netbins);
547 >    saveHistos(hAbsRspGenEt,netbins);
548 >    saveHistos(hRelRspGenEt,netbins);
549 >
550 >    saveHistos(hEt,         netbins);
551 >    saveHistos(hEtCorrEt,   netbins);
552 >    saveHistos(hAbsRspEt,   netbins);
553 >    saveHistos(hRelRspEt,   netbins);
554  
555 <  // response & resolution vs genpT
556 <  TCanvas* cRspVsGenPt=plotResponse(jetalgs,applyjes,gRspVsGenPt,
557 <                                    "GenPt","jet p_{T}^{gen} [GeV]");
558 <  TCanvas* cSgmVsGenPt=plotResolution(jetalgs,gSgmVsGenPt,
423 <                                      "GenPt","jet p_{T}^{gen} [GeV]");
424 <  cRspVsGenPt->cd();
425 <  cSgmVsGenPt->cd();
426 <
427 <  // response & resolution vs pT
428 <  TCanvas* cRspVsPt  = plotResponse(jetalgs,applyjes,gRspVsPt,"Pt","jet p_{T} [GeV]");
429 <  TCanvas* cSgmVsPt  = plotResolution(jetalgs,gSgmVsPt,"Pt","jet p_{T} [GeV]");
430 <  cRspVsPt->cd();
431 <  cSgmVsPt->cd();
432 <
433 <  // response & resolution vs eta
434 <  TCanvas* cRspVsEta = plotResponse(jetalgs,applyjes,gRspVsEta,"Eta","jet #eta");
435 <  TCanvas* cSgmVsEta = plotResolution(jetalgs,gSgmVsEta,"Eta","jet #eta");
436 <  cRspVsEta->cd();
437 <  cSgmVsEta->cd();
438 <
439 <  // response & resolution vs phi
440 <  TCanvas* cRspVsPhi = plotResponse(jetalgs,applyjes,gRspVsPhi,"Phi","jet #phi");
441 <  TCanvas* cSgmVsPhi = plotResolution(jetalgs,gSgmVsPhi,"Phi","jet #phi");
442 <  cRspVsPhi->cd();
443 <  cSgmVsPhi->cd();
555 >    saveHistos(hEta,        netabins);
556 >    saveHistos(hEtCorrEta,  netabins);
557 >    saveHistos(hAbsRspEta,  netabins);
558 >    saveHistos(hRelRspEta,  netabins);
559  
560 <  // response & resolution vs eta/pT
561 <  for (unsigned int ialg=0;ialg<jetalgs.size();ialg++) {
560 >    saveHistos(hPhi,        nphibins);
561 >    saveHistos(hEtCorrPhi,  nphibins);
562 >    saveHistos(hAbsRspPhi,  nphibins);
563 >    saveHistos(hRelRspPhi,  nphibins);
564      
565 <    TDirectory* d = plotfile->GetDirectory(jetalgs[ialg].c_str());
566 <    d->cd();
565 >    saveHistos(hEmf,        nemfbins);
566 >    saveHistos(hEtCorrEmf,  nemfbins);
567 >    saveHistos(hAbsRspEmf,  nemfbins);
568 >    saveHistos(hRelRspEmf,  nemfbins);
569      
570 <    TCanvas* cRspVsEtaPt = plotResponse(jetalgs[ialg],applyjes,fitcorr,
571 <                                        etabins,gRspVsEtaPt[ialg],
572 <                                        "Eta","#eta","Pt","jet p_{T} [GeV]");
573 <    TCanvas* cSgmVsEtaPt = plotResolution(jetalgs[ialg],etabins,gSgmVsEtaPt[ialg],
574 <                                          "Eta","#eta","Pt","jet p_{T} [GeV]");
570 >    saveHistos(hEtEtaEt,    netabins,netbins);
571 >    saveHistos(hEtCorrEtaEt,netabins,netbins);
572 >    saveHistos(hAbsRspEtaEt,netabins,netbins);
573 >    saveHistos(hRelRspEtaEt,netabins,netbins);
574 >
575 >    saveHistos(hEtEmfEt,    nemfbins,netbins);
576 >    saveHistos(hEtCorrEmfEt,nemfbins,netbins);
577 >    saveHistos(hAbsRspEmfEt,nemfbins,netbins);
578 >    saveHistos(hRelRspEmfEt,nemfbins,netbins);
579 >
580 >    saveHistos(hEtEmfEtaEt,    nemfbins,netabins,netbins);
581 >    saveHistos(hEtCorrEmfEtaEt,nemfbins,netabins,netbins);
582 >    saveHistos(hAbsRspEmfEtaEt,nemfbins,netabins,netbins);
583 >    saveHistos(hRelRspEmfEtaEt,nemfbins,netabins,netbins);
584  
585 <    cRspVsEtaPt->cd();
458 <    cSgmVsEtaPt->cd();
585 >    cout<<jetalgs[i]<<": histograms saved."<<endl;
586    }
587    
588    plotfile->Write();
589    plotfile->Close();
590    delete plotfile;
591 <  
592 <  app->Run();
591 >
592 >  //app->Run();
593    
594    return 0;
595   }
# Line 474 | Line 601 | int main(int argc,char**argv)
601   ////////////////////////////////////////////////////////////////////////////////
602  
603   //______________________________________________________________________________
604 < void fillRspAndSgm(TGraphErrors* gRsp,TGraphErrors* gSgm,float nsigma,
605 <                   const vector<float>& bins,TH1F** hVar,TH1F** hPt,TH1F** hRsp)
604 > void replaceHistos(int nbins,
605 >                   TH1F**& hVar,TH1F**& hEtCorr,TH1F**& hAbsRsp,TH1F**& hRelRsp,
606 >                   TTree**& tVar)
607   {
608 <  int nbins = bins.size()+1;
608 >  if (nbins==0) return;
609 >  
610 >  TH1::SetDefaultSumw2();
611    for (int ibin=0;ibin<nbins;ibin++) {
612      
613 <    float var = hVar[ibin]->GetMean();
614 <    float pt  = hPt[ibin]->GetMean();
615 <    TH1F* h   = hRsp[ibin];
616 <    
617 <    float de       = h->GetMean();
618 <    float errde    = h->GetMeanError();
619 <    float sgmde    = h->GetRMS();
620 <    float errsgmde = h->GetRMSError();
621 <    
622 <    if (h->GetEntries()>50) {
623 <      string fname = "fitfnc" + (string)h->GetName();
624 <      TF1* f = new TF1(fname.c_str(),"gaus",de-nsigma*sgmde,de+nsigma*sgmde);
625 <      f->SetLineColor(kRed);
626 <      f->SetLineWidth(1);
627 <      f->SetParameter(1,de);
628 <      f->SetParameter(2,sgmde);
629 <      h->Fit(f,"QR");
630 <      
631 <      de       = f->GetParameter(1);
632 <      errde    = f->GetParError(1);
633 <      sgmde    = f->GetParameter(2);
634 <      errsgmde = f->GetParError(2);
635 <    }
636 <    
637 <    float rsp    = pt/(pt+de);
638 <    float errrsp = pt/(pt+de)/(pt+de)*errde;
639 <    float sgm    = pt/(pt+de)/(pt+de)*sgmde;
640 <    float errsgm = 2*pt*sgmde/(pt+de)/(pt+de)/(pt+de)*errde;
641 <    
642 <    sgm    /= rsp;
643 <    errsgm /= rsp;
644 <    
645 <    gRsp->SetPoint(ibin,var,rsp);
646 <    gRsp->SetPointError(ibin,0.0,errrsp);
647 <    gSgm->SetPoint(ibin,var,sgm);
648 <    gSgm->SetPointError(ibin,0.0,errsgm);
613 >    TH1F* h(0);
614 >    
615 >    if (tVar[ibin]->GetEntries()==0) continue;
616 >
617 >    tVar[ibin]->Project("h(50)","val","weight*(1)");
618 >    h = (TH1F*)gROOT->FindObject("h");
619 >    h->SetNameTitle(hVar[ibin]->GetName(),hVar[ibin]->GetTitle());
620 >    h->SetXTitle(hVar[ibin]->GetXaxis()->GetTitle());
621 >    h->SetYTitle(hVar[ibin]->GetYaxis()->GetTitle());
622 >    delete hVar[ibin];
623 >    hVar[ibin] = h;
624 >
625 >    tVar[ibin]->Project("h(50)","etcorr","weight*(1)");
626 >    h = (TH1F*)gROOT->FindObject("h");
627 >    h->SetNameTitle(hEtCorr[ibin]->GetName(),hEtCorr[ibin]->GetTitle());
628 >    h->SetXTitle(hEtCorr[ibin]->GetXaxis()->GetTitle());
629 >    h->SetYTitle(hEtCorr[ibin]->GetYaxis()->GetTitle());
630 >    delete hEtCorr[ibin];
631 >    hEtCorr[ibin] = h;
632 >    
633 >    tVar[ibin]->Project("h(50)","absrsp","weight*(1)");
634 >    h = (TH1F*)gROOT->FindObject("h");
635 >    h->SetNameTitle(hAbsRsp[ibin]->GetName(),hAbsRsp[ibin]->GetTitle());
636 >    h->SetXTitle(hAbsRsp[ibin]->GetXaxis()->GetTitle());
637 >    h->SetYTitle(hAbsRsp[ibin]->GetYaxis()->GetTitle());
638 >    delete hAbsRsp[ibin];
639 >    hAbsRsp[ibin] = h;
640 >    
641 >    tVar[ibin]->Project("h(50)","relrsp","weight*(1)");
642 >    h = (TH1F*)gROOT->FindObject("h");
643 >    h->SetNameTitle(hRelRsp[ibin]->GetName(),hRelRsp[ibin]->GetTitle());
644 >    h->SetXTitle(hRelRsp[ibin]->GetXaxis()->GetTitle());
645 >    h->SetYTitle(hRelRsp[ibin]->GetYaxis()->GetTitle());
646 >    delete hRelRsp[ibin];
647 >    hRelRsp[ibin] = h;
648 >    
649 >    delete tVar[ibin];
650    }
651    
652 <  gRsp->Write();
522 <  gSgm->Write();
652 >  delete [] tVar;
653   }
654  
655  
656   //______________________________________________________________________________
657 < TCanvas* plotVarHistos(const vector<float>& bins,TH1F** hVar,bool logy)
657 > void replaceHistos(int nbins1,int nbins2,
658 >                   TH1F***& hVar,TH1F***& hEtCorr,TH1F***& hAbsRsp,TH1F***& hRelRsp,
659 >                   TTree***& tVar)
660   {
661 <  string::size_type pos;
530 <  string cname,jetalg;
531 <  
532 <  cname  = hVar[0]->GetName();
533 <  pos    = cname.find_last_of('_');
534 <  jetalg = cname.substr(pos+1);
535 <  cname  = cname.substr(0,pos);
536 <  pos    = cname.find_last_of(':');
537 <  cname  = cname.substr(0,pos);
538 <  cname += "_"+jetalg;
539 <  
540 <  TCanvas* c = new TCanvas(cname.c_str(),cname.c_str());
541 <  
542 <  int nbins = bins.size()+1;
543 <  int ncx= int(std::sqrt((float)nbins));
544 <  int ncy= ncx;
545 <  if (ncx*ncy<nbins) ncx++;
546 <  if (ncx*ncy<nbins) ncy++;
547 <  
548 <  c->Divide(ncx,ncy);
549 <  
550 <  for (int ibin=0;ibin<nbins;ibin++) {
551 <    c->cd(ibin+1);
552 <    gPad->SetLogy(logy);
553 <    hVar[ibin]->Draw("EHIST");
554 <    hVar[ibin]->Write();
555 <  }
556 <  
557 <  //c->Write();
661 >  if (nbins1==0) return;
662    
663 <  return c;
663 >  for (int i1=0;i1<nbins1;i1++)
664 >    replaceHistos(nbins2,hVar[i1],hEtCorr[i1],hAbsRsp[i1],hRelRsp[i1],tVar[i1]);
665 >  delete [] tVar;
666   }
667  
668  
669   //______________________________________________________________________________
670 < TCanvas* plotRspHistos(const vector<float>& bins,TH1F** hRsp)
670 > void replaceHistos(int nbins1,int nbins2,int nbins3,
671 >                   TH1F****& hVar,TH1F****& hEtCorr,TH1F****&hAbsRsp,TH1F****&hRelRsp,
672 >                   TTree****& tVar)
673   {
674 <  string::size_type pos;
675 <  string cname,jetalg;
676 <  
677 <
678 <  cname  = hRsp[0]->GetName();
571 <  pos    = cname.find_last_of('_');
572 <  jetalg = cname.substr(pos+1);
573 <  cname  = cname.substr(0,pos);
574 <  pos    = cname.find_last_of('_');
575 <  cname  = cname.substr(0,pos);
576 <  cname += "_"+jetalg;
577 <  
578 <  TCanvas* c = new TCanvas(cname.c_str(),cname.c_str());
579 <
580 <  int nbins = bins.size()+1;
581 <  int ncx= int(std::sqrt((float)nbins));
582 <  int ncy= ncx;
583 <  if (ncx*ncy<nbins) ncx++;
584 <  if (ncx*ncy<nbins) ncy++;
585 <  
586 <  c->Divide(ncx,ncy);
587 <  
588 <  for (int ibin=0;ibin<nbins;ibin++) {
589 <    c->cd(ibin+1);
590 <    if (hRsp[ibin]->GetEntries()>0) gPad->SetLogy();
591 <    hRsp[ibin]->Draw("EHIST");
592 <    string fname = "fitfnc" + (string)hRsp[ibin]->GetName();
593 <    TF1* fitfnc=hRsp[ibin]->GetFunction(fname.c_str());
594 <    if (0!=fitfnc) fitfnc->Draw("SAME");
595 <    hRsp[ibin]->Write();
596 <  }
597 <  
598 <  c->Write();
599 <  
600 <  return c;
674 >  if (nbins1==0) return;
675 >  for (int i1=0;i1<nbins1;i1++)
676 >    replaceHistos(nbins2,nbins3,
677 >                  hVar[i1],hEtCorr[i1],hAbsRsp[i1],hRelRsp[i1],tVar[i1]);
678 >  delete [] tVar;
679   }
680  
681  
682   //______________________________________________________________________________
683 < TCanvas* plotResponse(const vector<string>& jetalgs,short applyjes,
606 <                      const vector<TGraphErrors*> gRsp,
607 <                      const string& varname,const string& xtitle)
683 > void fitHistos(TH1F** histos,int n,double nsigma)
684   {
685 <  assert(jetalgs.size()==gRsp.size());
686 <  
687 <  string   cname = "RspVs"+varname;
688 <  TCanvas* c     = new TCanvas(cname.c_str(),cname.c_str());
689 <  TLegend* l     = new TLegend(0.7,0.25+jetalgs.size()*0.045,0.85,0.25);
690 <  
691 <  l->SetFillColor(10);
692 <  Color_t color=kBlack;
693 <  TMultiGraph* mg = new TMultiGraph();
694 <  
695 <  for (unsigned int i=0;i<gRsp.size();i++) {
620 <    TGraphErrors* g = gRsp[i];
621 <    l->AddEntry(g,jetalgs[i].c_str(),"l");
622 <    while (color==kYellow||color==kWhite||color==10) color++;
623 <    g->SetMarkerColor(color);
624 <    g->SetLineColor(color);
625 <    g->SetMarkerStyle(kFullCircle);
626 <    g->SetMarkerSize(0.5);
627 <    mg->Add(g,"P");
628 <
629 <    if (varname=="Pt") {
630 <      string fname = "fitfnc"+cname;
631 <      string fitfncAsString =
632 <        //(applyjes==0) ? "[0]+[1]*log(x)+[2]*log(x)*log(x)+[3]*x+[4]*x*x" : "pol0";
633 <        (applyjes==0) ? "[0]*(pow(log(x),[1])-1/x)" : "pol0";
634 <      
635 <      TF1* fitfnc = new TF1(fname.c_str(),fitfncAsString.c_str(),1,get_xmax(g));
636 <      
637 <      if (applyjes==0) {
638 <        fitfnc->SetParameter(0,0.4);
639 <        fitfnc->SetParameter(1,0.4);
640 <      }
641 <      else {
642 <        fitfnc->SetParameter(0,1.0);
643 <      }
644 <      
645 <      fitfnc->SetLineWidth(1);
646 <      fitfnc->SetLineColor(color);
647 <      g->Fit(fitfnc,"QR");
648 <    }
649 <    
650 <    g->Write();
651 <    
652 <    color++;
685 >  for (int i=0;i<n;i++) {
686 >    float nrm  = histos[i]->Integral();
687 >    float mean = histos[i]->GetMean();
688 >    float rms  = histos[i]->GetRMS();
689 >    TF1* f = new TF1("gaus","gaus",mean-nsigma*rms,mean+nsigma*rms);
690 >    f->SetParameter(0,nrm);
691 >    f->SetParameter(1,mean);
692 >    f->SetParameter(2,rms);
693 >    f->SetLineWidth(1);
694 >    f->SetLineColor(kRed);
695 >    histos[i]->Fit(f,"QR"); // Avoid TCanvas?
696    }
654  
655  string mgtitle = "Jet Response Vs "+varname;
656  mg->Draw("a");
657  mg->SetTitle(mgtitle.c_str());
658  mg->GetXaxis()->SetTitle(xtitle.c_str());
659  mg->GetYaxis()->SetTitle("E_{T}^{rec}/E_{T}^{gen}");
660  if (varname=="Pt") mg->GetXaxis()->SetLimits(0.0,mg->GetXaxis()->GetXmax());
661  mg->SetMinimum(0.0);
662  mg->SetMaximum(1.1);
663  l->Draw();
664  
665  c->Write();
666
667  return c;
697   }
698  
699  
700   //______________________________________________________________________________
701 < TCanvas* plotResponse(const string& jetalg,short applyjes,bool fitcorr,
673 <                      const vector<float>& bins1,TGraphErrors** gRsp,
674 <                      const string& varname1,const string vartitle1,
675 <                      const string& varname2,const string& xtitle)
676 <                      
701 > void fitHistos(TH1F*** histos,int n1,int n2,double nsigma)
702   {
703 <  int      nbins = bins1.size()+1;
704 <  
680 <  string   fitfncAsString = "pol0";
681 <  int      fitfncNbParams = 1;
682 <  
683 <  ofstream fout;
684 <  
685 <  if (applyjes==0) {
686 <    //fitfncAsString = "[0]+[1]*log(x)+[2]*log(x)*log(x)+[3]*x+[4]*x*x";
687 <    //fitfncNbParams = 5;
688 <    fitfncAsString = "[0]*(pow(log(x),[1])-1/x)";
689 <    fitfncNbParams = 2;
690 <  }
691 <  if (fitcorr) {
692 <    fout.open((jetalg+".txt").c_str()); assert(fout.is_open());
693 <    fout<<"#\n# jet correction parameters for "<<jetalg<<" algorithm\n#\n\n"
694 <        <<"# number of eta bins\n"<<nbins
695 <        <<"\n\n# fit function\n"<<fitfncAsString<<"\n"<<fitfncNbParams
696 <        <<endl<<endl;
697 <  }
698 <  
699 <  string   cname = "RspVs"+varname1+varname2+"_"+jetalg;
700 <  TCanvas* c = new TCanvas(cname.c_str(),cname.c_str());
701 <  TLegend* l = new TLegend(0.7,0.25+nbins*0.045,0.85,0.25);
702 <  
703 <  l->SetFillColor(10);
704 <  Color_t color=kBlack;
705 <  TMultiGraph* mg = new TMultiGraph();
706 <  
707 <  for (int ibin=0;ibin<nbins;ibin++) {
708 <    
709 <    TGraphErrors* g = gRsp[ibin];
710 <    l->AddEntry(g,hist::get_bin_title(vartitle1,ibin,bins1).c_str(),"l");
711 <    while (color==kYellow||color==kWhite||color==10) color++;
712 <    g->SetMarkerColor(color);
713 <    g->SetLineColor(color);
714 <    g->SetMarkerStyle(kFullCircle);
715 <    g->SetMarkerSize(0.5);
716 <    mg->Add(g,"P");
717 <    
718 <    if (varname2=="Pt") {
719 <      string fname = "fitfnc"+cname;
720 <      TF1* fitfnc = new TF1(fname.c_str(),fitfncAsString.c_str(),1,get_xmax(g));
721 <
722 <      if (applyjes==0) {
723 <        fitfnc->SetParameter(0,0.4);
724 <        fitfnc->SetParameter(1,0.4);
725 <      }
726 <      else {
727 <        fitfnc->SetParameter(0,1.0);
728 <      }
729 <      
730 <      fitfnc->SetLineWidth(1);
731 <      fitfnc->SetLineColor(color);
732 <      g->Fit(fitfnc,"QR");
733 <  
734 <      if (fitcorr==0) {
735 <        if (ibin==nbins-1)
736 <          fout<<"# eta > "<<bins1.back()<<"\n999.9"<<endl;
737 <        else
738 <          fout<<"# eta < "<<bins1[ibin]<<"\n"<<bins1[ibin]<<endl;
739 <        for (int ipar=0;ipar<fitfnc->GetNpar();ipar++)
740 <          fout<<setw(12)<<fitfnc->GetParameter(ipar)<<" "
741 <              <<setw(12)<<fitfnc->GetParError(ipar)<<endl;
742 <        fout<<endl;
743 <      }
744 <      
745 <    }
746 <    
747 <    g->Write();
748 <    
749 <    color++;
750 <  }
751 <  
752 <  if (fout.is_open()) fout.close();
753 <  
754 <  string mgtitle = jetalg+": Jet Response vs "+varname1+" and "+varname2;
755 <  mg->Draw("a");
756 <  mg->SetTitle(mgtitle.c_str());
757 <  mg->GetXaxis()->SetTitle(xtitle.c_str());
758 <  mg->GetYaxis()->SetTitle("E_{T}^{rec}/E_{T}^{gen}");
759 <  if (varname2=="Pt") mg->GetXaxis()->SetLimits(0.0,mg->GetXaxis()->GetXmax());
760 <  mg->SetMinimum(0.0);
761 <  mg->SetMaximum(1.1);
762 <  l->Draw();
763 <
764 <  c->Write();
765 <  
766 <  return c;
703 >  if (n1==0||n2==0) return;
704 >  for (int i1=0;i1<n1;i1++) fitHistos(histos[i1],n2,nsigma);
705   }
706  
707  
708   //______________________________________________________________________________
709 < TCanvas* plotResolution(const vector<string>& jetalgs,
772 <                        const vector<TGraphErrors*> gSgm,
773 <                        const string& varname,const string& xtitle)
709 > void fitHistos(TH1F**** histos,int n1,int n2,int n3,double nsigma)
710   {
711 <  string   cname = "SgmVs"+varname;
712 <  TCanvas* c     = new TCanvas(cname.c_str(),cname.c_str());
713 <  TLegend* l     = new TLegend(0.7,0.85-jetalgs.size()*0.045,0.85,0.85);
778 <  
779 <  l->SetFillColor(10);
780 <  Color_t color=kBlack;
781 <  TMultiGraph* mg = new TMultiGraph();
782 <  
783 <  for (unsigned int i=0;i<gSgm.size();i++) {
784 <    TGraphErrors* g = gSgm[i];
785 <    l->AddEntry(g,jetalgs[i].c_str(),"l");
786 <    while (color==kYellow||color==kWhite||color==10) color++;
787 <    g->SetMarkerColor(color);
788 <    g->SetLineColor(color);
789 <    g->SetMarkerStyle(kFullCircle);
790 <    g->SetMarkerSize(0.5);
791 <    mg->Add(g,"P");
792 <
793 <    if (varname=="Pt") {
794 <      string fname = "fitfnc"+cname;
795 <      TF1* fitfnc = new TF1(fname.c_str(),"sqrt([0]*[0]/x*x+[1]*[1]/x+[2]*[2])");
796 <      fitfnc->SetParameter(0,0.0);
797 <      fitfnc->SetParameter(1,0.5);
798 <      fitfnc->SetParameter(2,0.1);
799 <      fitfnc->SetLineWidth(1);
800 <      fitfnc->SetLineColor(color);
801 <      fitfnc->SetRange(2.0,get_xmax(g));
802 <      g->Fit(fitfnc,"QR");
803 <    }
711 >  if (n1==0||n2==0||n3==0) return;
712 >  for (int i1=0;i1<n1;i1++) fitHistos(histos[i1],n2,n3,nsigma);
713 > }
714  
715 <    g->Write();
716 <    
717 <    color++;
715 >
716 > //______________________________________________________________________________
717 > void saveHistos(TH1F**& histos,int n)
718 > {
719 >  if (n==0) return;
720 >  for (int i=0;i<n;i++) {
721 >    histos[i]->Write();
722 >    delete histos[i];
723    }
724 <  
810 <  string mgtitle = "Jet Resolution Vs "+varname;
811 <  mg->Draw("a");
812 <  mg->SetTitle(mgtitle.c_str());
813 <  mg->GetXaxis()->SetTitle(xtitle.c_str());
814 <  mg->GetYaxis()
815 <    ->SetTitle("#sigma(E_{T}^{rec}/E_{T}^{gen})/<E_{T}^{rec}/E_{T}^{gen}>");
816 <  if (varname=="Pt") mg->GetXaxis()->SetLimits(0.0,mg->GetXaxis()->GetXmax());
817 <  mg->SetMinimum(0.0);
818 <  mg->SetMaximum(0.5);
819 <  l->Draw();  
820 <  
821 <  c->Write();
822 <  
823 <  return c;
724 >  delete [] histos;
725   }
726  
727  
728   //______________________________________________________________________________
729 < TCanvas* plotResolution(const string& jetalg,const vector<float>& bins1,
829 <                        TGraphErrors** gSgm,
830 <                        const string& varname1,const string vartitle1,
831 <                        const string& varname2,const string& xtitle)
729 > void saveHistos(TH1F***& histos,int n1,int n2)
730   {
731 <  int      nbins = bins1.size()+1;
732 <  
733 <  string   cname = "SgmVs"+varname1+varname2+"_"+jetalg;
836 <  TCanvas* c     = new TCanvas(cname.c_str(),cname.c_str());
837 <  TLegend* l     = new TLegend(0.7,0.85-nbins*0.045,0.85,0.85);
838 <  
839 <  l->SetFillColor(10);
840 <  Color_t color=kBlack;
841 <  TMultiGraph* mg = new TMultiGraph();
842 <
843 <  for (int ibin=0;ibin<nbins;ibin++) {
844 <    TGraphErrors* g = gSgm[ibin];
845 <    l->AddEntry(g,hist::get_bin_title(vartitle1,ibin,bins1).c_str(),"l");
846 <    while (color==kYellow||color==kWhite||color==10) color++;
847 <    g->SetMarkerColor(color);
848 <    g->SetLineColor(color);
849 <    g->SetMarkerStyle(kFullCircle);
850 <    g->SetMarkerSize(0.5);
851 <    mg->Add(g,"P");
852 <  
853 <    if (varname2=="Pt") {
854 <      string fname = "fitfnc"+cname;
855 <      TF1* fitfnc = new TF1(fname.c_str(),"sqrt([0]*[0]/x*x+[1]*[1]/x+[2]*[2])");
856 <      fitfnc->SetParameter(0,0.0);
857 <      fitfnc->SetParameter(1,0.5);
858 <      fitfnc->SetParameter(2,0.1);
859 <      fitfnc->SetLineWidth(1);
860 <      fitfnc->SetLineColor(color);
861 <      fitfnc->SetRange(2.0,get_xmax(g));
862 <      g->Fit(fitfnc,"QR");
863 <    }
864 <  
865 <    g->Write();
866 <    
867 <    color++;
868 <  }
869 <  
870 <  string mgtitle=jetalg+": Jet Resolution Vs "+varname1+" and "+varname2;
871 <  mg->Draw("a");
872 <  mg->SetTitle(mgtitle.c_str());
873 <  mg->GetXaxis()->SetTitle(xtitle.c_str());
874 <  mg->GetYaxis()
875 <    ->SetTitle("#sigma(E_{T}^{rec}/E_{T}^{gen})/<E_{T}^{rec}/E_{T}^{gen}>");
876 <  if (varname2=="Pt") mg->GetXaxis()->SetLimits(0.0,mg->GetXaxis()->GetXmax());
877 <  mg->SetMinimum(0.0);
878 <  mg->SetMaximum(0.5);
879 <  l->Draw();
880 <  
881 <  c->Write();
882 <  
883 <  return c;
731 >  if (n1==0) return;
732 >  for (int i1=0;i1<n1;i1++) saveHistos(histos[i1],n2);
733 >  delete [] histos;
734   }
735  
736  
737   //______________________________________________________________________________
738 < float get_xmax(TGraph* g)
738 > void saveHistos(TH1F****& histos,int n1,int n2,int n3)
739   {
740 <  float result(-1e200);
741 <  for (int i=0;i<g->GetN();++i) if (g->GetX()[i]>result) result=g->GetX()[i];
742 <  return result;
740 >  if (n1==0) return;
741 >  for (int i1=0;i1<n1;i1++) saveHistos(histos[i1],n2,n3);
742 >  delete [] histos;
743   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines