ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/NonMCBackground/src/plot_emu_fakes.cc
Revision: 1.13
Committed: Tue Oct 23 11:22:53 2012 UTC (12 years, 6 months ago) by dkralph
Content type: text/plain
Branch: MAIN
CVS Tags: compiled, HEAD
Changes since 1.12: +7 -7 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.9 #include "PlotHeaders.h"
32 dkralph 1.1
33     #ifndef CMSSW_BASE
34     #define CMSSW_BASE "../../"
35     #endif
36    
37     using namespace std;
38     using namespace RooFit;
39 dkralph 1.11 using namespace mithep;
40 dkralph 1.1
41 dkralph 1.12
42     TH1D* hpu_2011;
43     TH1D* hpu_2012;
44     TH1D* hpu_f11_to_2012;
45     TH1D* hpu_f11_to_1112;
46     TH1D* hpu_s12_to_1112;
47    
48 dkralph 1.1 TCanvas *can;
49    
50 dkralph 1.9 void makeHTML(FOFlags &ctrl, TString plotLabel);
51 dkralph 1.1 map<TString,map<TString,TH1D*>* > init_hists(FOFlags &ctrl, TString str="");
52 dkralph 1.6 void fillHist(CSample *cs, TString type, TString var, double val, double wgt, double wgt_lo, double wgt_hi);
53 dkralph 1.12 void scaleHist(TH1D *hist, double scale)
54     {
55     for(int ibin=0; ibin<hist->GetNbinsX()+2; ibin++) {
56     hist->SetBinContent(ibin,hist->GetBinContent(ibin)*scale);
57     }
58     }
59 dkralph 1.1
60 dkralph 1.9 RooFitResult *fitMt(FOFlags &ctrl, TH1D *hmcZ, TH1D *hmcW, TH1D *hSSmet, double &total_mc_scale, TString plotLabel);
61 dkralph 1.1
62     //----------------------------------------------------------------------------------------
63     int main(int argc, char** argv)
64     {
65 dkralph 1.12 initPUWeights();
66     TFile *f_pu = TFile::Open("data/puWeights/f11_to_2012.root");
67     hpu_f11_to_2012 = (TH1D*)f_pu->Get("puWeights");
68     hpu_f11_to_2012->SetDirectory(0);
69     f_pu->Close();
70     f_pu = TFile::Open("data/puWeights/f11_to_1112.root");
71     hpu_f11_to_1112 = (TH1D*)f_pu->Get("puWeights");
72     hpu_f11_to_1112->SetDirectory(0);
73     f_pu->Close();
74     f_pu = TFile::Open("data/puWeights/s12-zjets_to_1112.root");
75     hpu_s12_to_1112 = (TH1D*)f_pu->Get("puWeights");
76     hpu_s12_to_1112->SetDirectory(0);
77     f_pu->Close();
78 dkralph 1.6
79 dkralph 1.1 FOFlags ctrl;
80     parse_foargs( argc, argv, ctrl );
81 dkralph 1.4 ctrl.dump();
82 dkralph 1.1
83     can = new TCanvas("can","can");
84    
85 dkralph 1.9 FR_struct fr2011 = initFRs(ctrl.mufakefile2011,ctrl.elefakefile2011);
86     FR_struct fr2012 = initFRs(ctrl.mufakefile2012,ctrl.elefakefile2012);
87 dkralph 1.1
88     TString cmsswpath(CMSSW_BASE + TString("/src"));
89     string xspath = (string(cmsswpath)+"/MitPhysics/data/xs.dat");
90     SimpleTable xstab(xspath.c_str());
91    
92     vector<CSample*> samplev;
93 dkralph 1.11 TString ntupledir(""),label(""),plotLabel(""),jsonFile("");
94 dkralph 1.12 int puTarget;
95     readConfigFile(ctrl.config, ntupledir, label, plotLabel, jsonFile, puTarget, samplev, &ctrl, init_hists);
96 dkralph 1.1
97 dkralph 1.8 vector<pair<unsigned,unsigned> > runEvtv; // vector to veto duplicate events
98 dkralph 1.9 UInt_t minRun=999999999,maxRun=0;
99 dkralph 1.1 for(unsigned ics=0; ics<samplev.size(); ics++) {
100     CSample *cs = samplev[ics];
101 dkralph 1.8 cout << cs->name << ": " << endl;
102 dkralph 1.1 for(unsigned ifs=0; ifs<(cs->fsv).size(); ifs++) {
103     filestuff *fs = (cs->fsv)[ifs];
104 dkralph 1.8 cout << "\t" << fs->fname_ << endl;
105 dkralph 1.9 FR_struct *fr = (fs->era_==2011) ? &fr2011 : &fr2012;
106     unsigned nDuplSkipped=0;
107 dkralph 1.8 for(unsigned ientry=0; ientry<fs->getentries("FOtree"); ientry++) {
108 dkralph 1.1 fs->getentry(ientry,"info");
109 dkralph 1.9
110     if(fs->isdata_) {
111     setMinMaxRun(fs->info->run, minRun, maxRun);
112     // bool dupl = takeCareOfDuplicateEvents(fs->info->run, fs->info->evt, runEvtv, nDuplSkipped);
113     // if(dupl) continue;
114 dkralph 1.11 pair<unsigned,unsigned> rl(fs->info->run, fs->info->lumi);
115     if(!fs->rlrm_.HasRunLumi(rl)) continue;
116 dkralph 1.9 }
117 dkralph 1.1
118     unsigned npass = fs->passingL->size();
119     unsigned nfail = fs->failingL->size();
120    
121 dkralph 1.6 // z veto
122     bool zVeto=false;
123     for(unsigned ilep1=0; ilep1<npass; ilep1++) {
124     SimpleLepton lep1 = (*fs->passingL)[ilep1];
125     for(unsigned ilep2=ilep1+1; ilep2<npass; ilep2++) {
126     SimpleLepton lep2 = (*fs->passingL)[ilep2];
127 dkralph 1.7 // if(abs(lep1.type) == abs(lep2.type)) zVeto = true;
128     // if(lep1.charge != lep2.charge) zVeto = true;
129     if(abs(lep1.type) == abs(lep2.type) &&
130     lep1.charge != lep2.charge) zVeto = true;
131 dkralph 1.6 }
132     }
133     if(zVeto) continue;
134    
135     for(unsigned iw=0; iw<npass; iw++) {
136 dkralph 1.1 SimpleLepton w_lep = (*fs->passingL)[iw];
137    
138 dkralph 1.7 double wgt=1;
139     if(!fs->isdata_) {
140 dkralph 1.12 double hackFactor = 1;//(ctrl.faketype=="mu") ? 1/*7.42*/ : 4.71;
141     double xsWgt = hackFactor*fs->lumi_*xstab.Get(fs->dataset_)/fs->total_entries_;
142     TH1D *hpu=0;
143     if(puTarget==1112 && fs->era_==2011) hpu = hpu_f11_to_1112;
144     else if(puTarget==1112 && fs->era_==2012) hpu = hpu_s12_to_1112;
145     else hpu = hpu_2012;
146     // if(puTarget==2011) hpu = hpu_2011;
147     // else if(puTarget==2012 && fs->era_==2012) hpu = hpu_2011;
148     // else if(puTarget==2012 && fs->era_==2011) hpu = hpu_f11_to_2012;
149     // else if(puTarget==1112 && fs->era_==2011) hpu = hpu_f11_to_1112;
150     // else if(puTarget==1112 && fs->era_==2012) hpu = hpu_s12_to_1112;
151     // else assert(0);
152    
153     double puWgt = 1;//hpu->GetBinContent(hpu->FindBin(fs->info->npu));
154 dkralph 1.7 wgt = xsWgt*puWgt;
155     }
156 dkralph 1.1
157 dkralph 1.6 // make fake ele and fake mu plots separately
158 dkralph 1.8 assert(ctrl.faketype=="mu" || ctrl.faketype=="ele");
159     assert(abs(w_lep.type)==11 || abs(w_lep.type)==13);
160     if( abs(w_lep.type)==11 && ctrl.faketype=="ele") continue;
161     if( abs(w_lep.type)==13 && ctrl.faketype=="mu" ) continue;
162 dkralph 1.7
163 dkralph 1.1 // selection on the W lepton
164 dkralph 1.13 // if(w_lep.status.isoPF04/w_lep.vec.Pt() > 0.025) continue;
165 dkralph 1.6 if(fs->info->met < 25) continue;
166 dkralph 1.2 if(!(w_lep.tightCutsApplied)) continue;
167 dkralph 1.6 if(w_lep.vec.Pt() < 25) continue;
168 dkralph 1.8 if(fabs(w_lep.ip3dSig) > 4) continue;
169     // if( fs->isdata_)
170     // if( !(w_lep.bdtfail & 2) ) continue;
171 dkralph 1.7 // if(!w_lep.isTight) continue; // this bit stores whether the w lepton is trigger matched
172 dkralph 1.1
173     // require an extra OF SS lepton
174     bool has_ssof=false;
175     vector<SimpleLepton> ssof_leps;
176     vector<TString> types;
177 dkralph 1.6 for(unsigned iextra=iw+1; iextra<npass; iextra++) {
178 dkralph 1.1 SimpleLepton ex_lep = (*fs->passingL)[iextra];
179     if(abs(ex_lep.type) == abs(w_lep.type)) continue;
180     if(ex_lep.charge != w_lep.charge) continue;
181     has_ssof = true;
182     ssof_leps.push_back(ex_lep);
183     types.push_back("pass");
184     }
185 dkralph 1.6 for(unsigned iextra=0; iextra<nfail; iextra++) {
186 dkralph 1.1 SimpleLepton ex_lep = (*fs->failingL)[iextra];
187     if(abs(ex_lep.type) == abs(w_lep.type)) continue;
188     if(ex_lep.charge != w_lep.charge) continue;
189     has_ssof = true;
190     ssof_leps.push_back(ex_lep);
191     types.push_back("fail");
192     }
193     if(!has_ssof) continue;
194    
195     for(unsigned iextra=0; iextra<ssof_leps.size(); iextra++) {
196     SimpleLepton ex_lep = ssof_leps[iextra];
197     double mass = (w_lep.vec + ex_lep.vec).M();
198 dkralph 1.2
199 dkralph 1.12 // fill hists with W + (fake or pass)
200     double mt = get_mt(w_lep.vec.Phi(),fs->info->metphi,w_lep.vec.Pt(),fs->info->met);
201     (*(cs->hists)["both"])["run"]->Fill( fs->info->run, wgt);
202     (*(cs->hists)["both"])["npv"]->Fill( fs->info->npv, wgt);
203     (*(cs->hists)["both"])["w_pt"]->Fill( w_lep.vec.Pt(), wgt);
204     (*(cs->hists)["both"])["l2_pt"]->Fill( ex_lep.vec.Pt(), wgt);
205     (*(cs->hists)["both"])["w_eta"]->Fill( w_lep.vec.Eta(), wgt);
206     (*(cs->hists)["both"])["l2_eta"]->Fill( ex_lep.vec.Eta(), wgt);
207 dkralph 1.13 (*(cs->hists)["both"])["w_isoPF04"]->Fill( w_lep.status.isoPF04/w_lep.vec.Pt(), wgt);
208     (*(cs->hists)["both"])["l2_isoPF04"]->Fill( ex_lep.status.isoPF04/w_lep.vec.Pt(), wgt);
209 dkralph 1.12 (*(cs->hists)["both"])["met"]->Fill( fs->info->met, wgt);
210     (*(cs->hists)["both"])["mt"]->Fill( mt, wgt);
211     (*(cs->hists)["both"])["npass"]->Fill( npass, wgt);
212     (*(cs->hists)["both"])["nfail"]->Fill( nfail, wgt);
213     (*(cs->hists)["both"])["w_ip3ds_lo"]->Fill( w_lep.ip3dSig, wgt);
214     (*(cs->hists)["both"])["l2_ip3ds_lo"]->Fill( ex_lep.ip3dSig, wgt);
215     (*(cs->hists)["both"])["w_ip3ds"]->Fill( w_lep.ip3dSig, wgt);
216     (*(cs->hists)["both"])["l2_ip3ds"]->Fill( ex_lep.ip3dSig, wgt);
217     (*(cs->hists)["both"])["w_d0"]->Fill( w_lep.d0, wgt);
218     (*(cs->hists)["both"])["l2_d0"]->Fill( ex_lep.d0, wgt);
219     (*(cs->hists)["both"])["w_dz"]->Fill( w_lep.dz, wgt);
220     (*(cs->hists)["both"])["l2_dz"]->Fill( ex_lep.dz, wgt);
221    
222     // if(fabs(ex_lep.ip3dSig) > 4) continue;
223 dkralph 1.1 if(types[iextra] == "fail") {
224    
225     // fill fake rate prediction
226 dkralph 1.5
227 dkralph 1.9 double fake_wgt = get_fake_weight("",ex_lep,*fr);
228     double fake_wgt_lo = get_fake_weight("lo",ex_lep,*fr);
229     double fake_wgt_hi = get_fake_weight("hi",ex_lep,*fr);
230 dkralph 1.5
231 dkralph 1.6 double centr=fake_wgt*wgt,lo=fake_wgt_lo*wgt,hi=fake_wgt_hi*wgt;
232 dkralph 1.12 fillHist( cs, "pred", "run", fs->info->run, centr, lo, hi);
233 dkralph 1.6 fillHist( cs, "pred", "w_pt", w_lep.vec.Pt(), centr, lo, hi);
234     fillHist( cs, "pred", "l2_pt", ex_lep.vec.Pt(), centr, lo, hi);
235     fillHist( cs, "pred", "w_eta", w_lep.vec.Eta(), centr, lo, hi);
236     fillHist( cs, "pred", "l2_eta", ex_lep.vec.Eta(), centr, lo, hi);
237 dkralph 1.13 fillHist( cs, "pred", "w_isoPF04", w_lep.status.isoPF04/w_lep.vec.Pt(), centr, lo, hi);
238     fillHist( cs, "pred", "l2_isoPF04", ex_lep.status.isoPF04/ex_lep.vec.Pt(), centr, lo, hi);
239 dkralph 1.8 fillHist( cs, "pred", "w_ip3ds", w_lep.ip3dSig, centr, lo, hi);
240     fillHist( cs, "pred", "l2_ip3ds", ex_lep.ip3dSig, centr, lo, hi);
241 dkralph 1.12 fillHist( cs, "pred", "w_ip3ds_lo", w_lep.ip3dSig, centr, lo, hi);
242     fillHist( cs, "pred", "l2_ip3ds_lo", ex_lep.ip3dSig, centr, lo, hi);
243     fillHist( cs, "pred", "w_d0", w_lep.d0, centr, lo, hi);
244     fillHist( cs, "pred", "l2_d0", ex_lep.d0, centr, lo, hi);
245     fillHist( cs, "pred", "w_dz", w_lep.dz, centr, lo, hi);
246     fillHist( cs, "pred", "l2_dz", ex_lep.dz, centr, lo, hi);
247 dkralph 1.6 fillHist( cs, "pred", "mass", mass, centr, lo, hi);
248     fillHist( cs, "pred", "met", fs->info->met, centr, lo, hi);
249     fillHist( cs, "pred", "mt", mt, centr, lo, hi);
250    
251 dkralph 1.3 (*(cs->hists)["pred"])["fr"]->Fill( fake_wgt );
252 dkralph 1.6
253 dkralph 1.1 } else if(types[iextra] == "pass") {
254    
255     // fill the observation
256 dkralph 1.12 fillHist( cs, "obs", "run", fs->info->run, wgt, 0, 0);
257 dkralph 1.6 fillHist( cs, "obs", "w_pt", w_lep.vec.Pt(), wgt, 0, 0);
258     fillHist( cs, "obs", "l2_pt", ex_lep.vec.Pt(), wgt, 0, 0);
259     fillHist( cs, "obs", "w_eta", w_lep.vec.Eta(), wgt, 0, 0);
260     fillHist( cs, "obs", "l2_eta", ex_lep.vec.Eta(), wgt, 0, 0);
261 dkralph 1.13 fillHist( cs, "obs", "w_isoPF04", w_lep.status.isoPF04/w_lep.vec.Pt(), wgt, 0, 0);
262     fillHist( cs, "obs", "l2_isoPF04", ex_lep.status.isoPF04/ex_lep.vec.Pt(), wgt, 0, 0);
263 dkralph 1.8 fillHist( cs, "obs", "w_ip3ds", w_lep.ip3dSig, wgt, 0, 0);
264     fillHist( cs, "obs", "l2_ip3ds", ex_lep.ip3dSig, wgt, 0, 0);
265 dkralph 1.12 fillHist( cs, "obs", "w_ip3ds_lo", w_lep.ip3dSig, wgt, 0, 0);
266     fillHist( cs, "obs", "l2_ip3ds_lo", ex_lep.ip3dSig, wgt, 0, 0);
267     fillHist( cs, "obs", "w_d0", w_lep.d0, wgt, 0, 0);
268     fillHist( cs, "obs", "l2_d0", ex_lep.d0, wgt, 0, 0);
269     fillHist( cs, "obs", "w_dz", w_lep.dz, wgt, 0, 0);
270     fillHist( cs, "obs", "l2_dz", ex_lep.dz, wgt, 0, 0);
271 dkralph 1.6 fillHist( cs, "obs", "mass", mass, wgt, 0, 0);
272     fillHist( cs, "obs", "met", fs->info->met, wgt, 0, 0);
273     fillHist( cs, "obs", "mt", mt, wgt, 0, 0);
274 dkralph 1.1 } else { cout << "ERROR: types: " << types[iextra] << endl; assert(0); }
275     }
276     }
277     }
278 dkralph 1.9 cout << "\t\tWARNING: not checking for duplicate events" << endl; //skipped " << nDuplSkipped << " duplicate events" << endl;
279 dkralph 1.1 }
280     }
281 dkralph 1.9 cout << "run range: " << setw(12) << minRun << setw(12) << maxRun << endl;
282 dkralph 1.1
283     assert(samplev.size()==3);
284     CSample *cs = samplev[0];
285     CSample *cs_mc = samplev[1];
286     CSample *cs_zmc = samplev[2];
287    
288     double total_mc_scale=-1;
289 dkralph 1.9 RooFitResult *fitres = fitMt(ctrl,(*(cs_zmc->hists)["both"])["mt"],(*(cs_mc->hists)["both"])["mt"],(*(cs->hists)["both"])["mt"],total_mc_scale,plotLabel);
290 dkralph 1.1
291 dkralph 1.6 cout << "\n\n\ntotal mc scale: " << total_mc_scale << endl << endl << endl;
292    
293 dkralph 1.1 // RooArgList float_pars = fitres->floatParsFinal();
294     // RooRealVar *nWandZ = (RooRealVar*)float_pars.at(float_pars.index("nWandZ"));
295     // cout << nWandZ->getVal() << endl;
296    
297 dkralph 1.9 gSystem->mkdir(ctrl.outdir+"/"+plotLabel+"/"+ctrl.faketype,true);
298     TFile runHistFile(ctrl.outdir+"/"+plotLabel+"/"+ctrl.faketype+"/runs.root","recreate");
299 dkralph 1.1 map<TString,TH1D*>::iterator it_v;
300 dkralph 1.12 for(it_v=(*(samplev[0]->hists)["both"]).begin(); it_v!=(*(samplev[0]->hists)["both"]).end(); it_v++) {
301 dkralph 1.1 TString var((*it_v).first);
302 dkralph 1.12 CPlot cplotBoth("both_"+var,"",(*(samplev[0]->hists)["both"])[var]->GetXaxis()->GetTitle(),"events",ctrl.outdir+"/"+plotLabel+"/"+ctrl.faketype+"/plots");
303     CPlot cplot(var,"",(*(samplev[0]->hists)["both"])[var]->GetXaxis()->GetTitle(),"events",ctrl.outdir+"/"+plotLabel+"/"+ctrl.faketype+"/plots");
304     double ymax=0,ymaxBoth=0;
305     for(unsigned isam=0; isam<samplev.size(); isam++) {
306     CSample *cs = samplev[isam];
307     TH1D *hBoth = (*(cs->hists)["both"])[var];
308     if(var!="mt") { // hack to avoid worrying about prescales {
309     // hBoth->Scale(total_mc_scale);
310     scaleHist(hBoth,total_mc_scale);
311     }
312     TH1D *hObs = (*(cs->hists)["obs"])[var];
313     TH1D *hPred = (*(cs->hists)["pred"])[var];
314     TH1D *hPred_lo = (*(cs->hists)["pred_lo"])[var];
315     TH1D *hPred_hi = (*(cs->hists)["pred_hi"])[var];
316     if((var.Contains("ip3d")&&(!var.Contains("_lo"))) || (var.Contains("d0")&&(!var.Contains("_lo"))) || var.Contains("dz") || var=="m4l" || var=="mZ2" || var.Contains("isopf",TString::kIgnoreCase)) {
317     shiftOverflows(hBoth);
318     shiftOverflows(hObs);
319     shiftOverflows(hPred);
320     shiftOverflows(hPred_lo);
321     shiftOverflows(hPred_hi);
322     }
323     if(cs->isdata) {
324     cplotBoth.AddHist1D(hBoth,"tight+fk'ble: "+integral_str(hBoth),"E");
325     ymaxBoth = max(ymaxBoth,hBoth->GetMaximum());
326     } else {
327     cplotBoth.AddToStack(hBoth,cs->legend+": "+integral_str(hBoth),cs->color);
328     ymaxBoth = max(ymaxBoth,hBoth->GetMaximum());
329     }
330     // if(var=="met") cplot_both.AddLine(25,0,25,18050,843,kDashed);
331 dkralph 1.1
332 dkralph 1.12 if(cs->isdata) {
333     double ratio = integrateHist(hObs) / integrateHist(hPred);
334     cplot.AddHist1D(hObs,cs->legend+": "+integral_str(hObs)+", ratio: "+f_to_a(ratio),"E");
335     cplot.AddHist1D(hPred,"FR predic: "+integral_str(hPred),"hist",kRed);
336     cplot.AddHist1D(hPred_lo,"stat lo: "+integral_str(hPred_lo),"hist",kRed,kDashed);
337     cplot.AddHist1D(hPred_hi,"stat hi: "+integral_str(hPred_hi),"hist",kRed,kDashed);
338     ymax = max(ymax,histMax(hObs,hPred,hPred_lo,hPred_hi));
339     } else {
340     TH1D *mchist = hObs;
341     cplot.AddToStack(mchist,cs->legend+": "+integral_str(mchist),cs->color);
342     // cplot.AddHist1D(mchist,cs->legend+": "+integral_str(mchist),"Ehist",cs->color);
343     ymax = max(ymax,mchist->GetMaximum());
344     }
345     if(cs->isdata && var=="run") {
346     assert(hBoth->GetBinContent(0)==0 && hBoth->GetBinContent(hBoth->GetXaxis()->GetNbins()+1)==0);
347     hBoth->Write("runs");
348     }
349 dkralph 1.1
350 dkralph 1.12 if(var.Contains("ip3d")) {
351     // double isoMin = 0.0005*ymax;
352     // double isoMax = 1.8*ymax;
353     // cplot.SetYRange(isoMin,isoMax);
354     // cplot.SetLogy();
355     cplot.SetYRange(0,1.2*ymax);
356     cplotBoth.SetYRange(0,1.2*ymaxBoth);
357     } else {
358     cplot.SetYRange(0,1.2*ymax);
359     cplotBoth.SetYRange(0,1.2*ymaxBoth);
360     }
361 dkralph 1.6 }
362 dkralph 1.12 cplotBoth.Draw(can,true,"png");
363 dkralph 1.1 cplot.Draw(can,true,"png");
364     }
365 dkralph 1.8 runHistFile.Close();
366 dkralph 1.1
367 dkralph 1.9 makeHTML(ctrl,plotLabel);
368 dkralph 1.1 }
369     //----------------------------------------------------------------------------------------
370     map<TString,map<TString,TH1D*>* > init_hists(FOFlags &ctrl, TString str)
371     {
372 dkralph 1.5 map<TString,map<TString,TH1D*>* > hists;
373     map<TString,TH1D*> *h_both = new map<TString,TH1D*>; hists["both"] = h_both;
374     map<TString,TH1D*> *h_pred = new map<TString,TH1D*>; hists["pred"] = h_pred;
375     map<TString,TH1D*> *h_pred_lo = new map<TString,TH1D*>; hists["pred_lo"] = h_pred_lo;
376     map<TString,TH1D*> *h_pred_hi = new map<TString,TH1D*>; hists["pred_hi"] = h_pred_hi;
377     map<TString,TH1D*> *h_obs = new map<TString,TH1D*>; hists["obs"] = h_obs;
378 dkralph 1.1 map<TString,map<TString,TH1D*>* >::iterator it_h;
379     double w_pt_max,l2_pt_max,mass_max;
380     if(ctrl.faketype=="mu") {
381     w_pt_max = 80;
382 dkralph 1.5 l2_pt_max = 35;
383 dkralph 1.1 mass_max = 120;
384     } else if(ctrl.faketype=="ele"){
385     w_pt_max = 90;
386     l2_pt_max = 60;
387     mass_max = 210;
388     } else assert(0);
389    
390     for(it_h=hists.begin(); it_h!=hists.end(); it_h++) {
391 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();
392     (*((*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();
393     (*((*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();
394     (*((*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();
395     (*((*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();
396 dkralph 1.12 (*((*it_h).second))["w_isoPF04"] = new TH1D(TString("w_isoPF04")+"_"+(*it_h).first+str,";#bf{real PF Iso};", 25,0,.2); (*((*it_h).second))["w_isoPF04"]->Sumw2();
397 dkralph 1.8 (*((*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();
398 dkralph 1.12
399     (*((*it_h).second))["w_ip3ds_lo"] = new TH1D(TString("w_ip3ds_lo")+"_"+(*it_h).first+str,";#bf{real ip3dSig};",35,-3,3); (*((*it_h).second))["w_ip3ds_lo"]->Sumw2();
400     (*((*it_h).second))["l2_ip3ds_lo"] = new TH1D(TString("l2_ip3ds_lo")+"_"+(*it_h).first+str,";#bf{fake ip3dSig};",35,-3,3); (*((*it_h).second))["l2_ip3ds_lo"]->Sumw2();
401     (*((*it_h).second))["w_ip3ds"] = new TH1D(TString("w_ip3ds")+"_"+(*it_h).first+str,";#bf{real ip3dSig};", 35,-3,3); (*((*it_h).second))["w_ip3ds"]->Sumw2();
402     (*((*it_h).second))["l2_ip3ds"] = new TH1D(TString("l2_ip3ds")+"_"+(*it_h).first+str,";#bf{fake ip3dSig};", 35,-6,6); (*((*it_h).second))["l2_ip3ds"]->Sumw2();
403     (*((*it_h).second))["w_d0"] = new TH1D(TString("w_d0")+"_"+(*it_h).first+str,";#bf{real d0};", 35,-.004,.004); (*((*it_h).second))["w_d0"]->Sumw2();
404     (*((*it_h).second))["l2_d0"] = new TH1D(TString("l2_d0")+"_"+(*it_h).first+str,";#bf{fake d0};", 35,-.03,.03); (*((*it_h).second))["l2_d0"]->Sumw2();
405     (*((*it_h).second))["w_dz"] = new TH1D(TString("w_dz")+"_"+(*it_h).first+str,";#bf{real dz};", 35,-.008,.008); (*((*it_h).second))["w_dz"]->Sumw2();
406     (*((*it_h).second))["l2_dz"] = new TH1D(TString("l2_dz")+"_"+(*it_h).first+str,";#bf{fake dz};", 35,-.03,.03); (*((*it_h).second))["l2_dz"]->Sumw2();
407    
408 dkralph 1.8 (*((*it_h).second))["mass"] = new TH1D(TString("mass") +"_"+(*it_h).first+str,";#bf{m_{ll}};", 25,0,mass_max); (*((*it_h).second))["mass"]->Sumw2();
409     (*((*it_h).second))["met"] = new TH1D(TString("met") +"_"+(*it_h).first+str,";#bf{MET};", 25,0,100); (*((*it_h).second))["met"]->Sumw2();
410     (*((*it_h).second))["mt"] = new TH1D(TString("mt") +"_"+(*it_h).first+str,";#bf{m_{T}};", 25,2,150); (*((*it_h).second))["mt"]->Sumw2();
411     (*((*it_h).second))["fr"] = new TH1D(TString("fr") +"_"+(*it_h).first+str,";#bf{FR};", 125,-.01,.1); (*((*it_h).second))["fr"]->Sumw2();
412     (*((*it_h).second))["npass"] = new TH1D(TString("npass") +"_"+(*it_h).first+str,";#bf{N Pass};", 125,-.01,4); (*((*it_h).second))["npass"]->Sumw2();
413     (*((*it_h).second))["nfail"] = new TH1D(TString("nfail") +"_"+(*it_h).first+str,";#bf{N Fail};", 125,-.01,4); (*((*it_h).second))["nfail"]->Sumw2();
414     (*((*it_h).second))["run"] = new TH1D(TString("run") +"_"+(*it_h).first+str,";#bf{Run};", 50,160405,196550); (*((*it_h).second))["run"]->Sumw2();
415 dkralph 1.1 }
416    
417     return hists;
418     }
419     //--------------------------------------------------------------------------------------------------
420 dkralph 1.9 void makeHTML(FOFlags &ctrl, TString plotLabel)
421 dkralph 1.1 {
422     // TString title(ctrl.inputdir);
423     // title = title(title.Last('/')+1,title.Length());
424     TString title("EMU closure test: Fake "+ctrl.faketype);
425 dkralph 1.9 TString htmlfname(ctrl.outdir+"/"+plotLabel+"/"+ctrl.faketype+"/plots.html");
426     ofstream htmlfile(htmlfname);;
427 dkralph 1.1
428     htmlfile << "<!DOCTYPE html" << endl;
429     htmlfile << " PUBLIC \"-//W3C//DTD HTML 3.2//EN\">" << endl;
430     htmlfile << "<html>" << endl;
431    
432     htmlfile << "<head><title>"+title+"</title></head>" << endl;
433     htmlfile << "<body bgcolor=\"EEEEEE\">" << endl;
434     htmlfile << "<h3 style=\"text-align:left; color:DD6600;\">"+title+"</h3>" << endl;
435    
436 dkralph 1.12 htmlfile << "</hr>" << endl;
437     htmlfile << "</hr>" << endl;
438     htmlfile << "Sample composition plots: " << endl;
439    
440 dkralph 1.1 htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;
441     htmlfile << "<tr>" << endl;
442 dkralph 1.12 htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/mtfitW.png\"><img src=\"plots/mtfitW.png\" alt=\"plots/mtfitW.png\" width=\"100%\"></a></td>" << endl;
443     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/both_met.png\"><img src=\"plots/both_met.png\" alt=\"plots/both_met.png\" width=\"100%\"></a></td>" << endl;
444     htmlfile << "<td width=\"25%\"><a></a></td>" << endl;
445     htmlfile << "<td width=\"25%\"><a></a></td>" << endl;
446 dkralph 1.1 htmlfile << "</tr>" << endl;
447     htmlfile << "</table>" << endl;
448    
449     htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;
450     htmlfile << "<tr>" << endl;
451 dkralph 1.12 htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/both_w_pt.png\"><img src=\"plots/both_w_pt.png\" alt=\"plots/both_w_pt.png\" width=\"100%\"></a></td>" << endl;
452     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/both_w_eta.png\"><img src=\"plots/both_w_eta.png\" alt=\"plots/both_w_eta.png\" width=\"100%\"></a></td>" << endl;
453     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;
454     htmlfile << "<td width=\"25%\"><a></a></td>" << endl;
455 dkralph 1.3 htmlfile << "</tr>" << endl;
456     htmlfile << "</table>" << endl;
457    
458     htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;
459     htmlfile << "<tr>" << endl;
460 dkralph 1.12 htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/both_l2_pt.png\"><img src=\"plots/both_l2_pt.png\" alt=\"plots/both_l2_pt.png\" width=\"100%\"></a></td>" << endl;
461     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/both_l2_eta.png\"><img src=\"plots/both_l2_eta.png\" alt=\"plots/both_l2_eta.png\" width=\"100%\"></a></td>" << endl;
462     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/both_l2_isoPF04.png\"><img src=\"plots/both_l2_isoPF04.png\" alt=\"plots/both_l2_isoPF04.png\" width=\"100%\"></a></td>" << endl;
463     htmlfile << "<td width=\"25%\"><a></a></td>" << endl;
464     htmlfile << "</tr>" << endl;
465     htmlfile << "</table>" << endl;
466    
467     htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;
468     htmlfile << "<tr>" << endl;
469     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/both_w_ip3ds_lo.png\"><img src=\"plots/both_w_ip3ds_lo.png\" alt=\"plots/both_w_ip3ds_lo.png\" width=\"100%\"></a></td>" << endl;
470     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/both_w_ip3ds.png\"><img src=\"plots/both_w_ip3ds.png\" alt=\"plots/both_w_ip3ds.png\" width=\"100%\"></a></td>" << endl;
471     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/both_w_d0.png\"><img src=\"plots/both_w_d0.png\" alt=\"plots/both_w_d0.png\" width=\"100%\"></a></td>" << endl;
472     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/both_w_dz.png\"><img src=\"plots/both_w_dz.png\" alt=\"plots/both_w_dz.png\" width=\"100%\"></a></td>" << endl;
473     htmlfile << "</tr>" << endl;
474     htmlfile << "</table>" << endl;
475    
476     htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;
477     htmlfile << "<tr>" << endl;
478     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/both_l2_ip3ds_lo.png\"><img src=\"plots/both_l2_ip3ds_lo.png\" alt=\"plots/both_l2_ip3ds_lo.png\" width=\"100%\"></a></td>" << endl;
479     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/both_l2_ip3ds.png\"><img src=\"plots/both_l2_ip3ds.png\" alt=\"plots/both_l2_ip3ds.png\" width=\"100%\"></a></td>" << endl;
480     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/both_l2_d0.png\"><img src=\"plots/both_l2_d0.png\" alt=\"plots/both_l2_d0.png\" width=\"100%\"></a></td>" << endl;
481     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/both_l2_dz.png\"><img src=\"plots/both_l2_dz.png\" alt=\"plots/both_l2_dz.png\" width=\"100%\"></a></td>" << endl;
482     htmlfile << "</tr>" << endl;
483     htmlfile << "</table>" << endl;
484    
485     htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;
486     htmlfile << "<tr>" << endl;
487     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/fr.png\"><img src=\"plots/fr.png\" alt=\"plots/fr.png\" width=\"100%\"></a></td>" << endl;
488     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/both_npv.png\"><img src=\"plots/both_npv.png\" alt=\"plots/both_npv.png\" width=\"100%\"></a></td>" << endl;
489     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;
490     htmlfile << "<td width=\"25%\"><a></a></td>" << endl;
491     htmlfile << "</tr>" << endl;
492     htmlfile << "</table>" << endl;
493    
494     htmlfile << "</hr>" << endl;
495     htmlfile << "</hr>" << endl;
496     htmlfile << "Fake rate closure plots: " << endl;
497    
498     htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;
499     htmlfile << "<tr>" << endl;
500 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;
501     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;
502     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;
503     htmlfile << "</tr>" << endl;
504     htmlfile << "</table>" << endl;
505    
506     htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;
507     htmlfile << "<tr>" << endl;
508 dkralph 1.12 htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/w_ip3ds_lo.png\"><img src=\"plots/w_ip3ds_lo.png\" alt=\"plots/w_ip3ds_lo.png\" width=\"100%\"></a></td>" << endl;
509     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;
510     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/w_d0.png\"><img src=\"plots/w_d0.png\" alt=\"plots/w_d0.png\" width=\"100%\"></a></td>" << endl;
511     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/w_dz.png\"><img src=\"plots/w_dz.png\" alt=\"plots/w_dz.png\" width=\"100%\"></a></td>" << endl;
512     htmlfile << "</tr>" << endl;
513     htmlfile << "</table>" << endl;
514    
515     htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;
516     htmlfile << "<tr>" << endl;
517     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/l2_ip3ds_lo.png\"><img src=\"plots/l2_ip3ds_lo.png\" alt=\"plots/l2_ip3ds_lo.png\" width=\"100%\"></a></td>" << endl;
518     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;
519     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/l2_d0.png\"><img src=\"plots/l2_d0.png\" alt=\"plots/l2_d0.png\" width=\"100%\"></a></td>" << endl;
520     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/l2_dz.png\"><img src=\"plots/l2_dz.png\" alt=\"plots/l2_dz.png\" width=\"100%\"></a></td>" << endl;
521     htmlfile << "</tr>" << endl;
522     htmlfile << "</table>" << endl;
523    
524     htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;
525     htmlfile << "<tr>" << endl;
526 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;
527     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;
528     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;
529 dkralph 1.12 htmlfile << "<td width=\"25%\"><a></a></td>" << endl;
530 dkralph 1.1 htmlfile << "</tr>" << endl;
531     htmlfile << "</table>" << endl;
532    
533     htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;
534     htmlfile << "<tr>" << endl;
535 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;
536     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;
537     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;
538 dkralph 1.12 htmlfile << "<td width=\"25%\"><a></a></td>" << endl;
539 dkralph 1.1 htmlfile << "</tr>" << endl;
540     htmlfile << "</table>" << endl;
541     htmlfile << "<hr />" << endl;
542    
543 dkralph 1.6 // htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;
544     // htmlfile << "<tr>" << endl;
545     // 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;
546     // htmlfile << "</tr>" << endl;
547     // htmlfile << "</table>" << endl;
548 dkralph 1.1
549     htmlfile << "<hr />" << endl;
550    
551     htmlfile << "<hr />" << endl;
552    
553     htmlfile << "</body>" << endl;
554     htmlfile << "</html>" << endl;
555     htmlfile.close();
556     }
557     //----------------------------------------------------------------------------------------
558 dkralph 1.9 RooFitResult *fitMt(FOFlags &ctrl, TH1D *hmcZ, TH1D *hmcW, TH1D *hdata, double &total_mc_scale, TString plotLabel)
559 dkralph 1.1 {
560 dkralph 1.2 double original_w_integral = integrateHist(hmcW);
561 dkralph 1.1
562 dkralph 1.2 RooRealVar mt("mt","MT [GeV]",hdata->GetXaxis()->GetXmin(),hdata->GetXaxis()->GetXmax());
563     RooBinning mtbins(hdata->GetNbinsX(),hdata->GetXaxis()->GetXmin(),hdata->GetXaxis()->GetXmax());
564     mt.setBinning(mtbins);
565 dkralph 1.1
566     // data
567 dkralph 1.2 RooDataHist obsMt("obsMt","obsMt",RooArgList(mt),hdata);
568 dkralph 1.1 // W mc pdf
569     hmcW->Smooth(2);
570 dkralph 1.2 RooDataHist hW("hW","hW",RooArgList(mt),hmcW);
571     RooHistPdf pW("pW","pW",RooArgSet(mt),hW,3);
572     // Z mc pdf
573     RooDataHist hZ("hZ","hZ",RooArgList(mt),hmcZ);
574     RooHistPdf pZ("pZ","pZ",RooArgSet(mt),hZ,3);
575 dkralph 1.1 // rayleigh pdf
576 dkralph 1.2 RooRealVar sig("sig","sig",400,15,585);
577     RooRealVar a1("a1","a1",-3,-10,10);
578     // RooGenericPdf rayl("rayl","rayl","mt*exp(-(mt)*(mt)/(2*sig*sig))",RooArgList(mt,sig));
579     RooGenericPdf rayl("rayl","rayl","mt*exp(-(mt)*(mt)/(sig*sig +2*sig*a1*(mt) + a1*a1*(mt)*(mt)))",RooArgList(mt,sig,a1));
580    
581     RooRealVar nObs("nObs","nObs",integrateHist(hdata));
582     RooRealVar nWandZ("nWandZ","nWandZ",0.98*nObs.getVal(),.95*nObs.getVal(),1.1*nObs.getVal());
583     RooRealVar zFrac("zFrac","zFrac",integrateHist(hmcZ) / (integrateHist(hmcZ) + integrateHist(hmcW)));
584     RooAddPdf modelWandZ("modelWandZ","modelWandZ",RooArgList(pZ,pW),zFrac);
585     RooRealVar nRayl("nRayl","nRayl",0.05*nObs.getVal(),0.0*nObs.getVal(),.05*nObs.getVal());
586    
587     RooExtendPdf extWandZ("extWandZ","extWandZ",modelWandZ,nWandZ);
588     RooExtendPdf extRayl("extRayl","extRayl",rayl,nRayl);
589    
590     RooAddPdf *model = new RooAddPdf("model","model",RooArgList(extWandZ,extRayl));
591    
592     RooFitResult *fitres = model->fitTo(obsMt,Extended(kTRUE),Save(kTRUE));
593     fitres->Print();
594    
595     TH1D *rayl_hist = (TH1D*)rayl.createHistogram("rayl_mt",mt,Binning(mtbins));
596     rayl_hist->Scale(nRayl.getVal()/integrateHist(rayl_hist));
597     TH1D *totalHist = (TH1D*)model->createHistogram("total_mt",mt,Binning(mtbins));
598     totalHist->Scale( (nWandZ.getVal() + nRayl.getVal()) / integrateHist(totalHist) );
599    
600     // // why the fuck do I need this?
601     // double tot = nWandZ.getVal() + nRayl.getVal();
602     // if(use_exp) tot += nExp.getVal();
603     // double hack_fac = integrateHist(hdata)/tot;
604     // rayl_hist->Scale(hack_fac);
605     // exp_hist->Scale(hack_fac);
606     // hmcZ->Scale(hack_fac);
607     // hmcW->Scale(hack_fac);
608    
609     // RooPlot *frame = mt.frame();
610     // obsMt.plotOn(frame);
611     // model->plotOn(frame);
612     // model->plotOn(frame,Components("pW"),LineStyle(kDashed),LineColor(kRed));
613     // model->plotOn(frame,Components("pZ"),LineColor(kRed));
614     // model->plotOn(frame,Components("rayl"),LineColor(kGreen));
615     // frame->Draw();
616     // can->SaveAs("~/public_html/foo.png");
617    
618     hmcW->Scale( (nWandZ.getVal()*(1-zFrac.getVal())) / integrateHist(hmcW) );
619     hmcZ->Scale( (nWandZ.getVal()*zFrac.getVal()) / integrateHist(hmcZ) );
620 dkralph 1.1
621 dkralph 1.9 CPlot plotMtFit("mtfitW","","MT [GeV]","Events",ctrl.outdir+"/"+plotLabel+"/"+ctrl.faketype+"/plots");
622 dkralph 1.2 plotMtFit.AddHist1D(hdata,"tight+fk'ble: " + integral_str(hdata),"E");
623     plotMtFit.AddHist1D(totalHist,"total MC: " + integral_str(totalHist),"hist",kRed);
624 dkralph 1.12 plotMtFit.AddToStack(rayl_hist,"jet: " + integral_str(rayl_hist),kBlue);
625 dkralph 1.2 plotMtFit.AddToStack(hmcZ,"Z+j MC: " + integral_str(hmcZ),kGray);
626     plotMtFit.AddToStack(hmcW,"W+j MC: " + integral_str(hmcW),kCyan-6);
627     plotMtFit.Draw(can,kTRUE,"png");
628 dkralph 1.1
629 dkralph 1.2 total_mc_scale = integrateHist(hmcW) / original_w_integral;
630 dkralph 1.1
631     return fitres;
632     }
633 dkralph 1.6 //----------------------------------------------------------------------------------------
634     void fillHist(CSample *cs, TString type, TString var, double val, double wgt, double wgt_lo, double wgt_hi)
635     {
636     (*(cs->hists)[type])[var]->Fill( val, wgt);
637     if(type=="pred") {
638     (*(cs->hists)[type+"_lo"])[var]->Fill( val, wgt_lo);
639     (*(cs->hists)[type+"_hi"])[var]->Fill( val, wgt_hi);
640     }
641     }