ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/NonMCBackground/src/compute_fakes.cc
Revision: 1.5
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.4: +360 -286 lines
Log Message:
*** empty log message ***

File Contents

# User Rev Content
1 dkralph 1.5 #ifndef CMSSW_BASE
2     #define CMSSW_BASE "../../"
3     #endif
4 dkralph 1.1 #include "compute_fakes.h"
5    
6     //=== MAIN =================================================================================================
7 dkralph 1.2 void makeHTML(FOFlags &ctrl, const TString outDir);
8     TGraphAsymmErrors* computeFakeRate1D(const TH1D* hpass, const TH1D* htotal);
9 dkralph 1.1 void computeFakeRate2D(const TH2D *hpass, const TH2D* htotal,
10     TH2D *hresult, TH2D* herrl, TH2D* herrh);
11 dkralph 1.5 map<TString,map<TString,TH1D*>* > init_hists(FOFlags &ctrl, TString str);
12     map<TString,map<TString,TH2D*>* > init_hists2d(FOFlags &ctrl, TString str);
13    
14     Double_t *ptbinning,*etabinning,*phibinning,*npvbinning,*massbinning;
15     Int_t nbinsPt,nbinsEta,nbinsPhi,nbinsNpv,nbinsMass;
16 dkralph 1.1
17     int main(int argc, char** argv)
18     {
19    
20 dkralph 1.5 double lumiVal = 4000;
21     // UInt_t minRun=999999999,maxRun=0;
22    
23 dkralph 1.1 TString format("png");
24     Bool_t doAbsEta = kFALSE;
25 dkralph 1.5 TString seletype("");//"ID mva + Iso mva");
26     TString tightness("");//"loose");
27 dkralph 1.1
28     //
29     // args
30     //--------------------------------------------------------------------------------------------------------------
31     FOFlags ctrl;
32     parse_foargs( argc, argv, ctrl );
33     ctrl.dump();
34 dkralph 1.5 assert(ctrl.faketype.Contains("mu_") || ctrl.faketype.Contains("ele_"));
35 dkralph 1.1
36     //--------------------------------------------------------------------------------------------------------------
37     // Settings
38     //==============================================================================================================
39    
40 dkralph 1.2 ifstream ifs;
41     string line;
42    
43 dkralph 1.5 TString cmsswpath(CMSSW_BASE + TString("/src"));
44     string xspath = (string(cmsswpath)+"/MitPhysics/data/xs.dat");
45     SimpleTable xstab(xspath.c_str());
46 dkralph 1.2
47 dkralph 1.5 // bin edges for kinematic variables
48     vector<Double_t> ptBinEdgesv,etaBinEdgesv,phiBinEdgesv,npvBinEdgesv,massBinEdgesv;
49 dkralph 1.2
50     // parse FR config info
51 dkralph 1.3 ifs.open(ctrl.inputfile.c_str());
52 dkralph 1.1 assert(ifs.is_open());
53     Int_t state=0;
54     while(getline(ifs,line)) {
55     if(line[0]=='#') continue;
56     if(line[0]=='%') {
57     state++;
58     continue;
59     }
60    
61 dkralph 1.5 Double_t edge;
62     stringstream ss(line);
63     ss >> edge;
64     if(state==0) { etaBinEdgesv.push_back(edge); }
65     else if(state==1) { ptBinEdgesv.push_back(edge); }
66     else if(state==2) { phiBinEdgesv.push_back(edge); }
67     else if(state==3) { npvBinEdgesv.push_back(edge); }
68     else if(state==4) { massBinEdgesv.push_back(edge); }
69     else assert(0);
70 dkralph 1.1 }
71     ifs.close();
72    
73 dkralph 1.5 nbinsPt = ptBinEdgesv.size()-1;
74     ptbinning = new Double_t[ptBinEdgesv.size()];
75 dkralph 1.1 for(UInt_t i=0; i<ptBinEdgesv.size(); i++) { ptbinning[i] = ptBinEdgesv[i]; }
76 dkralph 1.5
77     nbinsEta = etaBinEdgesv.size()-1;
78     etabinning = new Double_t[etaBinEdgesv.size()];
79 dkralph 1.1 for(UInt_t i=0; i<etaBinEdgesv.size(); i++) { etabinning[i] = etaBinEdgesv[i]; }
80 dkralph 1.5
81     nbinsPhi = phiBinEdgesv.size()-1;
82     phibinning = new Double_t[phiBinEdgesv.size()];
83 dkralph 1.1 for(UInt_t i=0; i<phiBinEdgesv.size(); i++) { phibinning[i] = phiBinEdgesv[i]; }
84 dkralph 1.2
85 dkralph 1.5 nbinsNpv = npvBinEdgesv.size()-1;
86     npvbinning = new Double_t[npvBinEdgesv.size()];
87     for(UInt_t i=0; i<npvBinEdgesv.size(); i++) { npvbinning[i] = npvBinEdgesv[i]; }
88    
89     nbinsMass = massBinEdgesv.size()-1;
90     massbinning = new Double_t[massBinEdgesv.size()];
91     for(UInt_t i=0; i<massBinEdgesv.size(); i++) { massbinning[i] = massBinEdgesv[i]; }
92    
93     // parse input information
94     vector<CSample*> samplev;
95     TString ntupledir,inputLabel;
96     ifs.open(ctrl.config); assert(ifs.is_open());
97     while(getline(ifs,line)) {
98     if(line[0]=='#') continue;
99 dkralph 1.1
100 dkralph 1.5 stringstream ss(line);
101     if(line[0]=='^') {
102     TString dummy;
103     if(TString(line).Contains("^ntupdir")) ss >> dummy >> ntupledir;
104     if(TString(line).Contains("^label")) ss >> dummy >> inputLabel;
105     // if(TString(line).Contains("^skim")) ss >> dummy >> ntupledir;
106     continue;
107 dkralph 1.1 }
108    
109 dkralph 1.5 if(line[0]=='$') {
110     samplev.push_back(new CSample());
111     stringstream ss(line);
112     string chr;
113     TString sname;
114     Int_t color;
115     ss >> chr >> sname >> color;
116     string legend = line.substr(line.find('@')+1);
117     samplev.back()->name = sname;
118     samplev.back()->legend = legend;
119     samplev.back()->color = color;
120     samplev.back()->hists = init_hists(ctrl,sname);
121     samplev.back()->hists2d = init_hists2d(ctrl,sname);
122     continue;
123 dkralph 1.1 }
124    
125 dkralph 1.5 TString dataset;
126     bool isdata;
127     double xsec;
128     ss >> dataset >> isdata >> xsec;
129     assert(ctrl.muSele==ctrl.eleSele);
130     TString fname = ntupledir+"/"+inputLabel+"/"+dataset+"/merged.root";
131     samplev.back()->isdata = isdata;
132     (samplev.back()->fsv).push_back(new filestuff(dataset,fname,dataset,isdata,false));
133 dkralph 1.1 }
134 dkralph 1.5 ifs.close();
135 dkralph 1.1
136 dkralph 1.5 //--------------------------------------------------------------------------------------------------------------
137     // Main analysis code
138     //==============================================================================================================
139 dkralph 1.1
140 dkralph 1.5 for(unsigned isam=0; isam<samplev.size(); isam++) {
141     CSample *cs = samplev[isam];
142     cout << cs->name << endl;
143     map<TString,map<TString,TH1D*>* > *hists = &(cs->hists);
144     map<TString,map<TString,TH2D*>* > *hists2d = &(cs->hists2d);
145     for(unsigned ifs=0; ifs<(cs->fsv).size(); ifs++) {
146     filestuff *fs = (cs->fsv)[ifs];
147     cout << "\t" << fs->fname_ << endl;
148     TFile *infile = TFile::Open(fs->fname_); assert(infile->IsOpen());
149     TTree *intree = (TTree*)infile->Get("FO");
150     foinfo foi;
151     intree->SetBranchAddress("FO",&foi);
152    
153     for(unsigned ientry=0; ientry<intree->GetEntries(); ientry++) {
154     intree->GetEntry(ientry);
155     // if(fs->isdata_) {
156     // if(foi.run < minRun) minRun = foi.run;
157     // if(foi.run > maxRun) maxRun = foi.run;
158     // }
159    
160     assert(abs(foi.type)==11 || abs(foi.type)==13);
161     if(abs(foi.type)==13 && ctrl.faketype.Contains("ele_")) continue;
162     if(abs(foi.type)==11 && ctrl.faketype.Contains("mu_")) continue;
163    
164     double etaVal = doAbsEta ? fabs(foi.foeta) : foi.foeta;
165     double ptVal = foi.fopt;//fs->isdata_ ? foi.pt : 0.99*foi.pt;
166     double massVal = foi.mass;//fs->isdata_ ? foi.mass : 0.993*foi.mass;
167     if(ptVal < ptBinEdgesv.front() || ptVal > ptBinEdgesv.back()) continue;
168     if(etaVal < etaBinEdgesv.front() || etaVal > etaBinEdgesv.back()) continue;
169     if(foi.npv < npvBinEdgesv.front() || foi.npv > npvBinEdgesv.back()) continue;
170     if(massVal < massBinEdgesv.front() || massVal > massBinEdgesv.back()) continue;
171    
172     if(foi.met>25) continue;
173     if(massVal < 80) continue;
174    
175     double wgt=1;
176     if(!fs->isdata_) {
177     double xsWgt = lumiVal*xstab.Get(fs->dataset_)/fs->total_entries_;
178     double puWgt = weightTrue2012(foi.npu);
179     wgt = xsWgt*puWgt;
180     }
181    
182     if(foi.l1pt < 20 && foi.l2pt < 20) continue;
183     if(foi.l1pt < 10 || foi.l2pt < 10) continue;
184    
185    
186     (*(*hists)["denom"])["run"]->Fill( foi.run, wgt);
187     (*(*hists)["denom"])["met"]->Fill( foi.met, wgt);
188     (*(*hists)["denom"])["mass"]->Fill( massVal, wgt);
189     (*(*hists)["denom"])["l1pt"]->Fill( foi.l1pt, wgt);
190     (*(*hists)["denom"])["l2pt"]->Fill( foi.l2pt, wgt);
191     (*(*hists)["denom"])["npv"]->Fill( foi.npv, wgt);
192     (*(*hists)["denom"])["Pt"]->Fill( ptVal, wgt);
193     (*(*hists)["denom"])["Eta"]->Fill( etaVal, wgt);
194     (*(*hists)["denom"])["foPt"]->Fill( ptVal, wgt);
195     (*(*hists)["denom"])["foEta"]->Fill( etaVal, wgt);
196     (*(*hists)["denom"])["foPhi"]->Fill( foi.fophi, wgt);
197     (*(*hists)["denom"])["foNpv"]->Fill( foi.npv, wgt);
198     (*(*hists)["denom"])["foMass"]->Fill( massVal, wgt);
199     (*(*hists2d)["denom"])["frEtaPt"]->Fill( etaVal, ptVal, wgt);
200    
201     if(foi.fopass) {
202     (*(*hists)["numer"])["run"]->Fill( foi.run, wgt);
203     (*(*hists)["numer"])["mass"]->Fill( massVal, wgt);
204     (*(*hists)["numer"])["met"]->Fill( foi.met, wgt);
205     (*(*hists)["numer"])["l1pt"]->Fill( foi.l1pt, wgt);
206     (*(*hists)["numer"])["l2pt"]->Fill( foi.l2pt, wgt);
207     (*(*hists)["numer"])["npv"]->Fill( foi.npv, wgt);
208     (*(*hists)["numer"])["Pt"]->Fill( ptVal, wgt);
209     (*(*hists)["numer"])["Eta"]->Fill( etaVal, wgt);
210     (*(*hists)["numer"])["foPt"]->Fill( ptVal, wgt);
211     (*(*hists)["numer"])["foEta"]->Fill( etaVal, wgt);
212     (*(*hists)["numer"])["foPhi"]->Fill( foi.fophi, wgt);
213     (*(*hists)["numer"])["foNpv"]->Fill( foi.npv, wgt);
214     (*(*hists)["numer"])["foMass"]->Fill( massVal, wgt);
215     (*(*hists2d)["numer"])["frEtaPt"]->Fill( etaVal, ptVal, wgt);
216     }
217     }
218     }
219    
220     map<TString,TGraphAsymmErrors*> *graphs = &(cs->graphs);
221     (*graphs)["frPt"] = computeFakeRate1D((*(*hists)["numer"])["foPt"],(*(*hists)["denom"])["foPt"]);
222     (*graphs)["frEta"] = computeFakeRate1D((*(*hists)["numer"])["foEta"],(*(*hists)["denom"])["foEta"]);
223     (*graphs)["frPhi"] = computeFakeRate1D((*(*hists)["numer"])["foPhi"],(*(*hists)["denom"])["foPhi"]);
224     (*graphs)["frNpv"] = computeFakeRate1D((*(*hists)["numer"])["foNpv"],(*(*hists)["denom"])["foNpv"]);
225     (*graphs)["frMass"] = computeFakeRate1D((*(*hists)["numer"])["foMass"],(*(*hists)["denom"])["foMass"]);
226    
227     computeFakeRate2D((*(*hists2d)["numer"])["frEtaPt"],(*(*hists2d)["denom"])["frEtaPt"],(*(*hists2d)["eff"])["frEtaPt"],(*(*hists2d)["errl"])["frEtaPt"],(*(*hists2d)["errh"])["frEtaPt"]);
228    
229     (*(*hists2d)["eff"])["frEtaPt"]->SetName("frEtaPt");
230     (*(*hists2d)["errl"])["frEtaPt"]->SetName("errlEtaPt");
231     (*(*hists2d)["errh"])["frEtaPt"]->SetName("errhEtaPt");
232     }
233     // cout << setw(12) << minRun << setw(12) << maxRun << endl;
234 dkralph 1.1
235     //--------------------------------------------------------------------------------------------------------------
236     // Make plots
237     //==============================================================================================================
238     TCanvas *c = new TCanvas("c","c",800,600);
239     char ylabel[100];
240     TString etalabel = (doAbsEta) ? "#scale[1.3]{|#eta|}" : "#scale[1.3]{#eta}";
241    
242 dkralph 1.5 TFile runHistFile(ctrl.outdir+"/runs.root","recreate");
243     assert(samplev.size()>0);
244     TH1D *hTmpData = (*(samplev[0]->hists)["denom"])["mass"];
245     TH1D *hTmpMc = (*(samplev[1]->hists)["denom"])["mass"];
246     double tmpMcScale = integrateHist(hTmpData)/integrateHist(hTmpMc);
247     map<TString,TH1D*>::iterator it_v;
248     for(it_v=(*(samplev[0]->hists)["numer"]).begin(); it_v!=(*(samplev[0]->hists)["numer"]).end(); it_v++) {
249     TString var((*it_v).first);
250     CPlot plotDenom("denom"+var,"","","Events",ctrl.outdir+"/plots");
251     CPlot plotNumer("numer"+var,"","","Events",ctrl.outdir+"/plots");
252     for(unsigned isam=0; isam<samplev.size(); isam++) {
253     CSample *cs = samplev[isam];
254     TH1D *hNumer = (*(cs->hists)["numer"])[var];
255     TH1D *hDenom = (*(cs->hists)["denom"])[var];
256     if(cs->isdata && var=="run") {
257     assert(hNumer->GetBinContent(0)==0 && hNumer->GetBinContent(hNumer->GetXaxis()->GetNbins()+1)==0);
258     assert(hDenom->GetBinContent(0)==0 && hDenom->GetBinContent(hDenom->GetXaxis()->GetNbins()+1)==0);
259     hNumer->Write("run_1");
260     hDenom->Write("run_2");
261     }
262     plotDenom.SetXTitle(TString("denominator ")+hDenom->GetXaxis()->GetTitle());
263     plotNumer.SetXTitle(TString("numerator ")+hNumer->GetXaxis()->GetTitle());
264     if(cs->isdata) {
265     plotDenom.AddHist1D(hDenom,cs->legend,"E");
266     plotNumer.AddHist1D(hNumer,cs->legend,"E");
267     } else if(var != "run") {
268     hDenom->Scale(tmpMcScale);
269     hNumer->Scale(tmpMcScale);
270     plotDenom.AddToStack(hDenom,cs->legend,cs->color);
271     plotNumer.AddToStack(hNumer,cs->legend,cs->color);
272     }
273     }
274     plotDenom.TransLegend(0.1,0);
275     plotNumer.TransLegend(0.1,0);
276     if(var.Contains("eta",TString::kIgnoreCase)) {
277     plotDenom.SetYRange(0,1.4*(plotDenom.GetStack()->GetMaximum()));
278     plotNumer.SetYRange(0,1.4*(plotNumer.GetStack()->GetMaximum()));
279     }
280     plotDenom.Draw(c,kTRUE,format);
281     plotNumer.Draw(c,kTRUE,format);
282 dkralph 1.1 }
283 dkralph 1.5 runHistFile.Close();
284    
285     map<TString,TGraphAsymmErrors*>::iterator it_gr;
286     for(it_gr=samplev[0]->graphs.begin(); it_gr!=samplev[0]->graphs.end(); it_gr++) {
287     TString var((*it_gr).first);
288     CPlot plotGr(var,"","","#scale[1.5]{#varepsilon_{fake}}",ctrl.outdir+"/plots");
289     for(unsigned isam=0; isam<samplev.size(); isam++) {
290     CSample *cs = samplev[isam];
291     TGraphAsymmErrors *gr = (cs->graphs)[var];
292     plotGr.SetXTitle(var);
293     plotGr.AddGraph(gr,cs->legend,"",cs->color);
294     double thisMax = TMath::MaxElement(gr->GetN(),gr->GetY());
295     if(1.4*thisMax > plotGr.GetYMax())
296     plotGr.SetYRange(0,1.4*thisMax);
297     // plotFRPt.AddTextBox(seletype+" "+tightness,0.25,0.75,0.45,0.85);
298     }
299     plotGr.Draw(c,kTRUE,format);
300 dkralph 1.1 }
301    
302     // gStyle->SetPalette(1);
303     c->SetRightMargin(0.15);
304 dkralph 1.5 c->SetLeftMargin(0.15);
305     TH2D *frEtaPt = (*(samplev[0]->hists2d)["eff"])["frEtaPt"];
306 dkralph 1.1 frEtaPt->SetTitleOffset(1.2,"Y");
307     // TPaletteAxis *paxis = (TPaletteAxis*)frEtaPt->GetListOfFunctions()->FindObject("palette");
308     // paxis->SetX1NDC(0.87);
309     // paxis->SetX2NDC(0.92);
310 dkralph 1.5 CPlot plotFRPtEta("frPtEta","",etalabel,"p_{T} [GeV/c]",ctrl.outdir+"/plots");
311 dkralph 1.1 plotFRPtEta.AddHist2D(frEtaPt,"COLZ");
312     plotFRPtEta.Draw(c,kTRUE,format);
313    
314     //--------------------------------------------------------------------------------------------------------------
315     // Summary print out
316     //==============================================================================================================
317 dkralph 1.5 for(unsigned isam=0; isam<samplev.size(); isam++) {
318     CSample *cs = samplev[isam];
319     TString outfilename(ctrl.outdir + TString("/fr.root"));
320     if(!cs->isdata) outfilename.ReplaceAll("fr.root","fr-"+cs->name+".root");
321     TFile outfile(outfilename,"RECREATE");
322     (cs->graphs)["frPt"]->Write();
323     (cs->graphs)["frEta"]->Write();
324     (*(cs->hists2d)["eff"])["frEtaPt"]->Write();
325     (*(cs->hists2d)["errl"])["frEtaPt"]->Write();
326     (*(cs->hists2d)["errh"])["frEtaPt"]->Write();
327     }
328    
329 dkralph 1.2 makeHTML(ctrl,ctrl.outdir);
330 dkralph 1.1
331 dkralph 1.2 cout << " <> Output saved in " << ctrl.outdir << "/" << endl;
332 dkralph 1.1 cout << endl;
333 dkralph 1.5
334     delete [] ptbinning;
335     delete [] etabinning;
336     delete [] phibinning;
337     delete [] npvbinning;
338 dkralph 1.1 }
339    
340    
341     //=== FUNCTION DEFINITIONS ======================================================================================
342    
343 dkralph 1.5 //----------------------------------------------------------------------------------------
344     map<TString,map<TString,TH1D*>* > init_hists(FOFlags &ctrl, TString str)
345     {
346     map<TString,map<TString,TH1D*>* > hists;
347     map<TString,TH1D*> *h_denom = new map<TString,TH1D*>; hists["denom"] = h_denom;
348     map<TString,TH1D*> *h_numer = new map<TString,TH1D*>; hists["numer"] = h_numer;
349     map<TString,map<TString,TH1D*>* >::iterator it_h;
350    
351     UInt_t minRun=160440,maxRun=196535;
352     int otherNmassBins = 20;
353     Double_t otherMassBinning[] = { 55,60,65,70,75,80,84,85.5,87,88.5,90,91.5,93,94.5,96,99,105,110,115,120,125 };
354     double ptMax,ptMin=4;
355     if(ctrl.faketype.Contains("mu_"))
356     ptMax = 15;
357     else
358     ptMax = 40;
359     for(it_h=hists.begin(); it_h!=hists.end(); it_h++) {
360     (*((*it_h).second))["run"] = new TH1D(TString("run") +"_"+(*it_h).first+str,";#bf{run};", 50,minRun,maxRun); (*((*it_h).second))["run"]->Sumw2();
361     (*((*it_h).second))["mass"] = new TH1D(TString("mass") +"_"+(*it_h).first+str,";#bf{m_{ll}};", otherNmassBins,otherMassBinning); (*((*it_h).second))["mass"]->Sumw2();
362     (*((*it_h).second))["l1pt"] = new TH1D(TString("l1pt") +"_"+(*it_h).first+str,";#bf{l1 p_{T} [GeV]};", 37,0,80); (*((*it_h).second))["l1pt"]->Sumw2();
363     (*((*it_h).second))["l2pt"] = new TH1D(TString("l2pt") +"_"+(*it_h).first+str,";#bf{l2 p_{T} [GeV]};", 37,0,80); (*((*it_h).second))["l2pt"]->Sumw2();
364     (*((*it_h).second))["npv"] = new TH1D(TString("npv") +"_"+(*it_h).first+str,";#bf{n PV};", 50,-.5,49.5); (*((*it_h).second))["npv"]->Sumw2();
365     (*((*it_h).second))["met"] = new TH1D(TString("met") +"_"+(*it_h).first+str,";#bf{MET};", 25,0,35); (*((*it_h).second))["met"]->Sumw2();
366     (*((*it_h).second))["Pt"] = new TH1D(TString("Pt") +"_"+(*it_h).first+str,";#bf{PT};", 25,ptMin,ptMax); (*((*it_h).second))["Pt"]->Sumw2();
367     (*((*it_h).second))["Eta"] = new TH1D(TString("Eta") +"_"+(*it_h).first+str,";#bf{ETA};", 15,0,2.5); (*((*it_h).second))["Eta"]->Sumw2();
368    
369     (*((*it_h).second))["foPt"] = new TH1D(TString("foPt") +"_"+(*it_h).first+str,";#bf{p_{T} [GeV]};", nbinsPt,ptbinning); (*((*it_h).second))["foPt"]->Sumw2();
370     (*((*it_h).second))["foEta"] = new TH1D(TString("foEta") +"_"+(*it_h).first+str,";#bf{#eta};", nbinsEta,etabinning); (*((*it_h).second))["foEta"]->Sumw2();
371     (*((*it_h).second))["foPhi"] = new TH1D(TString("foPhi") +"_"+(*it_h).first+str,";#bf{#phi};", nbinsPhi,phibinning); (*((*it_h).second))["foPhi"]->Sumw2();
372     (*((*it_h).second))["foNpv"] = new TH1D(TString("foNpv") +"_"+(*it_h).first+str,";#bf{NPV };", nbinsNpv,npvbinning); (*((*it_h).second))["foNpv"]->Sumw2();
373     (*((*it_h).second))["foMass"] = new TH1D(TString("foMass") +"_"+(*it_h).first+str,";#bf{m_{ll} };", nbinsMass,massbinning); (*((*it_h).second))["foMass"]->Sumw2();
374     }
375    
376     return hists;
377     }
378     //----------------------------------------------------------------------------------------
379     map<TString,map<TString,TH2D*>* > init_hists2d(FOFlags &ctrl, TString str)
380     {
381     map<TString,map<TString,TH2D*>* > hists;
382     map<TString,TH2D*> *h_denom = new map<TString,TH2D*>; hists["denom"] = h_denom;
383     map<TString,TH2D*> *h_numer = new map<TString,TH2D*>; hists["numer"] = h_numer;
384     map<TString,TH2D*> *h_eff = new map<TString,TH2D*>; hists["eff"] = h_eff;
385     map<TString,TH2D*> *h_errl = new map<TString,TH2D*>; hists["errl"] = h_errl;
386     map<TString,TH2D*> *h_errh = new map<TString,TH2D*>; hists["errh"] = h_errh;
387     map<TString,map<TString,TH2D*>* >::iterator it_h;
388    
389     for(it_h=hists.begin(); it_h!=hists.end(); it_h++) {
390     (*((*it_h).second))["frEtaPt"] = new TH2D(TString("frEtaPt") +"_"+(*it_h).first+str,";#bf{#eta};#bf{p_{T}}", nbinsEta,etabinning,nbinsPt,ptbinning); (*((*it_h).second))["frEtaPt"]->Sumw2();
391     }
392    
393     return hists;
394     }
395 dkralph 1.1 //--------------------------------------------------------------------------------------------------
396 dkralph 1.2 TGraphAsymmErrors* computeFakeRate1D(const TH1D* hpass, const TH1D* htotal)
397 dkralph 1.1 {
398     const Int_t nbins = htotal->GetNbinsX();
399     Double_t xval[nbins], xerr[nbins], yval[nbins], yerrl[nbins], yerrh[nbins];
400     for(Int_t ibin=1; ibin<=nbins; ibin++) {
401     xval[ibin-1] = htotal->GetBinCenter(ibin);
402     xerr[ibin-1] = 0.5*htotal->GetBinWidth(ibin);
403    
404     Int_t total = htotal->GetBinContent(ibin);
405     Int_t passed = hpass->GetBinContent(ibin);
406     yval[ibin-1] = (total>0) ? (Double_t)passed/(Double_t)total : 0;
407     yerrl[ibin-1] = (total>0) ? yval[ibin-1] - TEfficiency::ClopperPearson(total,passed,0.68269,kFALSE) : 0;
408     yerrh[ibin-1] = (total>0) ? TEfficiency::ClopperPearson(total,passed,0.68269,kTRUE) - yval[ibin-1] : 0;
409     }
410    
411     return new TGraphAsymmErrors(nbins,xval,yval,xerr,xerr,yerrl,yerrh);
412     }
413    
414     //--------------------------------------------------------------------------------------------------
415     void computeFakeRate2D(const TH2D *hpass, const TH2D* htotal,
416     TH2D *hresult, TH2D* herrl, TH2D* herrh)
417     {
418     assert(hresult);
419     assert(herrl);
420     assert(herrh);
421    
422     const Int_t nbinsx = htotal->GetNbinsX();
423     const Int_t nbinsy = htotal->GetNbinsY();
424     for(Int_t ix=1; ix<=nbinsx; ix++) {
425     for(Int_t iy=1; iy<=nbinsy; iy++) {
426     Int_t total = htotal->GetCellContent(ix,iy);
427     Int_t passed = hpass->GetCellContent(ix,iy);
428     Double_t eff = (total>0) ? (Double_t)passed/(Double_t)total : 0;
429     Double_t errl = (total>0) ? eff - TEfficiency::ClopperPearson(total,passed,0.68269,kFALSE) : 0;
430     Double_t errh = (total>0) ? TEfficiency::ClopperPearson(total,passed,0.68269,kTRUE) - eff : 0;
431     hresult->SetCellContent(ix,iy,eff);
432     herrl ->SetCellContent(ix,iy,errl);
433     herrh ->SetCellContent(ix,iy,errh);
434     }
435     }
436     }
437    
438     //--------------------------------------------------------------------------------------------------
439 dkralph 1.2 void makeHTML(FOFlags &ctrl, const TString outDir)
440 dkralph 1.1 {
441 dkralph 1.5 TString title;
442     if(ctrl.faketype.Contains("mu_")) title = "Muon fakes: " + ctrl.faketype;
443     else if(ctrl.faketype.Contains("ele_")) title = "Electron fakes: " + ctrl.faketype;
444     else assert(0);
445 dkralph 1.1 ofstream htmlfile;
446     char htmlfname[100];
447     sprintf(htmlfname,"%s/fr.html",outDir.Data());
448     htmlfile.open(htmlfname);
449 dkralph 1.2
450 dkralph 1.1 htmlfile << "<!DOCTYPE html" << endl;
451     htmlfile << " PUBLIC \"-//W3C//DTD HTML 3.2//EN\">" << endl;
452     htmlfile << "<html>" << endl;
453    
454 dkralph 1.2 htmlfile << "<head><title>"+title+"</title></head>" << endl;
455 dkralph 1.1 htmlfile << "<body bgcolor=\"EEEEEE\">" << endl;
456 dkralph 1.2 htmlfile << "<h3 style=\"text-align:left; color:DD6600;\">"+title+"</h3>" << endl;
457 dkralph 1.1
458     htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;
459     htmlfile << "<tr>" << endl;
460 dkralph 1.5 htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/frPtEta.png\"><img src=\"plots/frPtEta.png\" alt=\"plots/frPtEta.png\" width=\"100%\"></a></td>" << endl;
461     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/denomrun.png\"><img src=\"plots/denomrun.png\" alt=\"plots/denomrun.png\" width=\"100%\"></a></td>" << endl;
462     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/numerrun.png\"><img src=\"plots/numerrun.png\" alt=\"plots/numerrun.png\" width=\"100%\"></a></td>" << endl;
463     htmlfile << "<td width=\"25%\"><a></a></td>" << endl;
464     htmlfile << "<td width=\"25%\"><a></a></td>" << endl;
465     htmlfile << "<tr>" << endl;
466     htmlfile << "<td width=\"33.3%\"><a target=\"_blank\" href=\"plots/frMass.png\"><img src=\"plots/frMass.png\" alt=\"plots/frMass.png\" width=\"100%\"></a></td>" << endl;
467     htmlfile << "<td width=\"33.3%\"><a target=\"_blank\" href=\"plots/frPt.png\"><img src=\"plots/frPt.png\" alt=\"plots/frPt.png\" width=\"100%\"></a></td>" << endl;
468     htmlfile << "<td width=\"33.3%\"><a target=\"_blank\" href=\"plots/frEta.png\"><img src=\"plots/frEta.png\" alt=\"plots/frEta.png\" width=\"100%\"></a></td>" << endl;
469     htmlfile << "</tr>" << endl;
470     htmlfile << "<tr>" << endl;
471     htmlfile << "<td width=\"33.3%\"><a target=\"_blank\" href=\"plots/frPhi.png\"><img src=\"plots/frPhi.png\" alt=\"plots/frPhi.png\" width=\"100%\"></a></td>" << endl;
472     htmlfile << "<td width=\"33.3%\"><a target=\"_blank\" href=\"plots/frNpv.png\"><img src=\"plots/frNpv.png\" alt=\"plots/frNpv.png\" width=\"100%\"></a></td>" << endl;
473     htmlfile << "<td width=\"33.3%\"><a></a></td>" << endl;
474 dkralph 1.1 htmlfile << "</tr>" << endl;
475     htmlfile << "<tr>" << endl;
476 dkralph 1.2 // htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/frnpv.png\"><img src=\"plots/frnpv.png\" alt=\"plots/frnpv.png\" width=\"100%\"></a></td>" << endl;
477     // htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/frptnpv.png\"><img src=\"plots/frptnpv.png\" alt=\"plots/frptnpv.png\" width=\"100%\"></a></td>" << endl;
478     // htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/fretanpv.png\"><img src=\"plots/fretanpv.png\" alt=\"plots/fretanpv.png\" width=\"100%\"></a></td>" << endl;
479     // htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/frphinpv.png\"><img src=\"plots/frphinpv.png\" alt=\"plots/frphinpv.png\" width=\"100%\"></a></td>" << endl;
480 dkralph 1.1 htmlfile << "</tr>" << endl;
481     htmlfile << "</table>" << endl;
482     htmlfile << "<hr />" << endl;
483    
484     htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;
485 dkralph 1.5
486 dkralph 1.1 htmlfile << "<tr>" << endl;
487 dkralph 1.5 htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/denomPt.png\"><img src=\"plots/denomPt.png\" alt=\"plots/denomPt.png\" width=\"100%\"></a></td>" << endl;
488     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/denomEta.png\"><img src=\"plots/denomEta.png\" alt=\"plots/denomEta.png\" width=\"100%\"></a></td>" << endl;
489 dkralph 1.2 // htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/denomphi.png\"><img src=\"plots/denomphi.png\" alt=\"plots/denomphi.png\" width=\"100%\"></a></td>" << endl;
490 dkralph 1.5 htmlfile << "<td width=\"25%\"><a></a></td>" << endl;
491     htmlfile << "<td width=\"25%\"><a></a></td>" << endl;
492 dkralph 1.1 htmlfile << "</tr>" << endl;
493 dkralph 1.5
494 dkralph 1.1 htmlfile << "<tr>" << endl;
495 dkralph 1.5 htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/numerPt.png\"><img src=\"plots/numerPt.png\" alt=\"plots/numerPt.png\" width=\"100%\"></a></td>" << endl;
496     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/numerEta.png\"><img src=\"plots/numerEta.png\" alt=\"plots/numerEta.png\" width=\"100%\"></a></td>" << endl;
497     htmlfile << "<td width=\"25%\"><a></a></td>" << endl;
498     htmlfile << "<td width=\"25%\"><a></a></td>" << endl;
499 dkralph 1.2 // htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/numerphi.png\"><img src=\"plots/numerphi.png\" alt=\"plots/numerphi.png\" width=\"100%\"></a></td>" << endl;
500 dkralph 1.1 htmlfile << "</tr>" << endl;
501 dkralph 1.5
502     htmlfile << "<tr>" << endl;
503     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/denommass.png\"><img src=\"plots/denommass.png\" alt=\"plots/denommass.png\" width=\"100%\"></a></td>" << endl;
504     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/denoml1pt.png\"><img src=\"plots/denoml1pt.png\" alt=\"plots/denoml1pt.png\" width=\"100%\"></a></td>" << endl;
505     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/denoml2pt.png\"><img src=\"plots/denoml2pt.png\" alt=\"plots/denoml2pt.png\" width=\"100%\"></a></td>" << endl;
506     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/denomnpv.png\"><img src=\"plots/denomnpv.png\" alt=\"plots/denomnpv.png\" width=\"100%\"></a></td>" << endl;
507     htmlfile << "</tr>" << endl;
508    
509 dkralph 1.1 htmlfile << "<tr>" << endl;
510 dkralph 1.5 htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/numermass.png\"><img src=\"plots/numermass.png\" alt=\"plots/numermass.png\" width=\"100%\"></a></td>" << endl;
511     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/numerl1pt.png\"><img src=\"plots/numerl1pt.png\" alt=\"plots/numerl1pt.png\" width=\"100%\"></a></td>" << endl;
512     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/numerl2pt.png\"><img src=\"plots/numerl2pt.png\" alt=\"plots/numerl2pt.png\" width=\"100%\"></a></td>" << endl;
513     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/numernpv.png\"><img src=\"plots/numernpv.png\" alt=\"plots/numernpv.png\" width=\"100%\"></a></td>" << endl;
514 dkralph 1.1 htmlfile << "</tr>" << endl;
515     htmlfile << "</table>" << endl;
516     htmlfile << "<hr />" << endl;
517    
518 dkralph 1.5 htmlfile << "<tr>" << endl;
519     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/denommet.png\"><img src=\"plots/denommet.png\" alt=\"plots/denommet.png\" width=\"25%\"></a></td>" << endl;
520     htmlfile << "<td width=\"25%\"><a></a></td>" << endl;
521     htmlfile << "<td width=\"25%\"><a></a></td>" << endl;
522     htmlfile << "<td width=\"25%\"><a></a></td>" << endl;
523     htmlfile << "</tr>" << endl;
524    
525     htmlfile << "<tr>" << endl;
526     htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/numermet.png\"><img src=\"plots/numermet.png\" alt=\"plots/numermet.png\" width=\"25%\"></a></td>" << endl;
527     htmlfile << "<td width=\"25%\"><a></a></td>" << endl;
528     htmlfile << "<td width=\"25%\"><a></a></td>" << endl;
529     htmlfile << "<td width=\"25%\"><a></a></td>" << endl;
530     htmlfile << "</tr>" << endl;
531    
532 dkralph 1.1 htmlfile << "<table border=\"0\" cellspacing=\"5\" width=\"100%\">" << endl;
533 dkralph 1.2 // htmlfile << "<tr>" << endl;
534     // htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/jetpt.png\"><img src=\"plots/jetpt.png\" alt=\"plots/jetpt.png\" width=\"100%\"></a></td>" << endl;
535     // htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/jetptlog.png\"><img src=\"plots/jetptlog.png\" alt=\"plots/jetptlog.png\" width=\"100%\"></a></td>" << endl;
536     // htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/frjetpt.png\"><img src=\"plots/frjetpt.png\" alt=\"plots/frjetpt.png\" width=\"100%\"></a></td>" << endl;
537     // htmlfile << "<td width=\"25%\"></td>" << endl;
538     // htmlfile << "</tr>" << endl;
539     // htmlfile << "<tr>" << endl;
540     // htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/fojetpt.png\"><img src=\"plots/fojetpt.png\" alt=\"plots/fojetpt.png\" width=\"100%\"></a></td>" << endl;
541     // htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/fojetptlog.png\"><img src=\"plots/fojetptlog.png\" alt=\"plots/fojetptlog.png\" width=\"100%\"></a></td>" << endl;
542     // htmlfile << "<td width=\"25%\"><a target=\"_blank\" href=\"plots/frfojetpt.png\"><img src=\"plots/frfojetpt.png\" alt=\"plots/frfojetpt.png\" width=\"100%\"></a></td>" << endl;
543     // htmlfile << "<td width=\"25%\"></td>" << endl;
544     // htmlfile << "</tr>" << endl;
545 dkralph 1.1 htmlfile << "</table>" << endl;
546     htmlfile << "<hr />" << endl;
547    
548     htmlfile << "</body>" << endl;
549     htmlfile << "</html>" << endl;
550     htmlfile.close();
551     }