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