ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/NonMCBackground/src/plot_emu_fakes.cc
Revision: 1.6
Committed: Thu Jun 21 21:40:30 2012 UTC (12 years, 10 months ago) by dkralph
Content type: text/plain
Branch: MAIN
Changes since 1.5: +88 -59 lines
Log Message: 
Lots of stuff

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     if(abs(lep1.type) == abs(lep2.type)) zVeto = true;
134     if(lep1.charge != lep2.charge) zVeto = true;
135     }
136     }
137     if(zVeto) continue;
138    
139     for(unsigned iw=0; iw<npass; iw++) {
140 dkralph 1.1 SimpleLepton w_lep = (*fs->passingL)[iw];
141    
142     double wgt = fs->isdata_ ? 1 : lumi*xstab.Get(fs->dataset_)/fs->total_entries_;
143    
144 dkralph 1.6 // make fake ele and fake mu plots separately
145     if( abs(w_lep.type)==11 && ctrl.faketype=="ele" ||
146     abs(w_lep.type)==13 && ctrl.faketype=="mu" )
147     continue;
148    
149 dkralph 1.1 // selection on the W lepton
150     // if(w_lep.isoPF04/w_lep.vec.Pt() > 0.025) continue;
151 dkralph 1.6 if(fs->info->met < 25) continue;
152 dkralph 1.2 if(!(w_lep.tightCutsApplied)) continue;
153 dkralph 1.6 if(w_lep.vec.Pt() < 25) continue;
154 dkralph 1.1
155     // require an extra OF SS lepton
156     bool has_ssof=false;
157     vector<SimpleLepton> ssof_leps;
158     vector<TString> types;
159 dkralph 1.6 for(unsigned iextra=iw+1; iextra<npass; iextra++) {
160 dkralph 1.1 SimpleLepton ex_lep = (*fs->passingL)[iextra];
161     if(abs(ex_lep.type) == abs(w_lep.type)) continue;
162     if(ex_lep.charge != w_lep.charge) continue;
163     has_ssof = true;
164     ssof_leps.push_back(ex_lep);
165     types.push_back("pass");
166     }
167 dkralph 1.6 for(unsigned iextra=0; iextra<nfail; iextra++) {
168 dkralph 1.1 SimpleLepton ex_lep = (*fs->failingL)[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("fail");
174     }
175     if(!has_ssof) continue;
176    
177     // fill hists with W + (fake or pass)
178     double mt = get_mt(w_lep.vec.Phi(),fs->info->metphi,w_lep.vec.Pt(),fs->info->met);
179     (*(cs->hists)["both"])["w_pt"]->Fill( w_lep.vec.Pt(), wgt);
180     (*(cs->hists)["both"])["w_eta"]->Fill( w_lep.vec.Eta(), wgt);
181     (*(cs->hists)["both"])["w_isoPF04"]->Fill( w_lep.isoPF04/w_lep.vec.Pt(), wgt);
182     (*(cs->hists)["both"])["met"]->Fill( fs->info->met, wgt);
183     (*(cs->hists)["both"])["mt"]->Fill( mt, wgt);
184 dkralph 1.6 (*(cs->hists)["both"])["npass"]->Fill( npass, wgt);
185     (*(cs->hists)["both"])["nfail"]->Fill( nfail, wgt);
186     nBoth->Fill(npass,nfail,wgt);
187 dkralph 1.1
188     for(unsigned iextra=0; iextra<ssof_leps.size(); iextra++) {
189     SimpleLepton ex_lep = ssof_leps[iextra];
190     double mass = (w_lep.vec + ex_lep.vec).M();
191 dkralph 1.2
192 dkralph 1.1 if(types[iextra] == "fail") {
193    
194     // fill fake rate prediction
195 dkralph 1.5
196     double fake_wgt = get_fake_weight("",ex_lep,fr);
197     double fake_wgt_lo = get_fake_weight("lo",ex_lep,fr);
198     double fake_wgt_hi = get_fake_weight("hi",ex_lep,fr);
199    
200 dkralph 1.6 double centr=fake_wgt*wgt,lo=fake_wgt_lo*wgt,hi=fake_wgt_hi*wgt;
201     fillHist( cs, "pred", "w_pt", w_lep.vec.Pt(), centr, lo, hi);
202     fillHist( cs, "pred", "l2_pt", ex_lep.vec.Pt(), centr, lo, hi);
203     fillHist( cs, "pred", "w_eta", w_lep.vec.Eta(), centr, lo, hi);
204     fillHist( cs, "pred", "l2_eta", ex_lep.vec.Eta(), centr, lo, hi);
205     fillHist( cs, "pred", "w_isoPF04", w_lep.isoPF04/w_lep.vec.Pt(), centr, lo, hi);
206     fillHist( cs, "pred", "l2_isoPF04", ex_lep.isoPF04/ex_lep.vec.Pt(), centr, lo, hi);
207     fillHist( cs, "pred", "mass", mass, centr, lo, hi);
208     fillHist( cs, "pred", "met", fs->info->met, centr, lo, hi);
209     fillHist( cs, "pred", "mt", mt, centr, lo, hi);
210    
211 dkralph 1.3 (*(cs->hists)["pred"])["fr"]->Fill( fake_wgt );
212 dkralph 1.6
213 dkralph 1.1 } else if(types[iextra] == "pass") {
214    
215     // fill the observation
216    
217 dkralph 1.6 fillHist( cs, "obs", "w_pt", w_lep.vec.Pt(), wgt, 0, 0);
218     fillHist( cs, "obs", "l2_pt", ex_lep.vec.Pt(), wgt, 0, 0);
219     fillHist( cs, "obs", "w_eta", w_lep.vec.Eta(), wgt, 0, 0);
220     fillHist( cs, "obs", "l2_eta", ex_lep.vec.Eta(), wgt, 0, 0);
221     fillHist( cs, "obs", "w_isoPF04", w_lep.isoPF04/w_lep.vec.Pt(), wgt, 0, 0);
222     fillHist( cs, "obs", "l2_isoPF04", ex_lep.isoPF04/ex_lep.vec.Pt(), wgt, 0, 0);
223     fillHist( cs, "obs", "mass", mass, wgt, 0, 0);
224     fillHist( cs, "obs", "met", fs->info->met, wgt, 0, 0);
225     fillHist( cs, "obs", "mt", mt, wgt, 0, 0);
226 dkralph 1.1 } else { cout << "ERROR: types: " << types[iextra] << endl; assert(0); }
227     }
228     }
229     }
230     }
231     }
232    
233     assert(samplev.size()==3);
234     CSample *cs = samplev[0];
235     CSample *cs_mc = samplev[1];
236     CSample *cs_zmc = samplev[2];
237    
238     double total_mc_scale=-1;
239 dkralph 1.2 RooFitResult *fitres = fitMt(ctrl,(*(cs_zmc->hists)["both"])["mt"],(*(cs_mc->hists)["both"])["mt"],(*(cs->hists)["both"])["mt"],total_mc_scale);
240 dkralph 1.1
241 dkralph 1.6 cout << "\n\n\ntotal mc scale: " << total_mc_scale << endl << endl << endl;
242    
243 dkralph 1.1 // RooArgList float_pars = fitres->floatParsFinal();
244     // RooRealVar *nWandZ = (RooRealVar*)float_pars.at(float_pars.index("nWandZ"));
245     // cout << nWandZ->getVal() << endl;
246    
247     map<TString,TH1D*>::iterator it_v;
248     for(it_v=(*(cs->hists)["pred"]).begin(); it_v!=(*(cs->hists)["pred"]).end(); it_v++) {
249     TString var((*it_v).first);
250    
251 dkralph 1.6 // plot passing and failing fakes together
252 dkralph 1.2 if(var!="mt") {
253 dkralph 1.1 (*(cs_mc->hists)["both"])[var]->Scale(total_mc_scale);
254     (*(cs_zmc->hists)["both"])[var]->Scale(total_mc_scale);
255     }
256     CPlot cplot_both("both_"+var,"",(*(cs->hists)["both"])[var]->GetXaxis()->GetTitle(),"events",ctrl.outdir+"/plots");
257     cplot_both.AddHist1D((*(cs->hists)["both"])[var],"tight+fk'ble: "+integral_str((*(cs->hists)["both"])[var]),"E");
258     cplot_both.AddToStack((*(cs_zmc->hists)["both"])[var],"Z+j MC: "+integral_str((*(cs_zmc->hists)["both"])[var]),920);
259     cplot_both.AddToStack((*(cs_mc->hists)["both"])[var],"W+j MC: "+integral_str((*(cs_mc->hists)["both"])[var]),426);
260    
261     if(var.Contains("eta")) cplot_both.SetYRange(0,1.4*(*(cs->hists)["both"])[var]->GetMaximum());
262     if(var.Contains("isoPF")) cplot_both.SetLogy();
263 dkralph 1.2 if(var=="met") cplot_both.AddLine(25,0,25,18050,843,kDashed);
264 dkralph 1.1 cplot_both.Draw(can,true,"png");
265    
266     // plot observation vs. prediction
267     (*(cs_mc->hists)["obs"])[var]->Scale(total_mc_scale);
268     (*(cs_zmc->hists)["obs"])[var]->Scale(total_mc_scale);
269     CPlot cplot(var,"",(*(cs->hists)["pred"])[var]->GetXaxis()->GetTitle(),"events",ctrl.outdir+"/plots");
270     cplot.AddHist1D((*(cs->hists)["obs"])[var],"Obs: "+integral_str((*(cs->hists)["obs"])[var]),"E");
271     cplot.AddHist1D((*(cs->hists)["pred"])[var],"FR predic: "+integral_str((*(cs->hists)["pred"])[var]),"hist",kRed);
272 dkralph 1.6 cplot.AddHist1D((*(cs->hists)["pred_lo"])[var],"stat lo: "+integral_str((*(cs->hists)["pred_lo"])[var]),"hist",kRed,kDashed);
273     cplot.AddHist1D((*(cs->hists)["pred_hi"])[var],"stat hi: "+integral_str((*(cs->hists)["pred_hi"])[var]),"hist",kRed,kDashed);
274    
275     // stack on the mc
276 dkralph 1.5 // cplot.AddToStack((*(cs_zmc->hists)["obs"])[var],"Z+j MC: "+integral_str((*(cs_zmc->hists)["obs"])[var]),920);
277     // cplot.AddToStack((*(cs_mc->hists)["obs"])[var],"W+j MC: "+integral_str((*(cs_mc->hists)["obs"])[var]),426);
278 dkralph 1.1
279     if(var.Contains("eta")) cplot.SetYRange(0,1.4*max((*(cs->hists)["obs"])[var]->GetMaximum(),(*(cs->hists)["pred"])[var]->GetMaximum()));
280 dkralph 1.6 if(var.Contains("isoPF")) {
281     cplot.SetLogy();
282     TH1D *hobs = (*(cs->hists)["obs"])[var];
283     TH1D *hpred = (*(cs->hists)["pred"])[var];
284     double isoMin = 0.008*min(hobs->GetBinContent(1), hpred->GetBinContent(1));
285     double isoMax = 1.8*max(hobs->GetBinContent(1), hpred->GetBinContent(1));
286    
287     cplot.SetYRange(isoMin,isoMax);
288     // cplot.SetYRange(0,0.8*min((*(cs->hists)["obs"])[var]->GetMinimum(),(*(cs->hists)["pred"])[var]->GetMinimum()));
289     }
290 dkralph 1.1 cplot.Draw(can,true,"png");
291     }
292    
293     makeHTML(ctrl);
294 dkralph 1.6 nBoth->Draw("surf");
295     can->SaveAs("~/public_html/foo.png");
296     can->SaveAs("~/public_html/foo.C");
297     nBoth->Draw("colz");
298     can->SaveAs("~/public_html/fooc.png");
299     can->SaveAs("~/public_html/fooc.C");
300 dkralph 1.1 }
301     //----------------------------------------------------------------------------------------
302 dkralph 1.5 double get_fake_weight(TString uncert, SimpleLepton fake_lep, FR_struct fr)
303 dkralph 1.1 {
304     double rate=0;
305 dkralph 1.5 assert(uncert=="" || uncert=="lo" || uncert=="hi");
306 dkralph 1.1 if(abs(fake_lep.type) == 13) { // fake muons:
307 dkralph 1.5 double x = fake_lep.vec.Eta();
308     if(fr.absEtaMu) x = fabs(x);
309     double y = fake_lep.vec.Pt();
310     rate = fr.mufr.getEff(x,y,true);
311     if(uncert == "lo") rate -= fr.mufr.getErrLow(x,y,true);
312     if(uncert == "hi") rate += fr.mufr.getErrHigh(x,y,true);
313 dkralph 1.1 } else if(abs(fake_lep.type) == 11) {
314 dkralph 1.5 double x = fake_lep.vec.Eta();
315     if(fr.absEtaEle) x = fabs(x);
316     double y = fake_lep.vec.Pt();
317     rate = fr.elefr.getEff(x,y,true);
318     if(uncert == "lo") rate -= fr.elefr.getErrLow(x,y,true);
319     if(uncert == "hi") rate += fr.elefr.getErrHigh(x,y,true);
320 dkralph 1.1 } else assert(0);
321    
322     return rate/(1-rate);
323     }
324     //----------------------------------------------------------------------------------------
325     map<TString,map<TString,TH1D*>* > init_hists(FOFlags &ctrl, TString str)
326     {
327 dkralph 1.5 map<TString,map<TString,TH1D*>* > hists;
328     map<TString,TH1D*> *h_both = new map<TString,TH1D*>; hists["both"] = h_both;
329     map<TString,TH1D*> *h_pred = new map<TString,TH1D*>; hists["pred"] = h_pred;
330     map<TString,TH1D*> *h_pred_lo = new map<TString,TH1D*>; hists["pred_lo"] = h_pred_lo;
331     map<TString,TH1D*> *h_pred_hi = new map<TString,TH1D*>; hists["pred_hi"] = h_pred_hi;
332     map<TString,TH1D*> *h_obs = new map<TString,TH1D*>; hists["obs"] = h_obs;
333 dkralph 1.1 map<TString,map<TString,TH1D*>* >::iterator it_h;
334     double w_pt_max,l2_pt_max,mass_max;
335     if(ctrl.faketype=="mu") {
336     w_pt_max = 80;
337 dkralph 1.5 l2_pt_max = 35;
338 dkralph 1.1 mass_max = 120;
339     } else if(ctrl.faketype=="ele"){
340     w_pt_max = 90;
341     l2_pt_max = 60;
342     mass_max = 210;
343     } else assert(0);
344    
345     for(it_h=hists.begin(); it_h!=hists.end(); it_h++) {
346 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();
347     (*((*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();
348     (*((*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();
349     (*((*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();
350 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();
351     (*((*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();
352 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();
353     (*((*it_h).second))["met"] = new TH1D(TString("met") +"_"+(*it_h).first+str,";#bf{MET};", 25,0,100); (*((*it_h).second))["met"]->Sumw2();
354     (*((*it_h).second))["mt"] = new TH1D(TString("mt") +"_"+(*it_h).first+str,";#bf{m_{T}};", 25,2,150); (*((*it_h).second))["mt"]->Sumw2();
355     (*((*it_h).second))["fr"] = new TH1D(TString("fr") +"_"+(*it_h).first+str,";#bf{FR};", 125,-.01,.1); (*((*it_h).second))["fr"]->Sumw2();
356 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();
357     (*((*it_h).second))["nfail"] = new TH1D(TString("nfail") +"_"+(*it_h).first+str,";#bf{N Fail};", 125,-.01,4); (*((*it_h).second))["nfail"]->Sumw2();
358 dkralph 1.1 }
359    
360     return hists;
361     }
362     //--------------------------------------------------------------------------------------------------
363     void makeHTML(FOFlags &ctrl)
364     {
365     // TString title(ctrl.inputdir);
366     // title = title(title.Last('/')+1,title.Length());
367     TString title("EMU closure test: Fake "+ctrl.faketype);
368     ofstream htmlfile;
369     char htmlfname[100];
370     sprintf(htmlfname,"%s/"+ctrl.faketype+".html",ctrl.outdir.Data());
371     htmlfile.open(htmlfname);
372    
373     htmlfile << "<!DOCTYPE html" << endl;
374     htmlfile << " PUBLIC \"-//W3C//DTD HTML 3.2//EN\">" << endl;
375     htmlfile << "<html>" << endl;
376    
377     htmlfile << "<head><title>"+title+"</title></head>" << endl;
378     htmlfile << "<body bgcolor=\"EEEEEE\">" << endl;
379     htmlfile << "<h3 style=\"text-align:left; color:DD6600;\">"+title+"</h3>" << endl;
380    
381     htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;
382     htmlfile << "<tr>" << endl;
383 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;
384     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;
385 dkralph 1.1 htmlfile << "</tr>" << endl;
386     htmlfile << "</table>" << endl;
387    
388     htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;
389     htmlfile << "<tr>" << endl;
390 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;
391     htmlfile << "<td width=\"33.3%\"><a></a></td>" << endl;
392     htmlfile << "<td width=\"33.3%\"><a></a></td>" << endl;
393     htmlfile << "</tr>" << endl;
394     htmlfile << "</table>" << endl;
395    
396     htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;
397     htmlfile << "<tr>" << endl;
398 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;
399     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;
400     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;
401     htmlfile << "</tr>" << endl;
402     htmlfile << "</table>" << endl;
403    
404     htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;
405     htmlfile << "<tr>" << endl;
406     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;
407     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;
408     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;
409     htmlfile << "</tr>" << endl;
410     htmlfile << "</table>" << endl;
411    
412     htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;
413     htmlfile << "<tr>" << endl;
414     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;
415     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;
416     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;
417     htmlfile << "</tr>" << endl;
418     htmlfile << "</table>" << endl;
419     htmlfile << "<hr />" << endl;
420    
421 dkralph 1.6 // htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;
422     // htmlfile << "<tr>" << endl;
423     // 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;
424     // htmlfile << "</tr>" << endl;
425     // htmlfile << "</table>" << endl;
426 dkralph 1.1
427     htmlfile << "<hr />" << endl;
428    
429     htmlfile << "<hr />" << endl;
430    
431     htmlfile << "</body>" << endl;
432     htmlfile << "</html>" << endl;
433     htmlfile.close();
434     }
435     //----------------------------------------------------------------------------------------
436     TVector3 get_corrected_met(double met, double metphi, TLorentzVector lepton)
437     {
438     TLorentzVector corr_met;
439     corr_met.SetPtEtaPhiM(met,0,metphi,0);
440     TVector3 w_lep_trans(lepton.Px(),lepton.Py(),0);
441     return corr_met.Vect() + w_lep_trans;
442     }
443     //----------------------------------------------------------------------------------------
444     TString entries_str(TH1D *hist)
445     {
446     int entries = hist->GetEntries();
447     stringstream ss;
448     ss << entries;
449     TString str;
450     ss >> str;
451     return str;
452     }
453     //----------------------------------------------------------------------------------------
454     TString integral_str(TH1D *hist)
455     {
456     stringstream ss;
457 dkralph 1.2 ss << fixed << setprecision(1) << integrateHist(hist);
458 dkralph 1.1 TString str;
459     ss >> str;
460     return str;
461     }
462     //----------------------------------------------------------------------------------------
463 dkralph 1.2 RooFitResult *fitMt(FOFlags &ctrl, TH1D *hmcZ, TH1D *hmcW, TH1D *hdata, double &total_mc_scale)
464 dkralph 1.1 {
465 dkralph 1.2 double original_w_integral = integrateHist(hmcW);
466 dkralph 1.1
467 dkralph 1.2 RooRealVar mt("mt","MT [GeV]",hdata->GetXaxis()->GetXmin(),hdata->GetXaxis()->GetXmax());
468     RooBinning mtbins(hdata->GetNbinsX(),hdata->GetXaxis()->GetXmin(),hdata->GetXaxis()->GetXmax());
469     mt.setBinning(mtbins);
470 dkralph 1.1
471     // data
472 dkralph 1.2 RooDataHist obsMt("obsMt","obsMt",RooArgList(mt),hdata);
473 dkralph 1.1 // W mc pdf
474     hmcW->Smooth(2);
475 dkralph 1.2 RooDataHist hW("hW","hW",RooArgList(mt),hmcW);
476     RooHistPdf pW("pW","pW",RooArgSet(mt),hW,3);
477     // Z mc pdf
478     RooDataHist hZ("hZ","hZ",RooArgList(mt),hmcZ);
479     RooHistPdf pZ("pZ","pZ",RooArgSet(mt),hZ,3);
480 dkralph 1.1 // rayleigh pdf
481 dkralph 1.2 RooRealVar sig("sig","sig",400,15,585);
482     RooRealVar a1("a1","a1",-3,-10,10);
483     // RooGenericPdf rayl("rayl","rayl","mt*exp(-(mt)*(mt)/(2*sig*sig))",RooArgList(mt,sig));
484     RooGenericPdf rayl("rayl","rayl","mt*exp(-(mt)*(mt)/(sig*sig +2*sig*a1*(mt) + a1*a1*(mt)*(mt)))",RooArgList(mt,sig,a1));
485    
486     RooRealVar nObs("nObs","nObs",integrateHist(hdata));
487     RooRealVar nWandZ("nWandZ","nWandZ",0.98*nObs.getVal(),.95*nObs.getVal(),1.1*nObs.getVal());
488     RooRealVar zFrac("zFrac","zFrac",integrateHist(hmcZ) / (integrateHist(hmcZ) + integrateHist(hmcW)));
489     RooAddPdf modelWandZ("modelWandZ","modelWandZ",RooArgList(pZ,pW),zFrac);
490     RooRealVar nRayl("nRayl","nRayl",0.05*nObs.getVal(),0.0*nObs.getVal(),.05*nObs.getVal());
491    
492     RooExtendPdf extWandZ("extWandZ","extWandZ",modelWandZ,nWandZ);
493     RooExtendPdf extRayl("extRayl","extRayl",rayl,nRayl);
494    
495     RooAddPdf *model = new RooAddPdf("model","model",RooArgList(extWandZ,extRayl));
496    
497     RooFitResult *fitres = model->fitTo(obsMt,Extended(kTRUE),Save(kTRUE));
498     fitres->Print();
499    
500     TH1D *rayl_hist = (TH1D*)rayl.createHistogram("rayl_mt",mt,Binning(mtbins));
501     rayl_hist->Scale(nRayl.getVal()/integrateHist(rayl_hist));
502     TH1D *totalHist = (TH1D*)model->createHistogram("total_mt",mt,Binning(mtbins));
503     totalHist->Scale( (nWandZ.getVal() + nRayl.getVal()) / integrateHist(totalHist) );
504    
505     // // why the fuck do I need this?
506     // double tot = nWandZ.getVal() + nRayl.getVal();
507     // if(use_exp) tot += nExp.getVal();
508     // double hack_fac = integrateHist(hdata)/tot;
509     // rayl_hist->Scale(hack_fac);
510     // exp_hist->Scale(hack_fac);
511     // hmcZ->Scale(hack_fac);
512     // hmcW->Scale(hack_fac);
513    
514     // RooPlot *frame = mt.frame();
515     // obsMt.plotOn(frame);
516     // model->plotOn(frame);
517     // model->plotOn(frame,Components("pW"),LineStyle(kDashed),LineColor(kRed));
518     // model->plotOn(frame,Components("pZ"),LineColor(kRed));
519     // model->plotOn(frame,Components("rayl"),LineColor(kGreen));
520     // frame->Draw();
521     // can->SaveAs("~/public_html/foo.png");
522    
523     hmcW->Scale( (nWandZ.getVal()*(1-zFrac.getVal())) / integrateHist(hmcW) );
524     hmcZ->Scale( (nWandZ.getVal()*zFrac.getVal()) / integrateHist(hmcZ) );
525 dkralph 1.1
526 dkralph 1.2 CPlot plotMtFit("mtfitW","","MT [GeV]","Events",ctrl.outdir+"/plots");
527     plotMtFit.AddHist1D(hdata,"tight+fk'ble: " + integral_str(hdata),"E");
528     plotMtFit.AddHist1D(totalHist,"total MC: " + integral_str(totalHist),"hist",kRed);
529     plotMtFit.AddToStack(rayl_hist,"jet (Rayleigh): " + integral_str(rayl_hist),kBlue);
530     plotMtFit.AddToStack(hmcZ,"Z+j MC: " + integral_str(hmcZ),kGray);
531     plotMtFit.AddToStack(hmcW,"W+j MC: " + integral_str(hmcW),kCyan-6);
532     plotMtFit.Draw(can,kTRUE,"png");
533 dkralph 1.1
534 dkralph 1.2 total_mc_scale = integrateHist(hmcW) / original_w_integral;
535 dkralph 1.1
536     return fitres;
537     }
538 dkralph 1.6 //----------------------------------------------------------------------------------------
539     void fillHist(CSample *cs, TString type, TString var, double val, double wgt, double wgt_lo, double wgt_hi)
540     {
541     (*(cs->hists)[type])[var]->Fill( val, wgt);
542     if(type=="pred") {
543     (*(cs->hists)[type+"_lo"])[var]->Fill( val, wgt_lo);
544     (*(cs->hists)[type+"_hi"])[var]->Fill( val, wgt_hi);
545     }
546     }