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