ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/NonMCBackground/src/plot_emu_fakes.cc
Revision: 1.7
Committed: Mon Jun 25 18:08:04 2012 UTC (12 years, 10 months ago) by dkralph
Content type: text/plain
Branch: MAIN
Changes since 1.6: +15 -5 lines
Log Message: 
*** empty log message ***

File Contents

# User Rev Content
1 dkralph 1.1 #include <iostream>
2     #include <algorithm>
3     #include <iomanip>
4    
5     #include "RooHistPdf.h"
6     #include "RooGlobalFunc.h"
7     #include "RooPlot.h"
8     #include "RooRealVar.h"
9     #include "RooDataHist.h"
10     #include "RooGenericPdf.h"
11     #include "RooDataSet.h"
12     #include "RooAddPdf.h"
13     #include "RooFormulaVar.h"
14     #include "RooFitResult.h"
15     #include "RooCurve.h"
16     #include "RooBinning.h"
17     #include "RooExponential.h"
18 dkralph 1.2 #include "RooExtendPdf.h"
19 dkralph 1.1
20     #include "TCanvas.h"
21     #include "TChain.h"
22     #include "TString.h"
23     #include "TStyle.h"
24     #include "TH1D.h"
25    
26     #include "CPlot.h"
27     #include "FOArgs.h"
28    
29     #include "fake_defs.h"
30 dkralph 1.5 #include "Various.h"
31 dkralph 1.1
32     #ifndef CMSSW_BASE
33     #define CMSSW_BASE "../../"
34     #endif
35    
36     using namespace std;
37     using namespace RooFit;
38    
39     TCanvas *can;
40    
41 dkralph 1.5 double get_fake_weight(TString uncert, SimpleLepton fake_lep, FR_struct fr);
42 dkralph 1.1 void makeHTML(FOFlags &ctrl);
43     map<TString,map<TString,TH1D*>* > init_hists(FOFlags &ctrl, TString str="");
44     TVector3 get_corrected_met(double met, double metphi, TLorentzVector lepton);
45     TString entries_str(TH1D *hist);
46     TString integral_str(TH1D *hist);
47 dkralph 1.2 double integrateHist(TH1D *hist) { return hist->Integral(0,hist->GetNbinsX()+1); }
48 dkralph 1.6 void fillHist(CSample *cs, TString type, TString var, double val, double wgt, double wgt_lo, double wgt_hi);
49 dkralph 1.1
50 dkralph 1.2 RooFitResult *fitMt(FOFlags &ctrl, TH1D *hmcZ, TH1D *hmcW, TH1D *hSSmet, double &total_mc_scale);
51 dkralph 1.1
52     //----------------------------------------------------------------------------------------
53     int main(int argc, char** argv)
54     {
55     double lumi = 1600;
56    
57 dkralph 1.6 TH2D *nBoth = new TH2D("nBoth","nBoth",100,-1,4,100,-1,4);
58    
59 dkralph 1.1 FOFlags ctrl;
60     parse_foargs( argc, argv, ctrl );
61 dkralph 1.4 ctrl.dump();
62 dkralph 1.1
63     can = new TCanvas("can","can");
64    
65     // FR_struct fr = initFRs(ctrl.mufakedir+"/fr.root",ctrl.elefakedir+"/fr.root");
66     FR_struct fr = initFRs(ctrl.mufakedir,ctrl.elefakedir);
67    
68     TString cmsswpath(CMSSW_BASE + TString("/src"));
69     string xspath = (string(cmsswpath)+"/MitPhysics/data/xs.dat");
70     SimpleTable xstab(xspath.c_str());
71    
72     vector<CSample*> samplev;
73    
74     ifstream instream;
75     instream.open(ctrl.config); assert(instream.is_open());
76     TString ntupledir;
77     string line;
78     while(getline(instream,line)) {
79     if(line[0]=='#') continue;
80    
81     stringstream ss(line);
82     if(line[0]=='^') {
83     TString dummy;
84     if(TString(line).Contains("^ntupdir")) ss >> dummy >> ntupledir;
85     // if(TString(line).Contains("^skim")) ss >> dummy >> ntupledir;
86     continue;
87     }
88    
89     if(line[0]=='$') {
90     samplev.push_back(new CSample());
91     stringstream ss(line);
92     string chr;
93     TString sname;
94     Int_t color;
95     ss >> chr >> sname >> color;
96     string legend = line.substr(line.find('@')+1);
97     samplev.back()->name = sname;
98     samplev.back()->legend = legend;
99     samplev.back()->color = color;
100     samplev.back()->hists = init_hists(ctrl,sname);
101     continue;
102     }
103    
104     TString dset;
105     bool isdata;
106     double xsec;
107     ss >> dset >> isdata >> xsec;
108 dkralph 1.4 assert(ctrl.muSele==ctrl.eleSele);
109     TString fname = ntupledir+"/"+ctrl.muSele+"/"+dset+"/merged.root";
110 dkralph 1.1 (samplev.back()->fsv).push_back(new filestuff(dset,fname,dset,isdata));
111     }
112     instream.close();
113    
114    
115     for(unsigned ics=0; ics<samplev.size(); ics++) {
116     CSample *cs = samplev[ics];
117    
118     for(unsigned ifs=0; ifs<(cs->fsv).size(); ifs++) {
119    
120     filestuff *fs = (cs->fsv)[ifs];
121     for(unsigned ientry=0; ientry<fs->getentries(); ientry++) {
122     fs->getentry(ientry,"info");
123    
124     unsigned npass = fs->passingL->size();
125     unsigned nfail = fs->failingL->size();
126    
127 dkralph 1.6 // z veto
128     bool zVeto=false;
129     for(unsigned ilep1=0; ilep1<npass; ilep1++) {
130     SimpleLepton lep1 = (*fs->passingL)[ilep1];
131     for(unsigned ilep2=ilep1+1; ilep2<npass; ilep2++) {
132     SimpleLepton lep2 = (*fs->passingL)[ilep2];
133 dkralph 1.7 // if(abs(lep1.type) == abs(lep2.type)) zVeto = true;
134     // if(lep1.charge != lep2.charge) zVeto = true;
135     if(abs(lep1.type) == abs(lep2.type) &&
136     lep1.charge != lep2.charge) zVeto = true;
137 dkralph 1.6 }
138     }
139     if(zVeto) continue;
140    
141     for(unsigned iw=0; iw<npass; iw++) {
142 dkralph 1.1 SimpleLepton w_lep = (*fs->passingL)[iw];
143    
144 dkralph 1.7 double wgt=1;
145     if(!fs->isdata_) {
146     double xsWgt = lumi*xstab.Get(fs->dataset_)/fs->total_entries_;
147     double puWgt = weightTrue2012(fs->info->npu);
148     wgt = xsWgt*puWgt;
149     }
150 dkralph 1.1
151 dkralph 1.6 // make fake ele and fake mu plots separately
152     if( abs(w_lep.type)==11 && ctrl.faketype=="ele" ||
153     abs(w_lep.type)==13 && ctrl.faketype=="mu" )
154     continue;
155 dkralph 1.7
156 dkralph 1.1 // selection on the W lepton
157     // if(w_lep.isoPF04/w_lep.vec.Pt() > 0.025) continue;
158 dkralph 1.6 if(fs->info->met < 25) continue;
159 dkralph 1.2 if(!(w_lep.tightCutsApplied)) continue;
160 dkralph 1.6 if(w_lep.vec.Pt() < 25) continue;
161 dkralph 1.7 // if(!w_lep.isTight) continue; // this bit stores whether the w lepton is trigger matched
162 dkralph 1.1
163     // require an extra OF SS lepton
164     bool has_ssof=false;
165     vector<SimpleLepton> ssof_leps;
166     vector<TString> types;
167 dkralph 1.6 for(unsigned iextra=iw+1; iextra<npass; iextra++) {
168 dkralph 1.1 SimpleLepton ex_lep = (*fs->passingL)[iextra];
169     if(abs(ex_lep.type) == abs(w_lep.type)) continue;
170     if(ex_lep.charge != w_lep.charge) continue;
171     has_ssof = true;
172     ssof_leps.push_back(ex_lep);
173     types.push_back("pass");
174     }
175 dkralph 1.6 for(unsigned iextra=0; iextra<nfail; iextra++) {
176 dkralph 1.1 SimpleLepton ex_lep = (*fs->failingL)[iextra];
177     if(abs(ex_lep.type) == abs(w_lep.type)) continue;
178     if(ex_lep.charge != w_lep.charge) continue;
179     has_ssof = true;
180     ssof_leps.push_back(ex_lep);
181     types.push_back("fail");
182     }
183     if(!has_ssof) continue;
184    
185     // fill hists with W + (fake or pass)
186     double mt = get_mt(w_lep.vec.Phi(),fs->info->metphi,w_lep.vec.Pt(),fs->info->met);
187 dkralph 1.7 (*(cs->hists)["both"])["npv"]->Fill( fs->info->npv, wgt);
188 dkralph 1.1 (*(cs->hists)["both"])["w_pt"]->Fill( w_lep.vec.Pt(), wgt);
189     (*(cs->hists)["both"])["w_eta"]->Fill( w_lep.vec.Eta(), wgt);
190     (*(cs->hists)["both"])["w_isoPF04"]->Fill( w_lep.isoPF04/w_lep.vec.Pt(), wgt);
191     (*(cs->hists)["both"])["met"]->Fill( fs->info->met, wgt);
192     (*(cs->hists)["both"])["mt"]->Fill( mt, wgt);
193 dkralph 1.6 (*(cs->hists)["both"])["npass"]->Fill( npass, wgt);
194     (*(cs->hists)["both"])["nfail"]->Fill( nfail, wgt);
195     nBoth->Fill(npass,nfail,wgt);
196 dkralph 1.1
197     for(unsigned iextra=0; iextra<ssof_leps.size(); iextra++) {
198     SimpleLepton ex_lep = ssof_leps[iextra];
199     double mass = (w_lep.vec + ex_lep.vec).M();
200 dkralph 1.2
201 dkralph 1.1 if(types[iextra] == "fail") {
202    
203     // fill fake rate prediction
204 dkralph 1.5
205     double fake_wgt = get_fake_weight("",ex_lep,fr);
206     double fake_wgt_lo = get_fake_weight("lo",ex_lep,fr);
207     double fake_wgt_hi = get_fake_weight("hi",ex_lep,fr);
208    
209 dkralph 1.6 double centr=fake_wgt*wgt,lo=fake_wgt_lo*wgt,hi=fake_wgt_hi*wgt;
210     fillHist( cs, "pred", "w_pt", w_lep.vec.Pt(), centr, lo, hi);
211     fillHist( cs, "pred", "l2_pt", ex_lep.vec.Pt(), centr, lo, hi);
212     fillHist( cs, "pred", "w_eta", w_lep.vec.Eta(), centr, lo, hi);
213     fillHist( cs, "pred", "l2_eta", ex_lep.vec.Eta(), centr, lo, hi);
214     fillHist( cs, "pred", "w_isoPF04", w_lep.isoPF04/w_lep.vec.Pt(), centr, lo, hi);
215     fillHist( cs, "pred", "l2_isoPF04", ex_lep.isoPF04/ex_lep.vec.Pt(), centr, lo, hi);
216     fillHist( cs, "pred", "mass", mass, centr, lo, hi);
217     fillHist( cs, "pred", "met", fs->info->met, centr, lo, hi);
218     fillHist( cs, "pred", "mt", mt, centr, lo, hi);
219    
220 dkralph 1.3 (*(cs->hists)["pred"])["fr"]->Fill( fake_wgt );
221 dkralph 1.6
222 dkralph 1.1 } else if(types[iextra] == "pass") {
223    
224     // fill the observation
225    
226 dkralph 1.6 fillHist( cs, "obs", "w_pt", w_lep.vec.Pt(), wgt, 0, 0);
227     fillHist( cs, "obs", "l2_pt", ex_lep.vec.Pt(), wgt, 0, 0);
228     fillHist( cs, "obs", "w_eta", w_lep.vec.Eta(), wgt, 0, 0);
229     fillHist( cs, "obs", "l2_eta", ex_lep.vec.Eta(), wgt, 0, 0);
230     fillHist( cs, "obs", "w_isoPF04", w_lep.isoPF04/w_lep.vec.Pt(), wgt, 0, 0);
231     fillHist( cs, "obs", "l2_isoPF04", ex_lep.isoPF04/ex_lep.vec.Pt(), wgt, 0, 0);
232     fillHist( cs, "obs", "mass", mass, wgt, 0, 0);
233     fillHist( cs, "obs", "met", fs->info->met, wgt, 0, 0);
234     fillHist( cs, "obs", "mt", mt, wgt, 0, 0);
235 dkralph 1.1 } else { cout << "ERROR: types: " << types[iextra] << endl; assert(0); }
236     }
237     }
238     }
239     }
240     }
241    
242     assert(samplev.size()==3);
243     CSample *cs = samplev[0];
244     CSample *cs_mc = samplev[1];
245     CSample *cs_zmc = samplev[2];
246    
247     double total_mc_scale=-1;
248 dkralph 1.2 RooFitResult *fitres = fitMt(ctrl,(*(cs_zmc->hists)["both"])["mt"],(*(cs_mc->hists)["both"])["mt"],(*(cs->hists)["both"])["mt"],total_mc_scale);
249 dkralph 1.1
250 dkralph 1.6 cout << "\n\n\ntotal mc scale: " << total_mc_scale << endl << endl << endl;
251    
252 dkralph 1.1 // RooArgList float_pars = fitres->floatParsFinal();
253     // RooRealVar *nWandZ = (RooRealVar*)float_pars.at(float_pars.index("nWandZ"));
254     // cout << nWandZ->getVal() << endl;
255    
256     map<TString,TH1D*>::iterator it_v;
257     for(it_v=(*(cs->hists)["pred"]).begin(); it_v!=(*(cs->hists)["pred"]).end(); it_v++) {
258     TString var((*it_v).first);
259    
260 dkralph 1.6 // plot passing and failing fakes together
261 dkralph 1.2 if(var!="mt") {
262 dkralph 1.1 (*(cs_mc->hists)["both"])[var]->Scale(total_mc_scale);
263     (*(cs_zmc->hists)["both"])[var]->Scale(total_mc_scale);
264     }
265     CPlot cplot_both("both_"+var,"",(*(cs->hists)["both"])[var]->GetXaxis()->GetTitle(),"events",ctrl.outdir+"/plots");
266     cplot_both.AddHist1D((*(cs->hists)["both"])[var],"tight+fk'ble: "+integral_str((*(cs->hists)["both"])[var]),"E");
267     cplot_both.AddToStack((*(cs_zmc->hists)["both"])[var],"Z+j MC: "+integral_str((*(cs_zmc->hists)["both"])[var]),920);
268     cplot_both.AddToStack((*(cs_mc->hists)["both"])[var],"W+j MC: "+integral_str((*(cs_mc->hists)["both"])[var]),426);
269    
270     if(var.Contains("eta")) cplot_both.SetYRange(0,1.4*(*(cs->hists)["both"])[var]->GetMaximum());
271     if(var.Contains("isoPF")) cplot_both.SetLogy();
272 dkralph 1.2 if(var=="met") cplot_both.AddLine(25,0,25,18050,843,kDashed);
273 dkralph 1.1 cplot_both.Draw(can,true,"png");
274    
275     // plot observation vs. prediction
276     (*(cs_mc->hists)["obs"])[var]->Scale(total_mc_scale);
277     (*(cs_zmc->hists)["obs"])[var]->Scale(total_mc_scale);
278     CPlot cplot(var,"",(*(cs->hists)["pred"])[var]->GetXaxis()->GetTitle(),"events",ctrl.outdir+"/plots");
279     cplot.AddHist1D((*(cs->hists)["obs"])[var],"Obs: "+integral_str((*(cs->hists)["obs"])[var]),"E");
280     cplot.AddHist1D((*(cs->hists)["pred"])[var],"FR predic: "+integral_str((*(cs->hists)["pred"])[var]),"hist",kRed);
281 dkralph 1.6 cplot.AddHist1D((*(cs->hists)["pred_lo"])[var],"stat lo: "+integral_str((*(cs->hists)["pred_lo"])[var]),"hist",kRed,kDashed);
282     cplot.AddHist1D((*(cs->hists)["pred_hi"])[var],"stat hi: "+integral_str((*(cs->hists)["pred_hi"])[var]),"hist",kRed,kDashed);
283    
284     // stack on the mc
285 dkralph 1.5 // cplot.AddToStack((*(cs_zmc->hists)["obs"])[var],"Z+j MC: "+integral_str((*(cs_zmc->hists)["obs"])[var]),920);
286     // cplot.AddToStack((*(cs_mc->hists)["obs"])[var],"W+j MC: "+integral_str((*(cs_mc->hists)["obs"])[var]),426);
287 dkralph 1.1
288     if(var.Contains("eta")) cplot.SetYRange(0,1.4*max((*(cs->hists)["obs"])[var]->GetMaximum(),(*(cs->hists)["pred"])[var]->GetMaximum()));
289 dkralph 1.6 if(var.Contains("isoPF")) {
290     cplot.SetLogy();
291     TH1D *hobs = (*(cs->hists)["obs"])[var];
292     TH1D *hpred = (*(cs->hists)["pred"])[var];
293     double isoMin = 0.008*min(hobs->GetBinContent(1), hpred->GetBinContent(1));
294     double isoMax = 1.8*max(hobs->GetBinContent(1), hpred->GetBinContent(1));
295    
296     cplot.SetYRange(isoMin,isoMax);
297     // cplot.SetYRange(0,0.8*min((*(cs->hists)["obs"])[var]->GetMinimum(),(*(cs->hists)["pred"])[var]->GetMinimum()));
298     }
299 dkralph 1.1 cplot.Draw(can,true,"png");
300     }
301    
302     makeHTML(ctrl);
303 dkralph 1.6 nBoth->Draw("surf");
304     can->SaveAs("~/public_html/foo.png");
305     can->SaveAs("~/public_html/foo.C");
306     nBoth->Draw("colz");
307     can->SaveAs("~/public_html/fooc.png");
308     can->SaveAs("~/public_html/fooc.C");
309 dkralph 1.1 }
310     //----------------------------------------------------------------------------------------
311 dkralph 1.5 double get_fake_weight(TString uncert, SimpleLepton fake_lep, FR_struct fr)
312 dkralph 1.1 {
313     double rate=0;
314 dkralph 1.5 assert(uncert=="" || uncert=="lo" || uncert=="hi");
315 dkralph 1.1 if(abs(fake_lep.type) == 13) { // fake muons:
316 dkralph 1.5 double x = fake_lep.vec.Eta();
317     if(fr.absEtaMu) x = fabs(x);
318     double y = fake_lep.vec.Pt();
319     rate = fr.mufr.getEff(x,y,true);
320     if(uncert == "lo") rate -= fr.mufr.getErrLow(x,y,true);
321     if(uncert == "hi") rate += fr.mufr.getErrHigh(x,y,true);
322 dkralph 1.1 } else if(abs(fake_lep.type) == 11) {
323 dkralph 1.5 double x = fake_lep.vec.Eta();
324     if(fr.absEtaEle) x = fabs(x);
325     double y = fake_lep.vec.Pt();
326     rate = fr.elefr.getEff(x,y,true);
327     if(uncert == "lo") rate -= fr.elefr.getErrLow(x,y,true);
328     if(uncert == "hi") rate += fr.elefr.getErrHigh(x,y,true);
329 dkralph 1.1 } else assert(0);
330    
331     return rate/(1-rate);
332     }
333     //----------------------------------------------------------------------------------------
334     map<TString,map<TString,TH1D*>* > init_hists(FOFlags &ctrl, TString str)
335     {
336 dkralph 1.5 map<TString,map<TString,TH1D*>* > hists;
337     map<TString,TH1D*> *h_both = new map<TString,TH1D*>; hists["both"] = h_both;
338     map<TString,TH1D*> *h_pred = new map<TString,TH1D*>; hists["pred"] = h_pred;
339     map<TString,TH1D*> *h_pred_lo = new map<TString,TH1D*>; hists["pred_lo"] = h_pred_lo;
340     map<TString,TH1D*> *h_pred_hi = new map<TString,TH1D*>; hists["pred_hi"] = h_pred_hi;
341     map<TString,TH1D*> *h_obs = new map<TString,TH1D*>; hists["obs"] = h_obs;
342 dkralph 1.1 map<TString,map<TString,TH1D*>* >::iterator it_h;
343     double w_pt_max,l2_pt_max,mass_max;
344     if(ctrl.faketype=="mu") {
345     w_pt_max = 80;
346 dkralph 1.5 l2_pt_max = 35;
347 dkralph 1.1 mass_max = 120;
348     } else if(ctrl.faketype=="ele"){
349     w_pt_max = 90;
350     l2_pt_max = 60;
351     mass_max = 210;
352     } else assert(0);
353    
354     for(it_h=hists.begin(); it_h!=hists.end(); it_h++) {
355 dkralph 1.7 (*((*it_h).second))["npv"] = new TH1D(TString("npv") +"_"+(*it_h).first+str,";#bf{n PV};", 50,-.5,49.5); (*((*it_h).second))["npv"]->Sumw2();
356 dkralph 1.3 (*((*it_h).second))["w_pt"] = new TH1D(TString("w_pt") +"_"+(*it_h).first+str,";#bf{real p_{T} [GeV]};", 37,15,w_pt_max); (*((*it_h).second))["w_pt"]->Sumw2();
357     (*((*it_h).second))["l2_pt"] = new TH1D(TString("l2_pt") +"_"+(*it_h).first+str,";#bf{fake p_{T} [GeV]};", 25,0,l2_pt_max); (*((*it_h).second))["l2_pt"]->Sumw2();
358     (*((*it_h).second))["w_eta"] = new TH1D(TString("w_eta")+"_"+(*it_h).first+str,";#bf{real #eta};", 25,-2.5,2.5); (*((*it_h).second))["w_eta"]->Sumw2();
359     (*((*it_h).second))["l2_eta"] = new TH1D(TString("l2_eta")+"_"+(*it_h).first+str,";#bf{fake #eta};", 25,-2.5,2.5); (*((*it_h).second))["l2_eta"]->Sumw2();
360 dkralph 1.1 (*((*it_h).second))["w_isoPF04"] = new TH1D(TString("w_isoPF04")+"_"+(*it_h).first+str,";#bf{real PF Iso};", 25,0,1); (*((*it_h).second))["w_isoPF04"]->Sumw2();
361     (*((*it_h).second))["l2_isoPF04"] = new TH1D(TString("l2_isoPF04")+"_"+(*it_h).first+str,";#bf{fake PF Iso};", 25,0,3); (*((*it_h).second))["l2_isoPF04"]->Sumw2();
362 dkralph 1.3 (*((*it_h).second))["mass"] = new TH1D(TString("mass") +"_"+(*it_h).first+str,";#bf{m_{ll}};", 25,0,mass_max); (*((*it_h).second))["mass"]->Sumw2();
363     (*((*it_h).second))["met"] = new TH1D(TString("met") +"_"+(*it_h).first+str,";#bf{MET};", 25,0,100); (*((*it_h).second))["met"]->Sumw2();
364     (*((*it_h).second))["mt"] = new TH1D(TString("mt") +"_"+(*it_h).first+str,";#bf{m_{T}};", 25,2,150); (*((*it_h).second))["mt"]->Sumw2();
365     (*((*it_h).second))["fr"] = new TH1D(TString("fr") +"_"+(*it_h).first+str,";#bf{FR};", 125,-.01,.1); (*((*it_h).second))["fr"]->Sumw2();
366 dkralph 1.6 (*((*it_h).second))["npass"] = new TH1D(TString("npass") +"_"+(*it_h).first+str,";#bf{N Pass};", 125,-.01,4); (*((*it_h).second))["npass"]->Sumw2();
367     (*((*it_h).second))["nfail"] = new TH1D(TString("nfail") +"_"+(*it_h).first+str,";#bf{N Fail};", 125,-.01,4); (*((*it_h).second))["nfail"]->Sumw2();
368 dkralph 1.1 }
369    
370     return hists;
371     }
372     //--------------------------------------------------------------------------------------------------
373     void makeHTML(FOFlags &ctrl)
374     {
375     // TString title(ctrl.inputdir);
376     // title = title(title.Last('/')+1,title.Length());
377     TString title("EMU closure test: Fake "+ctrl.faketype);
378     ofstream htmlfile;
379     char htmlfname[100];
380     sprintf(htmlfname,"%s/"+ctrl.faketype+".html",ctrl.outdir.Data());
381     htmlfile.open(htmlfname);
382    
383     htmlfile << "<!DOCTYPE html" << endl;
384     htmlfile << " PUBLIC \"-//W3C//DTD HTML 3.2//EN\">" << endl;
385     htmlfile << "<html>" << endl;
386    
387     htmlfile << "<head><title>"+title+"</title></head>" << endl;
388     htmlfile << "<body bgcolor=\"EEEEEE\">" << endl;
389     htmlfile << "<h3 style=\"text-align:left; color:DD6600;\">"+title+"</h3>" << endl;
390    
391     htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;
392     htmlfile << "<tr>" << endl;
393 dkralph 1.2 htmlfile << "<td width=\"33.3%\"><a target=\"_blank\" href=\"plots/mtfitW.png\"><img src=\"plots/mtfitW.png\" alt=\"plots/mtfitW.png\" width=\"100%\"></a></td>" << endl;
394     htmlfile << "<td width=\"33.3%\"><a target=\"_blank\" href=\"plots/both_met.png\"><img src=\"plots/both_met.png\" alt=\"plots/both_met.png\" width=\"100%\"></a></td>" << endl;
395 dkralph 1.1 htmlfile << "</tr>" << endl;
396     htmlfile << "</table>" << endl;
397    
398     htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;
399     htmlfile << "<tr>" << endl;
400 dkralph 1.3 htmlfile << "<td width=\"33.3%\"><a target=\"_blank\" href=\"plots/fr.png\"><img src=\"plots/fr.png\" alt=\"plots/fr.png\" width=\"100%\"></a></td>" << endl;
401 dkralph 1.7 htmlfile << "<td width=\"33.3%\"><a target=\"_blank\" href=\"plots/both_npv.png\"><img src=\"plots/both_npv.png\" alt=\"plots/both_npv.png\" width=\"100%\"></a></td>" << endl;
402 dkralph 1.3 htmlfile << "<td width=\"33.3%\"><a></a></td>" << endl;
403     htmlfile << "</tr>" << endl;
404     htmlfile << "</table>" << endl;
405    
406     htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;
407     htmlfile << "<tr>" << endl;
408 dkralph 1.1 htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/mass.png\"><img src=\"plots/mass.png\" alt=\"plots/mass.png\" width=\"100%\"></a></td>" << endl;
409     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/mt.png\"><img src=\"plots/mt.png\" alt=\"plots/mt.png\" width=\"100%\"></a></td>" << endl;
410     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/met.png\"><img src=\"plots/met.png\" alt=\"plots/met.png\" width=\"100%\"></a></td>" << endl;
411     htmlfile << "</tr>" << endl;
412     htmlfile << "</table>" << endl;
413    
414     htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;
415     htmlfile << "<tr>" << endl;
416     htmlfile << "<td width=\"33.3%\"><a target=\"_blank\" href=\"plots/w_pt.png\"><img src=\"plots/w_pt.png\" alt=\"plots/w_pt.png\" width=\"100%\"></a></td>" << endl;
417     htmlfile << "<td width=\"33.3%\"><a target=\"_blank\" href=\"plots/w_eta.png\"><img src=\"plots/w_eta.png\" alt=\"plots/w_eta.png\" width=\"100%\"></a></td>" << endl;
418     htmlfile << "<td width=\"33.3%\"><a target=\"_blank\" href=\"plots/w_isoPF04.png\"><img src=\"plots/w_isoPF04.png\" alt=\"plots/w_isoPF04.png\" width=\"100%\"></a></td>" << endl;
419     htmlfile << "</tr>" << endl;
420     htmlfile << "</table>" << endl;
421    
422     htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;
423     htmlfile << "<tr>" << endl;
424     htmlfile << "<td width=\"33.3%\"><a target=\"_blank\" href=\"plots/l2_pt.png\"><img src=\"plots/l2_pt.png\" alt=\"plots/l2_pt.png\" width=\"100%\"></a></td>" << endl;
425     htmlfile << "<td width=\"33.3%\"><a target=\"_blank\" href=\"plots/l2_eta.png\"><img src=\"plots/l2_eta.png\" alt=\"plots/l2_eta.png\" width=\"100%\"></a></td>" << endl;
426     htmlfile << "<td width=\"33.3%\"><a target=\"_blank\" href=\"plots/l2_isoPF04.png\"><img src=\"plots/l2_isoPF04.png\" alt=\"plots/l2_isoPF04.png\" width=\"100%\"></a></td>" << endl;
427     htmlfile << "</tr>" << endl;
428     htmlfile << "</table>" << endl;
429     htmlfile << "<hr />" << endl;
430    
431 dkralph 1.6 // htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;
432     // htmlfile << "<tr>" << endl;
433     // htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/both_w_isoPF04.png\"><img src=\"plots/both_w_isoPF04.png\" alt=\"plots/both_w_isoPF04.png\" width=\"100%\"></a></td>" << endl;
434     // htmlfile << "</tr>" << endl;
435     // htmlfile << "</table>" << endl;
436 dkralph 1.1
437     htmlfile << "<hr />" << endl;
438    
439     htmlfile << "<hr />" << endl;
440    
441     htmlfile << "</body>" << endl;
442     htmlfile << "</html>" << endl;
443     htmlfile.close();
444     }
445     //----------------------------------------------------------------------------------------
446     TVector3 get_corrected_met(double met, double metphi, TLorentzVector lepton)
447     {
448     TLorentzVector corr_met;
449     corr_met.SetPtEtaPhiM(met,0,metphi,0);
450     TVector3 w_lep_trans(lepton.Px(),lepton.Py(),0);
451     return corr_met.Vect() + w_lep_trans;
452     }
453     //----------------------------------------------------------------------------------------
454     TString entries_str(TH1D *hist)
455     {
456     int entries = hist->GetEntries();
457     stringstream ss;
458     ss << entries;
459     TString str;
460     ss >> str;
461     return str;
462     }
463     //----------------------------------------------------------------------------------------
464     TString integral_str(TH1D *hist)
465     {
466     stringstream ss;
467 dkralph 1.2 ss << fixed << setprecision(1) << integrateHist(hist);
468 dkralph 1.1 TString str;
469     ss >> str;
470     return str;
471     }
472     //----------------------------------------------------------------------------------------
473 dkralph 1.2 RooFitResult *fitMt(FOFlags &ctrl, TH1D *hmcZ, TH1D *hmcW, TH1D *hdata, double &total_mc_scale)
474 dkralph 1.1 {
475 dkralph 1.2 double original_w_integral = integrateHist(hmcW);
476 dkralph 1.1
477 dkralph 1.2 RooRealVar mt("mt","MT [GeV]",hdata->GetXaxis()->GetXmin(),hdata->GetXaxis()->GetXmax());
478     RooBinning mtbins(hdata->GetNbinsX(),hdata->GetXaxis()->GetXmin(),hdata->GetXaxis()->GetXmax());
479     mt.setBinning(mtbins);
480 dkralph 1.1
481     // data
482 dkralph 1.2 RooDataHist obsMt("obsMt","obsMt",RooArgList(mt),hdata);
483 dkralph 1.1 // W mc pdf
484     hmcW->Smooth(2);
485 dkralph 1.2 RooDataHist hW("hW","hW",RooArgList(mt),hmcW);
486     RooHistPdf pW("pW","pW",RooArgSet(mt),hW,3);
487     // Z mc pdf
488     RooDataHist hZ("hZ","hZ",RooArgList(mt),hmcZ);
489     RooHistPdf pZ("pZ","pZ",RooArgSet(mt),hZ,3);
490 dkralph 1.1 // rayleigh pdf
491 dkralph 1.2 RooRealVar sig("sig","sig",400,15,585);
492     RooRealVar a1("a1","a1",-3,-10,10);
493     // RooGenericPdf rayl("rayl","rayl","mt*exp(-(mt)*(mt)/(2*sig*sig))",RooArgList(mt,sig));
494     RooGenericPdf rayl("rayl","rayl","mt*exp(-(mt)*(mt)/(sig*sig +2*sig*a1*(mt) + a1*a1*(mt)*(mt)))",RooArgList(mt,sig,a1));
495    
496     RooRealVar nObs("nObs","nObs",integrateHist(hdata));
497     RooRealVar nWandZ("nWandZ","nWandZ",0.98*nObs.getVal(),.95*nObs.getVal(),1.1*nObs.getVal());
498     RooRealVar zFrac("zFrac","zFrac",integrateHist(hmcZ) / (integrateHist(hmcZ) + integrateHist(hmcW)));
499     RooAddPdf modelWandZ("modelWandZ","modelWandZ",RooArgList(pZ,pW),zFrac);
500     RooRealVar nRayl("nRayl","nRayl",0.05*nObs.getVal(),0.0*nObs.getVal(),.05*nObs.getVal());
501    
502     RooExtendPdf extWandZ("extWandZ","extWandZ",modelWandZ,nWandZ);
503     RooExtendPdf extRayl("extRayl","extRayl",rayl,nRayl);
504    
505     RooAddPdf *model = new RooAddPdf("model","model",RooArgList(extWandZ,extRayl));
506    
507     RooFitResult *fitres = model->fitTo(obsMt,Extended(kTRUE),Save(kTRUE));
508     fitres->Print();
509    
510     TH1D *rayl_hist = (TH1D*)rayl.createHistogram("rayl_mt",mt,Binning(mtbins));
511     rayl_hist->Scale(nRayl.getVal()/integrateHist(rayl_hist));
512     TH1D *totalHist = (TH1D*)model->createHistogram("total_mt",mt,Binning(mtbins));
513     totalHist->Scale( (nWandZ.getVal() + nRayl.getVal()) / integrateHist(totalHist) );
514    
515     // // why the fuck do I need this?
516     // double tot = nWandZ.getVal() + nRayl.getVal();
517     // if(use_exp) tot += nExp.getVal();
518     // double hack_fac = integrateHist(hdata)/tot;
519     // rayl_hist->Scale(hack_fac);
520     // exp_hist->Scale(hack_fac);
521     // hmcZ->Scale(hack_fac);
522     // hmcW->Scale(hack_fac);
523    
524     // RooPlot *frame = mt.frame();
525     // obsMt.plotOn(frame);
526     // model->plotOn(frame);
527     // model->plotOn(frame,Components("pW"),LineStyle(kDashed),LineColor(kRed));
528     // model->plotOn(frame,Components("pZ"),LineColor(kRed));
529     // model->plotOn(frame,Components("rayl"),LineColor(kGreen));
530     // frame->Draw();
531     // can->SaveAs("~/public_html/foo.png");
532    
533     hmcW->Scale( (nWandZ.getVal()*(1-zFrac.getVal())) / integrateHist(hmcW) );
534     hmcZ->Scale( (nWandZ.getVal()*zFrac.getVal()) / integrateHist(hmcZ) );
535 dkralph 1.1
536 dkralph 1.2 CPlot plotMtFit("mtfitW","","MT [GeV]","Events",ctrl.outdir+"/plots");
537     plotMtFit.AddHist1D(hdata,"tight+fk'ble: " + integral_str(hdata),"E");
538     plotMtFit.AddHist1D(totalHist,"total MC: " + integral_str(totalHist),"hist",kRed);
539     plotMtFit.AddToStack(rayl_hist,"jet (Rayleigh): " + integral_str(rayl_hist),kBlue);
540     plotMtFit.AddToStack(hmcZ,"Z+j MC: " + integral_str(hmcZ),kGray);
541     plotMtFit.AddToStack(hmcW,"W+j MC: " + integral_str(hmcW),kCyan-6);
542     plotMtFit.Draw(can,kTRUE,"png");
543 dkralph 1.1
544 dkralph 1.2 total_mc_scale = integrateHist(hmcW) / original_w_integral;
545 dkralph 1.1
546     return fitres;
547     }
548 dkralph 1.6 //----------------------------------------------------------------------------------------
549     void fillHist(CSample *cs, TString type, TString var, double val, double wgt, double wgt_lo, double wgt_hi)
550     {
551     (*(cs->hists)[type])[var]->Fill( val, wgt);
552     if(type=="pred") {
553     (*(cs->hists)[type+"_lo"])[var]->Fill( val, wgt_lo);
554     (*(cs->hists)[type+"_hi"])[var]->Fill( val, wgt_hi);
555     }
556     }