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