ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/NonMCBackground/src/plotZPlusFake.cc
Revision: 1.1
Committed: Fri Jul 6 22:27:04 2012 UTC (12 years, 10 months ago) by dkralph
Content type: text/plain
Branch: MAIN
Log Message:
*** empty log message ***

File Contents

# User Rev Content
1 dkralph 1.1 #include <iostream>
2     #include <algorithm>
3     #include <iomanip>
4     // #include <pair.h>
5    
6     #include "RooHistPdf.h"
7     #include "RooGlobalFunc.h"
8     #include "RooPlot.h"
9     #include "RooRealVar.h"
10     #include "RooDataHist.h"
11     #include "RooGenericPdf.h"
12     #include "RooDataSet.h"
13     #include "RooAddPdf.h"
14     #include "RooFormulaVar.h"
15     #include "RooFitResult.h"
16     #include "RooCurve.h"
17     #include "RooBinning.h"
18     #include "RooExponential.h"
19     #include "RooExtendPdf.h"
20    
21     #include "TCanvas.h"
22     #include "TChain.h"
23     #include "TString.h"
24     #include "TStyle.h"
25     #include "TH1D.h"
26    
27     #include "CPlot.h"
28     #include "FOArgs.h"
29    
30     #include "fake_defs.h"
31     #include "Various.h"
32     #include "CommonDefs.h"
33    
34     #ifndef CMSSW_BASE
35     #define CMSSW_BASE "../../"
36     #endif
37    
38     using namespace std;
39     using namespace RooFit;
40    
41     TCanvas *can;
42    
43     void fillHist(CSample *cs, TString channel, TString type, TString var, double val, double wgt, double wgt_lo=0, double wgt_hi=0);
44     void makeHTML(FOFlags &ctrl, TString type);
45     map<TString,map<TString,TH1D*>* > init_hists(FOFlags &ctrl, TString str="");
46     TString f_to_a(double val) { stringstream ss; ss << fixed << setprecision(2) << val; return TString(ss.str().c_str()); }
47    
48     //----------------------------------------------------------------------------------------
49     int main(int argc, char** argv)
50     {
51     double lumi = 5000;
52     // UInt_t minRun=999999999,maxRun=0;
53    
54     FOFlags ctrl;
55     parse_foargs( argc, argv, ctrl );
56     ctrl.dump();
57    
58     can = new TCanvas("can","can");
59    
60     // FR_struct fr = initFRs(ctrl.mufakedir+"/fr.root",ctrl.elefakedir+"/fr.root");
61     FR_struct fr = initFRs(ctrl.mufakedir,ctrl.elefakedir);
62    
63     TString cmsswpath(CMSSW_BASE + TString("/src"));
64     string xspath = (string(cmsswpath)+"/MitPhysics/data/xs.dat");
65     SimpleTable xstab(xspath.c_str());
66    
67     vector<CSample*> samplev;
68    
69     ifstream instream;
70     instream.open(ctrl.config); assert(instream.is_open());
71     TString ntupledir;
72     string line;
73     while(getline(instream,line)) {
74     if(line[0]=='#') continue;
75    
76     stringstream ss(line);
77     if(line[0]=='^') {
78     TString dummy;
79     if(TString(line).Contains("^ntupdir")) ss >> dummy >> ntupledir;
80     // if(TString(line).Contains("^skim")) ss >> dummy >> ntupledir;
81     continue;
82     }
83    
84     if(line[0]=='$') {
85     samplev.push_back(new CSample());
86     stringstream ss(line);
87     string chr;
88     TString sname;
89     Int_t color;
90     ss >> chr >> sname >> color;
91     string legend = line.substr(line.find('@')+1);
92     samplev.back()->name = sname;
93     samplev.back()->legend = legend;
94     samplev.back()->color = color;
95     samplev.back()->hists = init_hists(ctrl,sname);
96     continue;
97     }
98    
99     TString dset;
100     bool isdata;
101     double xsec;
102     ss >> dset >> isdata >> xsec;
103     assert(ctrl.muSele==ctrl.eleSele);
104     TString fname = ntupledir+"/"+ctrl.muSele+"/"+dset+"/merged.root";
105     samplev.back()->isdata = isdata;
106     (samplev.back()->fsv).push_back(new filestuff(dset,fname,dset,isdata));
107     }
108     instream.close();
109    
110    
111     vector<pair<unsigned,unsigned> > runEvtv; // vector to veto duplicate events
112     for(unsigned ics=0; ics<samplev.size(); ics++) {
113     CSample *cs = samplev[ics];
114     cout << cs->name << endl;
115     for(unsigned ifs=0; ifs<(cs->fsv).size(); ifs++) {
116     filestuff *fs = (cs->fsv)[ifs];
117     cout << "\t" << fs->fname_ << endl;
118     for(unsigned ientry=0; ientry<fs->getentries("FOtree"); ientry++) {
119     fs->getentry(ientry,"","FOtree");
120     fs->getentry(ientry,"info","zznt");
121     // if(fs->isdata_) {
122     // if(fs->info->run < minRun) minRun = fs->info->run;
123     // if(fs->info->run > maxRun) maxRun = fs->info->run;
124     // }
125     bool dupl=false;
126     for(unsigned ievt=0; ievt<runEvtv.size(); ievt++) {
127     if( (runEvtv[ievt].first == fs->info->run) && (runEvtv[ievt].second == fs->info->evt) ) {
128     // cout << "found duplicate event!" << setw(12) << fs->info->run << setw(12) << fs->info->evt << endl;
129     dupl = true;
130     }
131     }
132     if(dupl) continue;
133     runEvtv.push_back(pair<unsigned,unsigned> (fs->info->run,fs->info->evt));
134    
135     double wgt=1;
136     // if(!fs->isdata_) {
137     // double xsWgt = lumi*xstab.Get(fs->dataset_)/fs->total_entries_;
138     // double puWgt = weightTrue2012(fs->info->npu);
139     // wgt = xsWgt*puWgt;
140     // }
141    
142     unsigned npass = fs->passingL->size();
143     unsigned nfail = fs->failingL->size();
144     // if(npass > 2) continue; // exclude signal region (???)
145    
146     float best_z_mass = -99999;
147     pair<int,int> best_z_indices;
148     for(unsigned ilep1=0; ilep1<npass; ilep1++) {
149     SimpleLepton lep1 = (*fs->passingL)[ilep1];
150     for(unsigned ilep2=ilep1+1; ilep2<npass; ilep2++) {
151     SimpleLepton lep2 = (*fs->passingL)[ilep2];
152     if(abs(lep1.type) != abs(lep2.type)) continue;
153     if(lep1.charge == lep2.charge) continue;
154    
155     if(lep1.vec.Pt() < 20 && lep2.vec.Pt() < 20) continue;
156     if(lep1.vec.Pt() < 10 || lep2.vec.Pt() < 10) continue;
157    
158     float tmpM = ((lep1.vec)+(lep2.vec)).M();
159     if( fabs(tmpM-Z_MASS) < fabs(best_z_mass-Z_MASS) ) {
160     best_z_mass = tmpM;
161     best_z_indices = pair<int,int> (ilep1,ilep2);
162     }
163     }
164     }
165     assert( best_z_mass != 0 );
166     if( best_z_mass < 80 || best_z_mass > 100 ) continue;
167     EventData evtdata;
168     evtdata.Z1leptons.push_back((*fs->passingL)[best_z_indices.first]);
169     evtdata.Z1leptons.push_back((*fs->passingL)[best_z_indices.second]);
170    
171     //
172     // Look extra leptons
173     //
174     vector<SimpleLepton> ex_leps;
175     vector<TString> types;
176     for(unsigned ipass=0; ipass<npass; ipass++) {
177     SimpleLepton ex_lep = (*fs->passingL)[ipass];
178     if(ipass == best_z_indices.first) continue;
179     if(ipass == best_z_indices.second) continue;
180     ex_leps.push_back(ex_lep);
181     types.push_back("pass");
182     }
183     for(unsigned ifail=0; ifail<nfail; ifail++) {
184     SimpleLepton ex_lep = (*fs->failingL)[ifail];
185     ex_leps.push_back(ex_lep);
186     types.push_back("fail");
187     }
188    
189     for(unsigned iextra=0; iextra<ex_leps.size(); iextra++) {
190     SimpleLepton ex_lep = ex_leps[iextra];
191    
192     TString channel;
193     if(abs(ex_lep.type) == 13) channel = "mu";
194     if(abs(ex_lep.type) == 11) channel = "ele";
195    
196     // if( !(fabs(lep1.ip3dSig) > 3. &&
197     // fabs(lep2.ip3dSig) > 5.) &&
198     // !(fabs(lep2.ip3dSig) > 3. &&
199     // fabs(lep1.ip3dSig) > 5.) ) continue;
200    
201     evtdata.Z2leptons.clear();
202     //????????????????????????????????????????????????????????????????????????????????????????
203     // hack:
204     evtdata.Z2leptons.push_back((*fs->passingL)[best_z_indices.first]);
205     evtdata.Z2leptons.push_back((*fs->passingL)[best_z_indices.second]);
206     //????????????????????????????????????????????????????????????????????????????????????????
207     KinematicsStruct kine;
208     fillKinematics(evtdata,kine);
209    
210     if(types[iextra]=="pass") {
211     fillHist( cs, channel, "obs", "run", fs->info->run , wgt);
212     fillHist( cs, channel, "obs", "mZ1", kine.mZ1 , wgt);
213     fillHist( cs, channel, "obs", "Z1pt", kine.Z1pt , wgt);
214     fillHist( cs, channel, "obs", "ip3ds_l1", evtdata.Z1leptons[0].ip3dSig , wgt);
215     fillHist( cs, channel, "obs", "ip3ds_l2", evtdata.Z1leptons[1].ip3dSig , wgt);
216     fillHist( cs, channel, "obs", "ip3ds_l3", evtdata.Z2leptons[0].ip3dSig , wgt);
217     fillHist( cs, channel, "obs", "pt_l1", evtdata.Z1leptons[0].vec.Pt() , wgt);
218     fillHist( cs, channel, "obs", "pt_l2", evtdata.Z1leptons[1].vec.Pt() , wgt);
219     fillHist( cs, channel, "obs", "pt_l3", ex_lep.vec.Pt() , wgt);
220     fillHist( cs, channel, "obs", "eta_l1", evtdata.Z1leptons[0].vec.Eta(), wgt);
221     fillHist( cs, channel, "obs", "eta_l2", evtdata.Z1leptons[1].vec.Eta(), wgt);
222     fillHist( cs, channel, "obs", "eta_l3", ex_lep.vec.Eta(), wgt);
223     } else if(types[iextra]=="fail") {
224     double fwgt = get_fake_weight("",ex_lep,fr);
225     double fwgt_lo = get_fake_weight("lo",ex_lep,fr);
226     double fwgt_hi = get_fake_weight("hi",ex_lep,fr);
227     double centr = wgt*fwgt;
228     double lo = wgt*fwgt_lo;
229     double hi = wgt*fwgt_hi;
230     fillHist( cs, channel, "pred", "run", fs->info->run, centr, lo, hi);
231     fillHist( cs, channel, "pred", "mZ1", kine.mZ1, centr, lo, hi);
232     fillHist( cs, channel, "pred", "Z1pt", kine.Z1pt, centr, lo, hi);
233     fillHist( cs, channel, "pred", "ip3ds_l1", evtdata.Z1leptons[0].ip3dSig, centr, lo, hi);
234     fillHist( cs, channel, "pred", "ip3ds_l2", evtdata.Z1leptons[1].ip3dSig, centr, lo, hi);
235     fillHist( cs, channel, "pred", "ip3ds_l3", ex_lep.ip3dSig, centr, lo, hi);
236     fillHist( cs, channel, "pred", "pt_l1", evtdata.Z1leptons[0].vec.Pt(), centr, lo, hi);
237     fillHist( cs, channel, "pred", "pt_l2", evtdata.Z1leptons[1].vec.Pt(), centr, lo, hi);
238     fillHist( cs, channel, "pred", "pt_l3", ex_lep.vec.Pt(), centr, lo, hi);
239     fillHist( cs, channel, "pred", "eta_l1", evtdata.Z1leptons[0].vec.Eta(), centr, lo, hi);
240     fillHist( cs, channel, "pred", "eta_l2", evtdata.Z1leptons[1].vec.Eta(), centr, lo, hi);
241     fillHist( cs, channel, "pred", "eta_l3", ex_lep.vec.Eta(), centr, lo, hi);
242     } else assert(0);
243     }
244     }
245     }
246     }
247     // cout << setw(12) << minRun << setw(12) << maxRun << endl;
248    
249     assert(samplev.size()>0);
250     vector<TString> typev;
251     typev.push_back("both");
252     typev.push_back("mu");
253     typev.push_back("ele");
254     for(unsigned itype=0; itype<typev.size(); itype++) {
255     TString type(typev[itype]);
256     TFile runHistFile(ctrl.outdir+"/"+type+"/runs.root","recreate");
257     map<TString,TH1D*>::iterator it_v;
258     for(it_v=(*(samplev[0]->hists)[type+"_obs"]).begin(); it_v!=(*(samplev[0]->hists)[type+"_obs"]).end(); it_v++) {
259     TString var((*it_v).first);
260     CPlot cplot(var,"",(*(samplev[0]->hists)[type+"_obs"])[var]->GetXaxis()->GetTitle(),"events",ctrl.outdir+"/"+type+"/plots");
261     for(unsigned isam=0; isam<samplev.size(); isam++) {
262     CSample *cs = samplev[isam];
263     TH1D *hObs = (*(cs->hists)[type+"_obs"])[var];
264     TH1D *hPred = (*(cs->hists)[type+"_pred"])[var];
265     TH1D *hPred_lo = (*(cs->hists)[type+"_pred_lo"])[var];
266     TH1D *hPred_hi = (*(cs->hists)[type+"_pred_hi"])[var];
267     if(cs->isdata) {
268     cplot.AddHist1D(hObs,"data: "+integral_str(hObs),"E");
269     double ratio = integrateHist(hObs) / integrateHist(hPred);
270     cplot.AddHist1D(hPred,"FR predic: "+integral_str(hPred),"hist",kRed);
271     cplot.AddHist1D(hPred_lo,"stat lo: "+integral_str(hPred_lo),"hist",kRed,kDashed);
272     cplot.AddHist1D(hPred_hi,"stat hi: "+integral_str(hPred_hi),"hist",kRed,kDashed);
273     } else {
274     cplot.AddToStack(hObs,"mc: "+integral_str(hObs),kCyan-6);
275     }
276     if(cs->isdata && var=="run") {
277     assert(hObs->GetBinContent(0)==0 && hObs->GetBinContent(hObs->GetXaxis()->GetNbins()+1)==0);
278     hObs->Write("runs");
279     }
280     if(var.Contains("ip3d")) {
281     cplot.SetLogy();
282     double maxVal = max(hObs->GetMaximum(), hPred->GetMaximum());
283     double isoMin = 0.00005*maxVal;
284     double isoMax = 1.8*maxVal;
285     cplot.SetYRange(isoMin,isoMax);
286     }
287     }
288     cplot.Draw(can,true,"png");
289     }
290     runHistFile.Close();
291     makeHTML(ctrl,type);
292     }
293     }
294     //----------------------------------------------------------------------------------------
295     map<TString,map<TString,TH1D*>* > init_hists(FOFlags &ctrl, TString str)
296     {
297     map<TString,map<TString,TH1D*>* > hists;
298    
299     map<TString,TH1D*> *h_mu_obs = new map<TString,TH1D*>; hists["mu_obs"] = h_mu_obs;
300     map<TString,TH1D*> *h_mu_pred = new map<TString,TH1D*>; hists["mu_pred"] = h_mu_pred;
301     map<TString,TH1D*> *h_mu_pred_lo = new map<TString,TH1D*>; hists["mu_pred_lo"] = h_mu_pred_lo;
302     map<TString,TH1D*> *h_mu_pred_hi = new map<TString,TH1D*>; hists["mu_pred_hi"] = h_mu_pred_hi;
303    
304     map<TString,TH1D*> *h_ele_obs = new map<TString,TH1D*>; hists["ele_obs"] = h_ele_obs;
305     map<TString,TH1D*> *h_ele_pred = new map<TString,TH1D*>; hists["ele_pred"] = h_ele_pred;
306     map<TString,TH1D*> *h_ele_pred_lo = new map<TString,TH1D*>; hists["ele_pred_lo"] = h_ele_pred_lo;
307     map<TString,TH1D*> *h_ele_pred_hi = new map<TString,TH1D*>; hists["ele_pred_hi"] = h_ele_pred_hi;
308    
309     map<TString,TH1D*> *h_both_obs = new map<TString,TH1D*>; hists["both_obs"] = h_both_obs;
310     map<TString,TH1D*> *h_both_pred = new map<TString,TH1D*>; hists["both_pred"] = h_both_pred;
311     map<TString,TH1D*> *h_both_pred_lo = new map<TString,TH1D*>; hists["both_pred_lo"] = h_both_pred_lo;
312     map<TString,TH1D*> *h_both_pred_hi = new map<TString,TH1D*>; hists["both_pred_hi"] = h_both_pred_hi;
313    
314     map<TString,map<TString,TH1D*>* >::iterator it_h;
315     for(it_h=hists.begin(); it_h!=hists.end(); it_h++) {
316     (*((*it_h).second))["run"] = new TH1D(TString("run") +"_"+(*it_h).first+str,";#bf{run};", 50,160800,196535); (*((*it_h).second))["run"]->Sumw2();
317     (*((*it_h).second))["mZ1"] = new TH1D(TString("mZ1") +"_"+(*it_h).first+str,";#bf{mZ1 [GeV]};", 20,60,120); (*((*it_h).second))["mZ1"]->Sumw2();
318     (*((*it_h).second))["Z1pt"] = new TH1D(TString("Z1pt") +"_"+(*it_h).first+str,";#bf{Z1pt [GeV]};", 20,0,200); (*((*it_h).second))["Z1pt"]->Sumw2();
319    
320     (*((*it_h).second))["pt_l1"] = new TH1D(TString("pt_l1") +"_"+(*it_h).first+str,";#bf{l1 p_{T} [GeV]};", 37,0,120); (*((*it_h).second))["pt_l1"]->Sumw2();
321     (*((*it_h).second))["pt_l2"] = new TH1D(TString("pt_l2") +"_"+(*it_h).first+str,";#bf{l2 p_{T} [GeV]};", 37,0,65); (*((*it_h).second))["pt_l2"]->Sumw2();
322     (*((*it_h).second))["pt_l3"] = new TH1D(TString("pt_l3") +"_"+(*it_h).first+str,";#bf{l3 p_{T} [GeV]};", 37,0,50); (*((*it_h).second))["pt_l3"]->Sumw2();
323    
324     (*((*it_h).second))["eta_l1"] = new TH1D(TString("eta_l1")+"_"+(*it_h).first+str,";#bf{l1 #eta};", 25,-2.5,2.5); (*((*it_h).second))["eta_l1"]->Sumw2();
325     (*((*it_h).second))["eta_l2"] = new TH1D(TString("eta_l2")+"_"+(*it_h).first+str,";#bf{l2 #eta};", 25,-2.5,2.5); (*((*it_h).second))["eta_l2"]->Sumw2();
326     (*((*it_h).second))["eta_l3"] = new TH1D(TString("eta_l3")+"_"+(*it_h).first+str,";#bf{l3 #eta};", 25,-2.5,2.5); (*((*it_h).second))["eta_l3"]->Sumw2();
327    
328     (*((*it_h).second))["ip3ds_l1"] = new TH1D(TString("ip3ds_l1") +"_"+(*it_h).first+str,";#bf{ip3ds l1};",30,0,10); (*((*it_h).second))["ip3ds_l1"]->Sumw2();
329     (*((*it_h).second))["ip3ds_l2"] = new TH1D(TString("ip3ds_l2") +"_"+(*it_h).first+str,";#bf{ip3ds l2};",30,0,10); (*((*it_h).second))["ip3ds_l2"]->Sumw2();
330     (*((*it_h).second))["ip3ds_l3"] = new TH1D(TString("ip3ds_l3") +"_"+(*it_h).first+str,";#bf{ip3ds l3};",30,0,8); (*((*it_h).second))["ip3ds_l3"]->Sumw2();
331     }
332    
333     return hists;
334     }
335     //--------------------------------------------------------------------------------------------------
336     void makeHTML(FOFlags &ctrl, TString type)
337     {
338     TString title("Z Plus Fake: "+type);
339     ofstream htmlfile;
340     char htmlfname[100];
341     sprintf(htmlfname,"%s/plots.html",TString(ctrl.outdir+"/"+type).Data());
342     htmlfile.open(htmlfname);
343    
344     htmlfile << "<!DOCTYPE html" << endl;
345     htmlfile << " PUBLIC \"-//W3C//DTD HTML 3.2//EN\">" << endl;
346     htmlfile << "<html>" << endl;
347    
348     htmlfile << "<head><title>"+title+"</title></head>" << endl;
349     htmlfile << "<body bgcolor=\"EEEEEE\">" << endl;
350     htmlfile << "<h3 style=\"text-align:left; color:DD6600;\">"+title+"</h3>" << endl;
351    
352     htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;
353     htmlfile << "<tr>" << endl;
354     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/mZ1.png\"><img src=\"plots/mZ1.png\" alt=\"plots/mZ1.png\" width=\"100%\"></a></td>" << endl;
355     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/Z1pt.png\"><img src=\"plots/Z1pt.png\" alt=\"plots/Z1pt.png\" width=\"100%\"></a></td>" << endl;
356     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/run.png\"><img src=\"plots/run.png\" alt=\"plots/run.png\" width=\"100%\"></a></td>" << endl;
357     htmlfile << "</tr>" << endl;
358     htmlfile << "</table>" << endl;
359    
360     htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;
361     htmlfile << "<tr>" << endl;
362     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/ip3ds_l1.png\"><img src=\"plots/ip3ds_l1.png\" alt=\"plots/ip3ds_l1.png\" width=\"100%\"></a></td>" << endl;
363     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/ip3ds_l2.png\"><img src=\"plots/ip3ds_l2.png\" alt=\"plots/ip3ds_l2.png\" width=\"100%\"></a></td>" << endl;
364     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/ip3ds_l3.png\"><img src=\"plots/ip3ds_l3.png\" alt=\"plots/ip3ds_l3.png\" width=\"100%\"></a></td>" << endl;
365     htmlfile << "</tr>" << endl;
366     htmlfile << "</table>" << endl;
367    
368     htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;
369     htmlfile << "<tr>" << endl;
370     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/pt_l1.png\"><img src=\"plots/pt_l1.png\" alt=\"plots/pt_l1.png\" width=\"100%\"></a></td>" << endl;
371     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/pt_l2.png\"><img src=\"plots/pt_l2.png\" alt=\"plots/pt_l2.png\" width=\"100%\"></a></td>" << endl;
372     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/pt_l3.png\"><img src=\"plots/pt_l3.png\" alt=\"plots/pt_l3.png\" width=\"100%\"></a></td>" << endl;
373     htmlfile << "</tr>" << endl;
374     htmlfile << "</table>" << endl;
375    
376     htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;
377     htmlfile << "<tr>" << endl;
378     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/eta_l1.png\"><img src=\"plots/eta_l1.png\" alt=\"plots/eta_l1.png\" width=\"100%\"></a></td>" << endl;
379     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/eta_l2.png\"><img src=\"plots/eta_l2.png\" alt=\"plots/eta_l2.png\" width=\"100%\"></a></td>" << endl;
380     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/eta_l3.png\"><img src=\"plots/eta_l3.png\" alt=\"plots/eta_l3.png\" width=\"100%\"></a></td>" << endl;
381     htmlfile << "</tr>" << endl;
382     htmlfile << "</table>" << endl;
383    
384     htmlfile << "<hr />" << endl;
385    
386     htmlfile << "<hr />" << endl;
387    
388     htmlfile << "</body>" << endl;
389     htmlfile << "</html>" << endl;
390     htmlfile.close();
391     }
392     //----------------------------------------------------------------------------------------
393     void fillHist(CSample *cs, TString channel, TString type, TString var, double val, double wgt, double wgt_lo, double wgt_hi)
394     {
395     (*(cs->hists)["both_"+type])[var]->Fill( val, wgt);
396     (*(cs->hists)[channel+"_"+type])[var]->Fill( val, wgt);
397     if(type=="pred") {
398     (*(cs->hists)["both_"+type+"_lo"])[var]->Fill( val, wgt_lo);
399     (*(cs->hists)["both_"+type+"_hi"])[var]->Fill( val, wgt_hi);
400     (*(cs->hists)[channel+"_"+type+"_lo"])[var]->Fill( val, wgt_lo);
401     (*(cs->hists)[channel+"_"+type+"_hi"])[var]->Fill( val, wgt_hi);
402     }
403     }