ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/AnalysisFramework/Plotting/Modules/Plotting_Functions.C
Revision: 1.12
Committed: Fri Jul 8 14:53:46 2011 UTC (13 years, 10 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.11: +39 -27 lines
Log Message:
Added several plots to understand sidebands better; made sideband definition global (now in Setup.C)

File Contents

# User Rev Content
1 buchmann 1.1 #include <iostream>
2     #include <vector>
3     #include <sys/stat.h>
4    
5     #include <TCut.h>
6     #include <TROOT.h>
7     #include <TCanvas.h>
8     #include <TMath.h>
9     #include <TColor.h>
10     #include <TPaveText.h>
11     #include <TRandom.h>
12     #include <TH1.h>
13     #include <TH2.h>
14     #include <TF1.h>
15     #include <TSQLResult.h>
16    
17 buchmann 1.9 //#include "TTbar_stuff.C"
18 buchmann 1.1 using namespace std;
19    
20     using namespace PlottingSetup;
21    
22     void find_peaks(float &MCPeak,float &MCPeakError, float &DataPeak, float &DataPeakError, float &MCSigma, float &DataSigma, stringstream &result)
23     {
24 buchmann 1.8 TH1F *rawJZBeemmMC = allsamples.Draw("rawJZBeemmMC",jzbvariablemc,100,-50,50, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets,mc, luminosity);
25     TH1F *rawJZBeemmData = allsamples.Draw("rawJZBeemmData",jzbvariabledata,100, -50,50, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets,data, luminosity);
26     TH1F *rawJZBemMC = allsamples.Draw("rawJZBemMC",jzbvariablemc,100,-50,50, "JZB (GeV)", "events", cutmass&&cutOSOF&&cutnJets,mc, luminosity);
27     TH1F *rawJZBemData = allsamples.Draw("rawJZBemData",jzbvariabledata,100, -50,50, "JZB (GeV)", "events", cutmass&&cutOSOF&&cutnJets,data, luminosity);
28 buchmann 1.1 TH1F *rawttbarjzbeemmMC;
29    
30     if(method==doKM) {
31     //we only need this histo for the KM fitting...
32 buchmann 1.8 rawttbarjzbeemmMC = allsamples.Draw("rawttbarjzbeemmMC",jzbvariablemc,100, -50,50, "JZB (GeV)", "events",cutmass&&cutOSSF&&cutnJets,mc,luminosity,allsamples.FindSample("TTJet"));
33 buchmann 1.1 MCPeak=find_peak(rawJZBeemmMC, rawttbarjzbeemmMC, -40, 40, mc, MCPeakError,MCSigma,method);
34     DataPeak=find_peak(rawJZBeemmData, rawJZBeemmData, -40, 40, data, DataPeakError,DataSigma,method);
35     }
36     else {
37     TH1F *reducedMC = (TH1F*)rawJZBeemmMC->Clone();
38     TH1F *reducedData = (TH1F*)rawJZBeemmData->Clone();
39     reducedMC->Add(rawJZBemMC,-1);
40     reducedData->Add(rawJZBemData,-1);
41     //this is Kostas' way of doing it - we subtract em to get rid of some of the ttbar contribution (in reality, of flavor-symmetric contribution)
42     MCPeak=find_peak(reducedMC, rawttbarjzbeemmMC, -40, 40, mc, MCPeakError,MCSigma,method);
43     DataPeak=find_peak(reducedData, rawJZBeemmData, -40, 40, data, DataPeakError,DataSigma,method);
44 buchmann 1.11
45 buchmann 1.1 }
46 buchmann 1.11
47 buchmann 1.1
48     // MCPeak=find_peak(rawJZBeemmMC, rawttbarjzbeemmMC, -40, 40, mc, MCPeakError,MCSigma,method);
49     // DataPeak=find_peak(rawJZBeemmData, rawJZBeemmData, -40, 40, data, DataPeakError,DataSigma,method);
50     cout << "We have found the peak in Data at " << DataPeak << " +/- " << DataPeakError << " with sigma=" << DataSigma << " +/- ?? (not impl.)" << endl;
51     result << "We have found the peak in Data at " << DataPeak << " +/- " << DataPeakError << " with sigma=" << DataSigma << " +/- ?? (not impl.)" << endl;
52     cout << "We have found the peak in MC at " << MCPeak << " +/- " << MCPeakError << " with sigma=" << MCSigma << " +/- ?? (not impl.)" << endl;
53     result << "We have found the peak in MC at " << MCPeak << " +/- " << MCPeakError << " with sigma=" << MCSigma << " +/- ?? (not impl.)" << endl;
54     }
55    
56     void make_special_mll_plot(int nbins, float min, float max, bool logscale,string xlabel) {
57    
58     TCanvas *ckin = new TCanvas("ckin","Kinematic Plots (in the making)",600,600);
59    
60     TH1F *datahistoOSSF = allsamples.Draw("datahistoOSSF","mll",nbins,min,max, xlabel, "events",cutOSSF&&cutnJets&&basiccut,data,luminosity);
61     THStack mcstackOSSF = allsamples.DrawStack("mcstackOSSF","mll",nbins,min,max, xlabel, "events",cutOSSF&&cutnJets&&basiccut,mc,luminosity);
62     TH1F *datahistoOSOF = allsamples.Draw("datahistoOSSF","mll",nbins,min,max, xlabel, "events",cutOSOF&&cutnJets&&basiccut,data,luminosity);
63    
64     if(logscale) ckin->SetLogy(1);
65     datahistoOSSF->GetXaxis()->SetTitle(xlabel.c_str());
66     datahistoOSSF->GetXaxis()->CenterTitle();
67     datahistoOSSF->GetYaxis()->SetTitle("events");
68     datahistoOSSF->GetYaxis()->CenterTitle();
69    
70     datahistoOSSF->Draw();
71     mcstackOSSF.Draw("same");
72     datahistoOSSF->Draw("same");
73     datahistoOSOF->SetMarkerColor(kOrange);
74     datahistoOSOF->Draw("same");
75     TLegend *kinleg = allsamples.allbglegend();
76     kinleg->AddEntry(datahistoOSOF,"OSOF (data)","p");
77     kinleg->Draw();
78     CompleteSave(ckin,"kin/mll_ossf_osof_distribution");
79     }
80    
81    
82     void make_plot(string variable, int nbins, float min, float max, bool logscale,string xlabel, string filename,bool isPF=true) {
83     // TCut basiccut("(pfJetGoodNum>=2&&pfJetGoodID[0])&&(pfJetGoodNum>=2&&pfJetGoodID[1])&&(passed_triggers||!is_data)");
84     TCut ibasiccut=basiccut;
85    
86     if(isPF) ibasiccut=basiccut&&"pfjzb[0]>-998";
87     //Step 1: Adapt the variable (if we're dealing with PF we need to adapt the variable!)
88     if(isPF) {
89     if(variable=="mll") variable="pfmll[0]";
90     if(variable=="jetpt[0]") variable="pfJetGoodPt[0]";
91     if(variable=="jeteta[0]") variable="pfJetGoodEta[0]";
92     if(variable=="pt") variable="pfpt[0]";
93     if(variable=="pt1") variable="pfpt1[0]";
94     if(variable=="eta1") variable="pfeta1[0]";
95     if(variable=="jzb[1]") variable="pfjzb[0]";
96     //if(variable=="pfJetGoodNum") variable="pfJetGoodNum"; // pointless.
97     }
98    
99     //Step 2: Refine the cut
100     TCanvas *ckin = new TCanvas("ckin","Kinematic Plots (in the making)",600,600);
101     ckin->SetLogy(logscale);
102     TCut cut;
103     cut=cutmass&&cutOSSF&&cutnJets&&ibasiccut;
104     if(filename=="nJets") cut=cutmass&&cutOSSF&&ibasiccut;
105     if(filename=="nJets_nocuts_except_mll_ossf") cut=cutmass&&cutOSSF;
106     if(filename=="mll") cut=cutOSSF&&cutnJets&&ibasiccut;
107     if(filename=="mll_inclusive") cut=cutOSSF;
108     if(filename=="pfJetGoodEta_0") cut=cutOSSF&&cutmass&&ibasiccut;
109    
110    
111    
112     TH1F *datahisto = allsamples.Draw("datahisto",variable,nbins,min,max, xlabel, "events",cut,data,luminosity);
113     THStack mcstack = allsamples.DrawStack("mcstack",variable,nbins,min,max, xlabel, "events",cut,mc,luminosity);
114     if(variable=="pfJetGoodPt[0]") datahisto->SetMaximum(10*datahisto->GetMaximum());
115     if(variable=="pt") datahisto->SetMaximum(10*datahisto->GetMaximum());
116     if(filename=="mll_inclusive") datahisto->SetMinimum(1);
117     // if(filename=="Zpt") datahisto->SetMaximum(300);
118     datahisto->Draw("e1");
119     mcstack.Draw("same");
120     datahisto->Draw("same,e1");
121     TLegend *kinleg = allsamples.allbglegend();
122     kinleg->Draw();
123     TText* write_cut = write_title(decipher_cut(cut,basicqualitycut));
124     // cout << "The file " << filename << " will have the cut " << (const char*)cut << " which is interpreted as " << decipher_cut(cut) << endl;
125     write_cut->Draw();
126     TText* write_variable = write_text(0.99,0.01,variable);
127     write_variable->SetTextAlign(31);
128     write_variable->SetTextSize(0.02);
129     write_variable->Draw();
130     if(isPF) CompleteSave(ckin,"kin/"+filename+"__PF");
131     else CompleteSave(ckin,"kin/"+filename);
132     datahisto->Delete();
133     }
134    
135     void do_kinematic_plots(bool doPF=false)
136     {
137     bool dolog=true;
138     bool nolog=false;
139     make_plot("mll",50,50,150,dolog,"m_{ll} (GeV)","mll",doPF);
140     make_plot("mll",100,50,150,dolog,"m_{ll} (GeV)","mll_inclusive",doPF);
141     make_plot("jetpt[0]",40,0,200,dolog,"leading jet p_{T} (GeV)","pfJetGoodPt_0",doPF);
142     make_plot("jeteta[0]",40,-5,5,nolog,"leading jet #eta","pfJetGoodEta_0",doPF);
143     make_plot("pt",50,0,200,dolog,"Z p_{T} (GeV)","Zpt",doPF);
144     make_plot("pt1",50,0,100,nolog,"p_{T} (GeV)","pt",doPF);
145     make_plot("eta1",40,-5,5,nolog,"#eta_{l}","eta",doPF);
146     make_plot("jzb[1]",100,-150,150,dolog,"JZB (GeV)","jzb_ossf",doPF);
147     make_plot("pfJetGoodNum",8,0.5,8.5,dolog,"nJets","nJets",doPF);
148     make_plot("pfJetGoodNum",8,0.5,8.5,dolog,"nJets","nJets_nocuts_except_mll_ossf",doPF);
149     make_special_mll_plot(30,50,150,dolog,"m_{ll} (GeV)");
150     }
151    
152     void do_kinematic_PF_plots()
153     {
154     do_kinematic_plots(true);
155     }
156    
157 buchmann 1.11 void signal_bg_comparison()
158 buchmann 1.1 {
159 buchmann 1.11 TCanvas *can = new TCanvas("can","Signal Background Comparison Canvas");
160 buchmann 1.1 int sbg_nbins=100;
161     float sbg_min=-200;
162     float sbg_max=400;
163     can->SetLogy(1);
164 buchmann 1.8 TH1F *JZBplotZJETs = allsamples.Draw("JZBplotZJETs",jzbvariablemc,sbg_nbins,sbg_min,sbg_max, "JZB (GeV)", "events",cutmass&&cutOSSF&&cutnJets,mc,luminosity,allsamples.FindSample("DYJetsToLL"));
165     TH1F *JZBplotLM4 = allsamples.Draw("JZBplotLM4",jzbvariablemc,sbg_nbins,sbg_min,sbg_max, "JZB (GeV)", "events",cutmass&&cutOSSF&&cutnJets,mc,luminosity,allsamples.FindSample("LM4"));
166     TH1F *JZBplotTtbar = allsamples.Draw("JZBplotTtbar",jzbvariablemc,sbg_nbins,sbg_min,sbg_max, "JZB (GeV)", "events",cutmass&&cutOSSF&&cutnJets,mc,luminosity,allsamples.FindSample("TTJets"));
167 buchmann 1.1
168     JZBplotTtbar->SetLineColor(allsamples.GetColor("TTJet"));
169     JZBplotTtbar->SetLineColor(allsamples.GetColor("TTJet"));
170     JZBplotTtbar->SetLineWidth(2);
171     JZBplotZJETs->SetFillColor(allsamples.GetColor("DY"));
172     JZBplotZJETs->SetLineColor(kBlack);
173     JZBplotLM4->SetLineWidth(2);
174     JZBplotLM4->SetLineStyle(2);
175     JZBplotZJETs->SetMinimum(0.05);
176     JZBplotLM4->Scale(20.0);
177     JZBplotZJETs->Draw("histo");
178     JZBplotLM4->Draw("histo,same");
179     JZBplotTtbar->Draw("histo,same");
180     TLegend *signal_bg_comparison_leg = make_legend("MC:");
181     signal_bg_comparison_leg->AddEntry(JZBplotZJETs,"Z+Jets","f");
182     signal_bg_comparison_leg->AddEntry(JZBplotTtbar,"t#bar{t}","f");
183     signal_bg_comparison_leg->AddEntry(JZBplotLM4,"20xLM4","f");
184     signal_bg_comparison_leg->Draw();
185     CompleteSave(can,"jzb_bg_vs_signal_distribution");
186     JZBplotLM4->Scale(1.0/20);
187     JZBplotZJETs->SetMinimum(0.004);
188     JZBplotZJETs->Draw("histo");
189     JZBplotLM4->Draw("histo,same");
190     JZBplotTtbar->Draw("histo,same");
191     TLegend *signal_bg_comparison_leg2 = make_legend("MC:");
192     signal_bg_comparison_leg2->AddEntry(JZBplotZJETs,"Z+Jets","f");
193     signal_bg_comparison_leg2->AddEntry(JZBplotTtbar,"t#bar{t}","f");
194     signal_bg_comparison_leg2->AddEntry(JZBplotLM4,"LM4","f");
195     signal_bg_comparison_leg2->Draw();
196     CompleteSave(can,"jzb_bg_vs_signal_distribution_unscaled_LM4");
197     }
198    
199 buchmann 1.9 //TF1 *BpredFunc,*BpredFuncP,*BpredFuncN;
200 buchmann 1.1
201 buchmann 1.9 vector<TF1*> do_cb_fit_to_plot(TH1F *histo, float Sigma, float doingfitacrosstheboard=false) {
202 buchmann 1.1
203     // BpredFunc = new TF1("BpredFunc",InvCrystalBall,0,histo->GetBinLowEdge(histo->GetNbinsX())+histo->GetBinWidth(histo->GetNbinsX()),5);
204 buchmann 1.9 TF1 *BpredFunc = new TF1("BpredFunc",InvCrystalBall,0,1000,5);
205 buchmann 1.1 BpredFunc->SetParameter(0,histo->GetBinContent(1));
206 buchmann 1.6 if(doingfitacrosstheboard) BpredFunc->SetParameter(0,histo->GetMaximum());
207 buchmann 1.1 BpredFunc->SetParameter(1,0.);
208     if(method==1) BpredFunc->SetParameter(2,10*Sigma);//KM
209     else BpredFunc->SetParameter(2,Sigma);//Gaussian based methods
210     if(method==-99) BpredFunc->SetParameter(2,2.0*Sigma);//Kostas
211     BpredFunc->SetParameter(3,1.8);
212     BpredFunc->SetParameter(4,2.5);
213 buchmann 1.11 histo->Fit(BpredFunc,"QN0");
214 buchmann 1.1 BpredFunc->SetLineColor(kBlue);
215     BpredFunc->SetLineWidth(1);
216    
217 buchmann 1.9 TF1 *BpredFuncP = new TF1("BpredFuncP",InvCrystalBallP,-1000,histo->GetBinLowEdge(histo->GetNbinsX())+histo->GetBinWidth(histo->GetNbinsX()),5);
218     TF1 *BpredFuncN = new TF1("BpredFuncN",InvCrystalBallN,-1000,histo->GetBinLowEdge(histo->GetNbinsX())+histo->GetBinWidth(histo->GetNbinsX()),5);
219 buchmann 1.1
220     BpredFuncP->SetParameters(BpredFunc->GetParameters());
221     BpredFuncP->SetLineColor(kBlue);
222     BpredFuncP->SetLineWidth(1);
223     BpredFuncP->SetLineStyle(2);
224    
225     BpredFuncN->SetParameters(BpredFunc->GetParameters());
226     BpredFuncN->SetLineColor(kBlue);
227     BpredFuncN->SetLineWidth(1);
228     BpredFuncN->SetLineStyle(2);
229 buchmann 1.9
230     vector<TF1*> functions;
231     functions.push_back(BpredFuncN);
232     functions.push_back(BpredFunc);
233     functions.push_back(BpredFuncP);
234     return functions;
235 buchmann 1.1 }
236     void do_prediction_plot(string jzb, TCanvas *globalcanvas, float sigma, float high, bool is_data, bool do_fit=true)
237     {
238     int nbins=100;
239     if(is_data) nbins=50;
240     float low=0;
241     float hi=500;
242    
243     TH1F *RcorrJZBeemm = allsamples.Draw("RcorrJZBeemm",jzb.c_str(),nbins,low,hi, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets,is_data, luminosity);
244     TH1F *LcorrJZBeemm = allsamples.Draw("LcorrJZBeemm",("-"+jzb).c_str(),nbins,low,hi, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets,is_data, luminosity);
245     TH1F *RcorrJZBem = allsamples.Draw("RcorrJZBem",jzb.c_str(),nbins,low,hi, "JZB (GeV)", "events", cutmass&&cutOSOF&&cutnJets,is_data, luminosity);
246     TH1F *LcorrJZBem = allsamples.Draw("LcorrJZBem",("-"+jzb).c_str(),nbins,low,hi, "JZB (GeV)", "events", cutmass&&cutOSOF&&cutnJets,is_data, luminosity);
247    
248     TH1F *Bpred = (TH1F*)LcorrJZBeemm->Clone();
249     Bpred->Add(RcorrJZBem);
250     Bpred->Add(LcorrJZBem,-1);
251     globalcanvas->cd();
252     globalcanvas->SetLogy(1);
253     Bpred->SetLineColor(kRed);
254     Bpred->SetStats(0);
255     Bpred->GetXaxis()->SetRangeUser(0,high);
256     Bpred->Draw("hist");
257     RcorrJZBeemm->SetMarkerStyle(20);
258     RcorrJZBeemm->Draw("e1x0,same");
259 buchmann 1.12 TLegend *legBpred = make_legend("",0.6,0.5);
260 buchmann 1.1
261     if(is_data)
262     {
263 buchmann 1.12 // legBpred->SetY1(0.5);
264     // legBpred->SetX1(0.6);
265 buchmann 1.9 vector<TF1*> functions = do_cb_fit_to_plot(Bpred,sigma);
266     functions[0]->Draw("same");functions[1]->Draw("same");functions[2]->Draw("same");
267 buchmann 1.12 legBpred->AddEntry(RcorrJZBeemm,"observed (data)","p");
268 buchmann 1.1 legBpred->AddEntry(Bpred,"predicted (data)","l");
269 buchmann 1.9 legBpred->AddEntry(functions[0],"predicted fit (data)","l");
270     legBpred->AddEntry(functions[1],"stat. uncert.","l");
271 buchmann 1.1 legBpred->Draw();
272     CompleteSave(globalcanvas,"Bpred_Data");
273     }
274     else {
275     legBpred->AddEntry(RcorrJZBeemm,"MC true B","l");
276     legBpred->AddEntry(Bpred,"data-driven B","l");
277     legBpred->Draw();
278     CompleteSave(globalcanvas,"Bpred_MC");
279     }
280     }
281    
282 buchmann 1.11 void do_prediction_plots(string mcjzb, string datajzb, float DataSigma, float MCSigma) {
283     TCanvas *globalcanvas = new TCanvas("globalcanvas","Prediction Canvas");
284 buchmann 1.1 do_prediction_plot(datajzb,globalcanvas,DataSigma,350,data);
285     do_prediction_plot(mcjzb,globalcanvas,MCSigma,300,mc);
286     }
287    
288     void do_ratio_plot(int is_data,vector<float> binning, string jzb, TCanvas *can, float high=-9999) {
289     bool do_data=0;
290     bool dosignal=0;
291     if(is_data==1) do_data=1;
292     if(is_data==2) dosignal=1;
293     TH1F *RcorrJZBeemm = allsamples.Draw("RcorrJZBeemm",jzb.c_str(),binning, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets,do_data, luminosity,dosignal);
294     TH1F *LcorrJZBeemm = allsamples.Draw("LcorrJZBeemm",("-"+jzb).c_str(),binning, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets,do_data, luminosity,dosignal);
295     TH1F *RcorrJZBem = allsamples.Draw("RcorrJZBem",jzb.c_str(),binning, "JZB (GeV)", "events", cutmass&&cutOSOF&&cutnJets,do_data, luminosity,dosignal);
296     TH1F *LcorrJZBem = allsamples.Draw("LcorrJZBem",("-"+jzb).c_str(),binning, "JZB (GeV)", "events", cutmass&&cutOSOF&&cutnJets,do_data, luminosity,dosignal);
297    
298     TH1F *Bpred = (TH1F*)LcorrJZBeemm->Clone();
299     Bpred->Add(RcorrJZBem);
300     Bpred->Add(LcorrJZBem,-1);
301     can->cd();
302     can->SetLogy(0);
303     Bpred->SetLineColor(kRed);
304     Bpred->SetStats(0);
305     if(high>0) Bpred->GetXaxis()->SetRangeUser(0,high);
306     TGraphAsymmErrors *JZBratio = histRatio(RcorrJZBeemm,Bpred,is_data,binning);
307     TH1F *JZBratioforfitting=(TH1F*)RcorrJZBeemm->Clone();
308     JZBratioforfitting->Divide(Bpred);
309    
310     JZBratio->SetTitle("");
311     JZBratio->GetYaxis()->SetRangeUser(0.0,9.0);
312     if(is_data==1) JZBratio->GetXaxis()->SetRangeUser(0,350);
313    
314     TF1 *pol0 = new TF1("pol0","[0]",0,1000);
315     TF1 *pol0d = new TF1("pol0","[0]",0,1000);
316     // straightline_fit->SetParameter(0,1);
317     JZBratioforfitting->Fit(pol0,"Q0R","",0,30);
318     pol0d->SetParameter(0,pol0->GetParameter(0));
319    
320     JZBratio->GetYaxis()->SetTitle("Observed/Predicted");
321     JZBratio->GetXaxis()->SetTitle("JZB (GeV)");
322     JZBratio->GetYaxis()->SetNdivisions(519);
323     JZBratio->GetYaxis()->SetRangeUser(0.0,9.0);
324     JZBratio->Draw("AP");
325     /////----------------------------
326     TPaveText *writeresult = new TPaveText(0.2,0.78,0.49,0.91,"blNDC");
327     writeresult->SetFillStyle(4000);
328     writeresult->SetFillColor(kWhite);
329     writeresult->SetTextFont(42);
330     ostringstream jzb_agreement_data_text;
331     jzb_agreement_data_text<< setprecision(2) << "mean =" << pol0->GetParameter(0) << " #pm " << setprecision(1) << pol0->GetParError(0);
332     if(is_data==1) fitresultconstdata=pol0->GetParameter(0);// data
333     if(is_data==0) fitresultconstmc=pol0->GetParameter(0); // monte carlo, no signal
334     /* if(is_data) writeresult->AddText("Data closure test");
335     else writeresult->AddText("MC closure test");
336     */
337     writeresult->AddText(jzb_agreement_data_text.str().c_str());
338     writeresult->Draw("same");
339     pol0d->Draw("same");
340     TF1 *topline = new TF1("","1.2",0,1000);
341     TF1 *bottomline = new TF1("","0.6",0,1000);
342     topline->SetLineColor(kBlue);
343     topline->SetLineStyle(2);
344     bottomline->SetLineColor(kBlue);
345     bottomline->SetLineStyle(2);
346     topline->Draw("same");
347     bottomline->Draw("same");
348     TF1 *oneline = new TF1("","1.0",0,1000);
349     oneline->SetLineColor(kBlue);
350     oneline->SetLineStyle(1);
351     oneline->Draw("same");
352     TLegend *phony_leg = make_legend("ratio");//this line is just to have the default CMS Preliminary bla on the canvas as well.
353     TLegend *leg = new TLegend(0.60,0.5,0.94,0.77);
354     leg->SetTextFont(42);
355     if(is_data==1) leg->SetHeader("Ratio (data)");// = make_legend("Ratio (data)");
356     else leg->SetHeader("Ratio (MC)"); //leg = make_legend("Ratio (MC)");
357    
358     leg->SetFillStyle(4000);
359     leg->SetFillColor(kWhite);
360     leg->SetTextFont(42);
361     // leg->AddEntry(topline,"+20\% sys envelope","l");
362     leg->AddEntry(JZBratio,"obs/pred","p");
363     leg->AddEntry(oneline,"ratio = 1","l");
364     leg->AddEntry(pol0d,"fit in [0,30] GeV","l");
365     leg->AddEntry(bottomline,"sys. envelope","l");
366     leg->Draw("same");
367     if(is_data==1) CompleteSave(can, "jzb_ratio_data");
368     if(is_data==0) CompleteSave(can, "jzb_ratio_mc");
369     if(is_data==2) CompleteSave(can, "jzb_ratio_mc_BandS");//special case, MC with signal!
370    
371     delete RcorrJZBeemm;
372     delete LcorrJZBeemm;
373     delete RcorrJZBem;
374     delete LcorrJZBem;
375     }
376    
377 buchmann 1.11 void do_ratio_plots(string mcjzb,string datajzb) {
378     TCanvas *globalc = new TCanvas("globalc","Ratio Plot Canvas");
379 buchmann 1.1 globalc->SetLogy(0);
380     vector<float> ratio_binning;
381     ratio_binning.push_back(0);
382     ratio_binning.push_back(5);
383     ratio_binning.push_back(10);
384     ratio_binning.push_back(20);
385     ratio_binning.push_back(50);
386     ratio_binning.push_back(100);
387     ratio_binning.push_back(350);
388     ratio_binning.push_back(500);
389    
390     do_ratio_plot(mc,ratio_binning,mcjzb,globalc);
391     do_ratio_plot(data,ratio_binning,datajzb,globalc);
392     do_ratio_plot(mcwithsignal,ratio_binning,mcjzb,globalc);
393     }
394    
395 buchmann 1.8 string give_jzb_expression(float peak, int type) {
396 buchmann 1.1 stringstream val;
397 buchmann 1.8 if(type==data) {
398     if(peak<0) val << jzbvariabledata << "+" << TMath::Abs(peak);
399     if(peak>0) val << jzbvariabledata << "-" << TMath::Abs(peak);
400     if(peak==0) val << jzbvariabledata;
401     }
402     if(type==mc) {
403     if(peak<0) val << jzbvariablemc << "+" << TMath::Abs(peak);
404     if(peak>0) val << jzbvariablemc << "-" << TMath::Abs(peak);
405     if(peak==0) val << jzbvariablemc;
406     }
407 buchmann 1.1 return val.str();
408     }
409    
410     string centralupdownname(int num) {
411     if (num==0) return "central";
412     if (num==1) return "down";
413     if (num==2) return "up";
414     return "centralupdownnameERROR";
415     }
416    
417     string dataormc(int isdata) {
418     if(isdata) return "Data";
419     else return "MC";
420     }
421    
422     int histocounter=0;
423     string give_histo_number(string basename="", int isdata=-1, int centralcounter=-1) {
424     histocounter++;
425     stringstream histocounternum;
426     if(isdata==-1) {
427     if(basename=="") histocounternum << histocounter;
428     else histocounternum << basename << "_" << histocounter;
429     }
430     else {
431     histocounternum << basename << "_" << dataormc(isdata) << "_" << centralupdownname(centralcounter);
432     }
433     return histocounternum.str();
434     }
435    
436     int central=0;
437     int down=1;
438     int up=2;
439    
440     void crunch_the_numbers(TH1F *RcorrJZBeemmop[3][2],TH1F *RcorrJZBeeop[3][2],TH1F *RcorrJZBmmop[3][2],TH1F *LcorrJZBeemmop[3][2],TH1F *LcorrJZBeeop[3][2],TH1F *LcorrJZBmmop[3][2],TH1F *RcorrJZBemop[3][2],TH1F *LcorrJZBemop[3][2],float zjetsestimate[3][2][20],float ttbarestimate[3][2][2][20],float predicted[3][2][20],float observed[3][2][20],int isdata,vector<float> jzbcuts,float eemmmumu[2][3][20]) {
441     // cout << "Crunching the numbers for is_data=" << isdata << endl;
442     for(int k=0;k<20;k++) {for(int i=0;i<2;i++){for(int j=0;j<3;j++) {eemmmumu[i][j][k]=0;}}}
443     for(int icut=0;icut<jzbcuts.size()-1;icut++) {//do calculation for each JZB cut
444     for(int ibin=1;ibin<jzbcuts.size();ibin++) {//and to do that, we need to consider each bin.
445     if(RcorrJZBeemmop[0][isdata]->GetBinLowEdge(ibin)<jzbcuts[icut]) continue;
446     // if(icut==0) cout << "Predicted+=" << LcorrJZBeemmop[central][isdata]->GetBinContent(ibin) << " eemm,L" << endl;
447     predicted[central][isdata][icut]+=LcorrJZBeemmop[central][isdata]->GetBinContent(ibin);
448     predicted[up][isdata][icut]+=LcorrJZBeemmop[up][isdata]->GetBinContent(ibin);
449     predicted[down][isdata][icut]+=LcorrJZBeemmop[down][isdata]->GetBinContent(ibin);
450    
451     // if(icut==0) cout << "Predicted+=" << RcorrJZBemop[central][isdata]->GetBinContent(ibin) << " em, R" << endl;
452     predicted[central][isdata][icut]+=RcorrJZBemop[central][isdata]->GetBinContent(ibin);
453     predicted[up][isdata][icut]+=RcorrJZBemop[up][isdata]->GetBinContent(ibin);
454     predicted[down][isdata][icut]+=RcorrJZBemop[down][isdata]->GetBinContent(ibin);
455    
456     // if(icut==0) cout << "Predicted-=" << LcorrJZBemop[central][isdata]->GetBinContent(ibin) <<" em, L" << endl;
457     predicted[central][isdata][icut]-=LcorrJZBemop[central][isdata]->GetBinContent(ibin);
458     predicted[up][isdata][icut]-=LcorrJZBemop[up][isdata]->GetBinContent(ibin);
459     predicted[down][isdata][icut]-=LcorrJZBemop[down][isdata]->GetBinContent(ibin);
460    
461     // if(icut==0) cout << "Observed=" << RcorrJZBeemmop[central][isdata]->GetBinContent(ibin) << "eemm, R" << endl;
462     observed[central][isdata][icut]+=RcorrJZBeemmop[central][isdata]->GetBinContent(ibin);
463     observed[up][isdata][icut]+=RcorrJZBeemmop[up][isdata]->GetBinContent(ibin);
464     observed[down][isdata][icut]+=RcorrJZBeemmop[down][isdata]->GetBinContent(ibin);
465    
466     zjetsestimate[central][isdata][icut]+=LcorrJZBeemmop[central][isdata]->GetBinContent(ibin);
467     ttbarestimate[central][isdata][0][icut]+=RcorrJZBemop[central][isdata]->GetBinContent(ibin);
468     ttbarestimate[central][isdata][1][icut]+=LcorrJZBemop[central][isdata]->GetBinContent(ibin);
469    
470     eemmmumu[0][0][icut]+=LcorrJZBeeop[central][isdata]->GetBinContent(ibin);
471     eemmmumu[0][1][icut]+=LcorrJZBmmop[central][isdata]->GetBinContent(ibin);
472     eemmmumu[0][2][icut]+=LcorrJZBemop[central][isdata]->GetBinContent(ibin);
473    
474     eemmmumu[1][0][icut]+=RcorrJZBeeop[central][isdata]->GetBinContent(ibin);
475     eemmmumu[1][1][icut]+=RcorrJZBmmop[central][isdata]->GetBinContent(ibin);
476     eemmmumu[1][2][icut]+=RcorrJZBemop[central][isdata]->GetBinContent(ibin);
477    
478     }//end of ibin loop
479     }//end of icut loop
480     }
481    
482     void make_table(float eemmmumu[2][3][20],int icut,float JZBcut) {
483     vector< vector <string> > entries;
484     vector <string> line;
485     line.push_back("");
486     line.push_back("eemm (ee/mm)");
487     line.push_back("em");
488     entries.push_back(line);
489     line.clear();
490     line.push_back("JZB>"+any2string(JZBcut));
491     line.push_back(any2string(eemmmumu[1][0][icut]+eemmmumu[1][1][icut])+"("+any2string(eemmmumu[1][0][icut])+"/"+any2string(eemmmumu[1][1][icut])+")");
492     line.push_back(any2string(eemmmumu[1][2][icut]));
493     entries.push_back(line);
494     line.clear();
495     line.push_back("JZB<-"+any2string(JZBcut));
496     line.push_back(any2string(eemmmumu[0][0][icut]+eemmmumu[0][1][icut])+"("+any2string(eemmmumu[0][0][icut])+"/"+any2string(eemmmumu[0][1][icut])+")");
497     line.push_back(any2string(eemmmumu[0][2][icut]));
498     entries.push_back(line);
499     make_nice_jzb_table(entries);
500     /*
501     cout << " \t\t | \t eemm (ee/mm) \t | \t em" << endl;
502     cout << "JZB>peak\t | \t " << eemmmumu[1][0]+eemmmumu[1][1] << "(" << eemmmumu[1][0] << "/" << eemmmumu[1][1] << ")\t | \t" << eemmmumu[1][2] << endl;
503     cout << "JZB<peak\t | \t " << eemmmumu[0][0]+eemmmumu[0][1] << "(" << eemmmumu[0][0] << "/" << eemmmumu[0][1] << ")\t | \t" << eemmmumu[0][2] << endl;
504     cout << endl;
505     */
506     }
507    
508     void present_result(vector<float> &jzbcuts,float predicted[3][2][20],float observed[3][2][20],float zjetsestimate[3][2][20],float ttbarestimate[3][2][2][20],int icut, float eemmmumu[2][3][20]) {
509    
510     //blublu
511     cout << endl << endl;
512     cout << "--------------------------------------------------------------------------------------" << endl;
513     cout << "DATA: " << endl;
514     cout << " For JZB>" << jzbcuts[icut];
515     float max,downvar,upvar;
516     stringstream printtitle;
517     printtitle << "JZB>" << jzbcuts[icut] << " (data)";
518     ComputePoissonError(zjetsestimate[central][data][icut],ttbarestimate[central][data][0][icut],ttbarestimate[central][data][1][icut],max,downvar,upvar,printtitle.str());
519     //cout << " Predicted: " << print_range(predicted[central][data][icut],predicted[up][data][icut],predicted[down][data][icut]) << " (stat) +" << 0.2*zjetsestimate[central][data][icut] << "-" << 0.4*zjetsestimate[central][data][icut] << " (sys) " << " using Poisson: " << print_range(max,max+upvar,max-downvar) << ")" << endl;
520     float pred=predicted[central][data][icut];
521     float sysP=abs(predicted[up][data][icut]-pred);
522     float sysN=abs(pred-predicted[down][data][icut]);
523     sysN = sysN + pred*(1-fitresultconstdata);//fitresultconst comes from the fit in the 0-30 GeV range!
524     cout << " Predicted: " << pred << "+" << statErrorP(pred) << "-" << statErrorN(pred) << " (stat) +" << sysP << " - " << sysN << " (sys) " << " using Poisson: " << print_range(max,max+upvar,max-downvar) << ")" << endl;
525     //0.2*zjetsestimate[central][data][icut] << "-" << 0.4*zjetsestimate[central][data][icut] << " (sys) " << " using Poisson: " << print_range(max,max+upvar,max-downvar) << ")" << endl;
526     cout << " Details: ZJetsEstimate= " << zjetsestimate[central][data][icut] << ", TTbar estimate=" << ttbarestimate[central][data][0][icut]-ttbarestimate[central][data][1][icut] << " [" << ttbarestimate[central][data][0][icut] << " R -" << ttbarestimate[central][data][1][icut] << " L ]" << endl;
527     cout << " For JZB>" << jzbcuts[icut] << " Observed: " << observed[central][data][icut] << endl;
528     cout << "TABLE:" << endl;
529     make_table(eemmmumu,icut,jzbcuts[icut]);
530     cout << "MC: " << endl;
531     printtitle.str("");
532     printtitle<<"JZB>"<<jzbcuts[icut]<<" (MC)";
533     ComputePoissonError(zjetsestimate[central][mc][icut],ttbarestimate[central][mc][0][icut],ttbarestimate[central][mc][1][icut],max,downvar,upvar,printtitle.str());
534     pred=predicted[central][mc][icut];
535     sysP=abs(predicted[up][mc][icut]-pred);
536     sysN=abs(pred-predicted[down][mc][icut]);
537     sysN = sysN + pred*(1-fitresultconstmc);//fitresultconst comes from the fit in the 0-30 GeV range!
538     cout << " Predicted: " << pred << "+" << statErrorP(pred) << "-" << statErrorN(pred) << " (stat) +" << sysP << " - " << sysN << " (sys) " << " using Poisson: " << print_range(max,max+upvar,max-downvar) << ")" << endl;
539     // cout << " Predicted: " << print_range(pred,statErrorP(pred),statErrorN(pred)) << " (stat) +" << sysP << " - " << sysN << " (sys) " << " using Poisson: " << print_range(max,max+upvar,max-downvar) << ")" << endl;
540     //cout << " Predicted: " << print_range(predicted[central][mc][icut],predicted[down][mc][icut],predicted[up][mc][icut]) << " (stat) +" << 0.2*zjetsestimate[central][mc][icut] << "-" << 0.4*zjetsestimate[central][mc][icut] << " (sys) " << " using Poisson: " << print_range(max,max+upvar,max-downvar) << ")" << endl;
541     //print_range(predicted[central][mc][icut],predicted[central][mc][icut]+upvar,predicted[central][mc][icut]-downvar) << ")" << endl;
542     cout << " Details: ZJetsEstimate= " << zjetsestimate[central][mc][icut] << ", TTbar estimate=" << ttbarestimate[central][mc][0][icut]-ttbarestimate[central][mc][1][icut] << " [" << ttbarestimate[central][mc][0][icut] << " R -" << ttbarestimate[central][mc][1][icut] << " L ] )" << endl;
543     cout << " For JZB>" << jzbcuts[icut] << " Observed: " << observed[central][mc][icut] << endl;
544     cout << endl << endl;
545     }
546    
547     void calculate_predicted_and_observed_eemm(float MCPeak,float MCPeakError,float DataPeak,float DataPeakError,vector<float> jzbcuts) { // DO NOT CHANGE THIS TO ... &jzbcuts !!! we add one!
548     jzbcuts.push_back(14000);
549    
550     /*
551     * We want the following numbers: For all JZB cuts, we want: ee,mm,em; observed, predicted. we want final results with errors from the peak and statistical error.
552     * How to accomplish this; we draw histos for MC&data, once with the peak (to extract obs/pred), once with peak+sigma, once with peak-sigma. this gives us the peak error.
553     * Statistical error: Comes from Poissonian errors ... for this we have a special little macro that will give us the value for three given parameters (cool, right?)
554     */
555    
556     string jzbpeak[3][2];
557 buchmann 1.8 jzbpeak[central][data]=give_jzb_expression(DataPeak,data);
558     jzbpeak[up][data]=give_jzb_expression(DataPeak+DataPeakError,data);
559     jzbpeak[down][data]=give_jzb_expression(DataPeak-DataPeakError,data);
560     jzbpeak[central][mc]=give_jzb_expression(MCPeak,mc);
561     jzbpeak[up][mc]=give_jzb_expression(MCPeak+MCPeakError,mc);
562     jzbpeak[down][mc]=give_jzb_expression(MCPeak-MCPeakError,mc);
563 buchmann 1.1
564     TH1F *RcorrJZBeemmop[3][2];
565     TH1F *LcorrJZBeemmop[3][2];
566     TH1F *RcorrJZBeeop[3][2];
567     TH1F *LcorrJZBeeop[3][2];
568     TH1F *RcorrJZBmmop[3][2];
569     TH1F *LcorrJZBmmop[3][2];
570     TH1F *RcorrJZBemop[3][2];
571     TH1F *LcorrJZBemop[3][2];
572     float zjetsestimate[3][2][20];
573     float predicted[3][2][20];
574     float observed[3][2][20];
575     float ttbarestimate[3][2][2][20];
576     float eemmmumu[2][3][20];
577    
578    
579     for(int isdata=0;isdata<=1;isdata++) {
580     for(int centraldownup=0;centraldownup<=2;centraldownup++) {
581     RcorrJZBeemmop[centraldownup][isdata] = allsamples.Draw(give_histo_number("RcorrJZBeemmop",isdata,centraldownup),jzbpeak[centraldownup][isdata].c_str(),jzbcuts, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets,isdata, luminosity);
582     LcorrJZBeemmop[centraldownup][isdata] = allsamples.Draw(give_histo_number("LcorrJZBeemmop",isdata,centraldownup),("-"+jzbpeak[centraldownup][isdata]).c_str(),jzbcuts, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets,isdata, luminosity);
583     RcorrJZBemop[centraldownup][isdata] = allsamples.Draw(give_histo_number("RcorrJZBemop",isdata,centraldownup),jzbpeak[centraldownup][isdata].c_str(),jzbcuts, "JZB (GeV)", "events", cutmass&&cutOSOF&&cutnJets,isdata, luminosity);
584     LcorrJZBemop[centraldownup][isdata] = allsamples.Draw(give_histo_number("LcorrJZBemop",isdata,centraldownup),("-"+jzbpeak[centraldownup][isdata]).c_str(),jzbcuts, "JZB (GeV)", "events", cutmass&&cutOSOF&&cutnJets,isdata,luminosity);
585    
586     RcorrJZBeeop[centraldownup][isdata] = allsamples.Draw(give_histo_number("RcorrJZBeeop",isdata,centraldownup),jzbpeak[centraldownup][isdata].c_str(),jzbcuts, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets&&"id1==0",isdata, luminosity);
587     LcorrJZBeeop[centraldownup][isdata] = allsamples.Draw(give_histo_number("LcorrJZBeeop",isdata,centraldownup),("-"+jzbpeak[centraldownup][isdata]).c_str(),jzbcuts, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets&&"id1==0",isdata, luminosity);
588     RcorrJZBmmop[centraldownup][isdata] = allsamples.Draw(give_histo_number("RcorrJZBmmop",isdata,centraldownup),jzbpeak[centraldownup][isdata].c_str(),jzbcuts, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets&&"id1==1",isdata, luminosity);
589     LcorrJZBmmop[centraldownup][isdata] = allsamples.Draw(give_histo_number("LcorrJZBmmop",isdata,centraldownup),("-"+jzbpeak[centraldownup][isdata]).c_str(),jzbcuts, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets&&"id1==1",isdata, luminosity);
590     for(int icut=0;icut<=jzbcuts.size();icut++) {
591     zjetsestimate[centraldownup][isdata][icut]=0;
592     predicted[centraldownup][isdata][icut]=0;
593     observed[centraldownup][isdata][icut]=0;
594     ttbarestimate[centraldownup][isdata][0][icut]=0;
595     ttbarestimate[centraldownup][isdata][1][icut]=0;
596     }
597     }//end of for central,up,down loop
598     crunch_the_numbers(RcorrJZBeemmop,RcorrJZBeeop,RcorrJZBmmop,LcorrJZBeemmop,LcorrJZBeeop,LcorrJZBmmop,RcorrJZBemop,LcorrJZBemop,zjetsestimate,ttbarestimate,predicted,observed,isdata,jzbcuts,eemmmumu);
599    
600     // void crunch_the_numbers(TH1F *RcorrJZBeemmop[3][2],TH1F *RcorrJZBeeop[3][2],TH1F *RcorrJZBmmop[3][2],TH1F *LcorrJZBeemmop[3][2],TH1F *LcorrJZBeeop[3][2],TH1F *LcorrJZBmmop[3][2],TH1F *RcorrJZBemop[3][2],TH1F *LcorrJZBemop[3][2],float zjetsestimate[3][2][20],float ttbarestimate[3][2][2][20],float predicted[3][2][20],float observed[3][2][20],int isdata,vector<float> jzbcuts,float eemmmumu[2][3]) {
601    
602     }//end of is data loop
603    
604    
605     cout << "We obtain the following results: " << endl;
606     for(int icut=0;icut<jzbcuts.size()-1;icut++) {
607     present_result(jzbcuts,predicted,observed,zjetsestimate,ttbarestimate,icut,eemmmumu);
608     }
609    
610    
611     }
612    
613 buchmann 1.11 void lepton_comparison_plots() {
614     TCanvas *can = new TCanvas("can","Lepton Comparison Canvas");
615 buchmann 1.1 can->SetLogy(1);
616     TH1F *eemc = allsamples.Draw("eemc","mll",50,50,150, "mll (GeV)", "events", cutOSSF&&cutnJets&&"(id1==0)",mc, luminosity,allsamples.FindSample("DYJetsToLL"));
617     TH1F *mmmc = allsamples.Draw("mmmc","mll",50,50,150, "mll (GeV)", "events", cutOSSF&&cutnJets&&"(id1==1)",mc, luminosity,allsamples.FindSample("DYJetsToLL"));
618     eemc->SetLineColor(kBlue);
619     mmmc->SetLineColor(kRed);
620     eemc->SetMinimum(0.1);
621     eemc->SetMaximum(10*eemc->GetMaximum());
622     eemc->Draw("histo");
623     mmmc->Draw("histo,same");
624     TLegend *leg = make_legend();
625     leg->AddEntry(eemc,"ZJets->ee (MC)","f");
626     leg->AddEntry(mmmc,"ZJets->#mu#mu (MC)","f");
627 buchmann 1.11 leg->Draw("same");
628 buchmann 1.1 CompleteSave(can, "mll_effratio_mc");
629    
630     TH1F *eed = allsamples.Draw("eed","mll",50,50,150, "mll (GeV)", "events", cutOSSF&&cutnJets&&"(id1==0)",data, luminosity);
631     TH1F *mmd = allsamples.Draw("mmd","mll",50,50,150, "mll (GeV)", "events", cutOSSF&&cutnJets&&"(id1==1)",data, luminosity);
632     eed->SetLineColor(kBlue);
633     mmd->SetLineColor(kRed);
634     eed->SetMinimum(0.1);
635     eed->SetMaximum(10*eed->GetMaximum());
636     eed->Draw("histo");
637     mmd->Draw("histo,same");
638     TLegend *leg2 = make_legend();
639     leg2->AddEntry(eed,"ZJets->ee (data)","f");
640     leg2->AddEntry(mmd,"ZJets->#mu#mu (data)","f");
641     leg2->Draw();
642     CompleteSave(can, "mll_effratio_data");
643    
644 buchmann 1.8 TH1F *jeed = allsamples.Draw("jeed",jzbvariabledata, 90,-100,350, "JZB (GeV)", "events", cutOSSF&&cutnJets&&"(id1==0)",data, luminosity);
645     TH1F *jmmd = allsamples.Draw("jmmd",jzbvariabledata, 90,-100,350, "JZB (GeV)", "events", cutOSSF&&cutnJets&&"(id1==1)",data, luminosity);
646     TH1F *jeemmd = allsamples.Draw("jeemmd",jzbvariabledata,90,-100,350, "JZB (GeV)", "events", cutOSSF&&cutnJets,data, luminosity);
647 buchmann 1.11 cout << "ee : " << jeed->GetMean() << "+/-" << jeed->GetMeanError() << endl;
648     cout << "ee : " << jmmd->GetMean() << "+/-" << jmmd->GetMeanError() << endl;
649     cout << "eemd : " << jeemmd->GetMean() << "+/-" << jeemmd->GetMeanError() << endl;
650 buchmann 1.1 jeemmd->SetLineColor(kBlack);
651     jeemmd->SetMarkerStyle(25);
652     jeed->SetLineColor(kBlue);
653     jmmd->SetLineColor(kRed);
654     jeed->SetMinimum(0.1);
655     jeed->SetMaximum(10*eed->GetMaximum());
656     jeemmd->DrawNormalized();
657     jeed->DrawNormalized("histo,same");
658     jmmd->DrawNormalized("histo,same");
659     jeemmd->DrawNormalized("same");
660     TLegend *jleg2 = make_legend("Data");
661     jleg2->AddEntry(jeemmd,"ee and #mu#mu","p");
662     jleg2->AddEntry(jeed,"ee","l");
663     jleg2->AddEntry(jmmd,"#mu#mu","l");
664     jleg2->Draw();
665     CompleteSave(can,"jzb_effratio_data");
666    
667 buchmann 1.11 TH1F *zjeed = allsamples.Draw("zjeed",jzbvariabledata, 90,-100,350, "JZB (GeV)", "events", cutOSSF&&cutnJets&&"(id1==0)",mc, luminosity,allsamples.FindSample("DYJetsToLL_TuneZ2"));
668     TH1F *zjmmd = allsamples.Draw("zjmmd",jzbvariabledata, 90,-100,350, "JZB (GeV)", "events", cutOSSF&&cutnJets&&"(id1==1)",mc, luminosity,allsamples.FindSample("DYJetsToLL_TuneZ2"));
669     TH1F *zjeemmd = allsamples.Draw("zjeemmd",jzbvariabledata,90,-100,350, "JZB (GeV)", "events", cutOSSF&&cutnJets,mc, luminosity,allsamples.FindSample("DYJetsToLL_TuneZ2"));
670     cout << "Z+Jets ee : " << zjeed->GetMean() << "+/-" << zjeed->GetMeanError() << endl;
671     cout << "Z+Jets ee : " << zjmmd->GetMean() << "+/-" << zjmmd->GetMeanError() << endl;
672     cout << "Z+Jets eemd : " << zjeemmd->GetMean() << "+/-" << zjeemmd->GetMeanError() << endl;
673     zjeemmd->SetLineColor(kBlack);
674     zjeemmd->SetMarkerStyle(25);
675     zjeed->SetLineColor(kBlue);
676     zjmmd->SetLineColor(kRed);
677     zjeed->SetMinimum(0.1);
678     zjeed->SetMaximum(10*eed->GetMaximum());
679     zjeemmd->DrawNormalized();
680     zjeed->DrawNormalized("histo,same");
681     zjmmd->DrawNormalized("histo,same");
682     zjeemmd->DrawNormalized("same");
683     TLegend *zjleg2 = make_legend("Data");
684     zjleg2->AddEntry(jeemmd,"ee and #mu#mu","p");
685     zjleg2->AddEntry(jeed,"ee","l");
686     zjleg2->AddEntry(jmmd,"#mu#mu","l");
687     zjleg2->Draw();
688     CompleteSave(can,"jzb_effratio_ZJets");
689    
690 buchmann 1.1 TH1F *ld = allsamples.Draw("ld","mll",50,50,150, "mll (GeV)", "events", cutOSSF&&cutnJets,data, luminosity);
691     ld->DrawNormalized("e1");
692     eed->DrawNormalized("histo,same");
693     mmd->DrawNormalized("histo,same");
694     TLegend *leg3 = make_legend();
695     leg3->AddEntry(ld,"ZJets->ll (data)","p");
696     leg3->AddEntry(eed,"ZJets->ee (data)","f");
697     leg3->AddEntry(mmd,"ZJets->#mu#mu (data)","f");
698     leg3->Draw();
699     CompleteSave(can,"mll_effratio_data__all_compared");
700     /*
701     TH1F *jzbld = allsamples.Draw("jzbld",jzbvariable,75,-150,150, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets,data, luminosity);
702     TH1F *jzbemd = allsamples.Draw("jzbemd",jzbvariable,75,-150,150, "JZB (GeV)", "events", cutmass&&cutOSOF&&cutnJets,data, luminosity);
703     */
704 buchmann 1.8 TH1F *jzbld = allsamples.Draw("jzbld",jzbvariabledata,92,-110,350, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets,data, luminosity);
705     TH1F *jzbemd = allsamples.Draw("jzbemd",jzbvariabledata,92,-110,350, "JZB (GeV)", "events", cutmass&&cutOSOF&&cutnJets,data, luminosity);
706 buchmann 1.1 jzbld->SetMarkerColor(kBlack);
707     jzbld->SetMarkerStyle(26);
708     jzbemd->SetMarkerStyle(25);
709     jzbemd->SetMarkerColor(kRed);
710     jzbemd->SetLineColor(kRed);
711     jzbld->SetMinimum(0.35);
712     jzbld->Draw("e1");
713     jzbemd->Draw("e1,same");
714     TLegend *leg4 = make_legend("Data");
715     leg4->AddEntry(jzbld,"ee or #mu#mu","p");
716     leg4->AddEntry(jzbemd,"e#mu","p");
717     leg4->Draw();
718     CompleteSave(can,"jzb_eemumu_emu_data");
719    
720 buchmann 1.8 TH1F *ttbarjzbld = allsamples.Draw("ttbarjzbld",jzbvariablemc,110,-150,400, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets,data, luminosity,allsamples.FindSample("TTJet"));
721     TH1F *ttbarjzbemd = allsamples.Draw("ttbarjzbemd",jzbvariablemc,110,-150,400, "JZB (GeV)", "events", cutmass&&cutOSOF&&cutnJets,data, luminosity,allsamples.FindSample("TTJet"));
722 buchmann 1.1 ttbarjzbld->SetLineColor(allsamples.GetColor("TTJet"));
723     ttbarjzbemd->SetLineColor(allsamples.GetColor("TTJet"));
724     ttbarjzbld->SetLineWidth(2);
725     ttbarjzbemd->SetLineWidth(2);
726     ttbarjzbld->Draw("histo");
727     ttbarjzbemd->SetLineStyle(2);
728     ttbarjzbemd->Draw("histo,same");
729     TLegend *leg5 = make_legend();
730     leg5->AddEntry(ttbarjzbld,"t#bar{t}->(ee or #mu#mu)","f");
731     leg5->AddEntry(ttbarjzbemd,"t#bar{t}->e#mu","f");
732     leg5->Draw();
733     CompleteSave(can,"ttbar_emu_mc");
734    
735    
736    
737     }
738    
739     bool is_OF(TCut cut) {
740     string scut = (const char*) cut;
741     if((int)scut.find("id1!=id2")>-1) return true;
742     if((int)scut.find("id1==id2")>-1) return false;
743     return false;
744     }
745    
746     void draw_pure_jzb_histo(TCut cut,string variable,string savename, TCanvas *can,float min,float max,int nbins) {
747     can->cd();
748     can->SetLogy(1);
749     string xlabel="JZB (GeV)";
750    
751     TH1F *datahisto = allsamples.Draw("datahisto",variable,nbins,min,max, xlabel, "events",cut,data,luminosity);
752     THStack mcstack = allsamples.DrawStack("mcstack",variable,nbins,min,max, xlabel, "events",cut,mc,luminosity);
753    
754     datahisto->SetMinimum(0.1);
755 buchmann 1.12 //if(savename=="jzb_OSOF") datahisto->SetMaximum(10);
756 buchmann 1.1 datahisto->Draw("e1");
757     mcstack.Draw("same");
758     datahisto->Draw("same,e1");
759    
760     TLegend *leg;
761     if(is_OF(cut)) leg = allsamples.allbglegend("Opposite Flavor",datahisto);
762     else leg = allsamples.allbglegend("Same Flavor",datahisto);
763     leg->Draw();
764     string write_cut = decipher_cut(cut,"");
765     TText *writeline1 = write_text(0.5,0.965,write_cut.c_str());
766     writeline1->SetTextSize(0.035);
767     writeline1->Draw();
768     CompleteSave(can, ("jzb/"+savename));
769    
770     datahisto->Delete();
771     mcstack.Delete();
772     }
773    
774     Double_t GausR(Double_t *x, Double_t *par) {
775     return gRandom->Gaus(x[0],par[0]);
776     }
777    
778 buchmann 1.11 void jzb_plots(string mcjzb, string datajzb) {
779     TCanvas *can = new TCanvas("can","JZB Plots Canvas");
780 buchmann 1.1 float max=350;
781     float min=-120;
782     int nbins=(max-min)/5.0; // we want 5 GeV/bin
783    
784     stringstream ss;
785     // ss << "GausRandom(" << datajzb << ",0)";
786     draw_pure_jzb_histo(cutOSSF&&cutnJets&&cutmass,datajzb,"jzb_OSSF",can,min,max,nbins);
787     draw_pure_jzb_histo(cutOSOF&&cutnJets&&cutmass,datajzb,"jzb_OSOF",can,min,max,nbins);
788 buchmann 1.12 draw_pure_jzb_histo(cutOSSF&&cutnJets&&sidebandcut,datajzb,"jzb_OSSF_SB",can,min,max,nbins);
789     draw_pure_jzb_histo(cutOSOF&&cutnJets&&sidebandcut,datajzb,"jzb_OSOF_SB",can,min,max,nbins);
790 buchmann 1.1 draw_pure_jzb_histo(cutOSSF&&cutnJets&&cutmass,datajzb,"jzb_OSSF_rebin",can,min,max,nbins/4.0);
791     draw_pure_jzb_histo(cutOSOF&&cutnJets&&cutmass,datajzb,"jzb_OSOF_rebin",can,min,max,nbins/4.0);
792 buchmann 1.12 draw_pure_jzb_histo(cutOSSF&&cutnJets&&sidebandcut,datajzb,"jzb_OSSF_SB_rebin",can,min,max,nbins/4.0);
793     draw_pure_jzb_histo(cutOSOF&&cutnJets&&sidebandcut,datajzb,"jzb_OSOF_SB_rebin",can,min,max,nbins/4.0);
794 buchmann 1.1 }
795    
796     void calculate_all_yields(string mcdrawcommand,vector<float> jzb_cuts) {
797     cout << "Calculating background yields in MC:" << endl;
798     jzb_cuts.push_back(14000);
799 buchmann 1.8 TH1F *allbgs = allsamples.Draw("allbgs",jzbvariablemc,jzb_cuts, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets,mc,luminosity);
800 buchmann 1.1 float cumulative=0;
801     for(int i=allbgs->GetNbinsX();i>=1;i--) {
802     cumulative+=allbgs->GetBinContent(i);
803     cout << "Above " << allbgs->GetBinLowEdge(i) << " GeV/c : " << cumulative << endl;
804     }
805     }
806    
807     void rediscover_the_top(string mcjzb, string datajzb) {
808     cout << "Hi! today we are going to (try to) rediscover the top!" << endl;
809     TCanvas *c3 = new TCanvas("c3","c3");
810     c3->SetLogy(1);
811     vector<float> binning;
812     //binning=allsamples.get_optimal_binsize(mcjzb,cutmass&&cutOSSF&&cutnJets,20,50,800);
813     /*
814     binning.push_back(50);
815     binning.push_back(100);
816     binning.push_back(150);
817     binning.push_back(200);
818     binning.push_back(500);
819    
820    
821     TH1F *dataprediction = allsamples.Draw("dataprediction", "-"+datajzb, binning, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets,data, luminosity);
822     TH1F *puresignal = allsamples.Draw("puresignal", datajzb, binning, "JZB (GeV)", "events", cutmass&&cutOSOF&&cutnJets,data, luminosity);
823     // TH1F *puresignal = allsamples.Draw("puresignal", mcjzb, binning, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets,mc, luminosity,allsamples.FindSample("TTJets"));
824     TH1F *observed = allsamples.Draw("observed", datajzb,binning, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
825     /*
826     ofstream myfile;
827     TH1F *ratio = (TH1F*)observed->Clone();
828     ratio->Divide(dataprediction);
829     ratio->GetYaxis()->SetTitle("Ratio obs/pred");
830     ratio->Draw();
831     c3->SaveAs("testratio.png");
832     myfile.open ("ShapeFit_log.txt");
833     establish_upper_limits(observed,dataprediction,puresignal,"LM4",myfile);
834     myfile.close();
835     */
836    
837    
838 buchmann 1.6 int nbins=100;
839 buchmann 1.1 float low=0;
840 buchmann 1.6 float hi=500;
841 buchmann 1.1 TCanvas *c4 = new TCanvas("c4","c4",900,900);
842     c4->Divide(2,2);
843     c4->cd(1);
844     c4->cd(1)->SetLogy(1);
845     TH1F *datapredictiont = allsamples.Draw("datapredictiont", "-"+datajzb, nbins,low,hi, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets,data, luminosity);
846 buchmann 1.6 TH1F *datapredictiono = allsamples.Draw("datapredictiono", "-"+datajzb, nbins,low,hi, "JZB (GeV)", "events", cutmass&&cutOSOF&&cutnJets,data, luminosity);
847     datapredictiont->Add(datapredictiono,-1);
848 buchmann 1.1 cout << "Second way of doing this !!!! Analytical shape to the left :-D" << endl;
849 buchmann 1.9 vector<TF1*> functions = do_cb_fit_to_plot(datapredictiont,10);
850 buchmann 1.1 datapredictiont->SetMarkerColor(kRed);
851     datapredictiont->SetLineColor(kRed);
852     datapredictiont->Draw();
853 buchmann 1.9 functions[1]->Draw("same");
854 buchmann 1.6 TText *title1 = write_title("Top Background Prediction (JZB<0, with osof subtr)");
855 buchmann 1.1 title1->Draw();
856    
857     c4->cd(2);
858     c4->cd(2)->SetLogy(1);
859     TH1F *observedt = allsamples.Draw("observedt", datajzb, nbins,low,hi, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets,data, luminosity);
860     observedt->Draw();
861     datapredictiont->Draw("histo,same");
862 buchmann 1.9 functions[1]->Draw("same");
863 buchmann 1.6 TText *title2 = write_title("Observed and predicted background");
864 buchmann 1.1 title2->Draw();
865    
866     c4->cd(3);
867     c4->cd(3)->SetLogy(1);
868     // TH1F *ratio = (TH1F*)observedt->Clone();
869 buchmann 1.6
870 buchmann 1.1 TH1F *analytical_background_prediction= new TH1F("analytical_background_prediction","",nbins,low,hi);
871     for(int i=0;i<=nbins;i++) {
872 buchmann 1.9 analytical_background_prediction->SetBinContent(i+1,functions[1]->Eval(((hi-low)/((float)nbins))*(i+0.5)));
873     analytical_background_prediction->SetBinError(i+1,TMath::Sqrt(functions[1]->Eval(((hi-low)/((float)nbins))*(i+0.5))));
874 buchmann 1.1 }
875     analytical_background_prediction->GetYaxis()->SetTitle("JZB (GeV)");
876     analytical_background_prediction->GetYaxis()->CenterTitle();
877     TH1F *analyticaldrawonly = (TH1F*)analytical_background_prediction->Clone();
878     analytical_background_prediction->SetFillColor(TColor::GetColor("#3399FF"));
879     analytical_background_prediction->SetMarkerSize(0);
880     analytical_background_prediction->Draw("e5");
881     analyticaldrawonly->Draw("histo,same");
882 buchmann 1.9 functions[1]->Draw("same");
883 buchmann 1.1 TText *title3 = write_title("Analytical bg pred histo");
884     title3->Draw();
885    
886     c4->cd(4);
887     // c4->cd(4)->SetLogy(1);
888     vector<float> ratio_binning;
889     ratio_binning.push_back(0);
890     ratio_binning.push_back(5);
891     ratio_binning.push_back(10);
892     ratio_binning.push_back(20);
893     ratio_binning.push_back(50);
894 buchmann 1.6 // ratio_binning.push_back(60);
895     /*
896     ratio_binning.push_back(51);
897     ratio_binning.push_back(52);
898     ratio_binning.push_back(53);
899     ratio_binning.push_back(54);
900     ratio_binning.push_back(55);
901     ratio_binning.push_back(56);
902     ratio_binning.push_back(57);
903     ratio_binning.push_back(58);
904     ratio_binning.push_back(59);
905     ratio_binning.push_back(60);
906     // ratio_binning.push_back(70);*/
907     // ratio_binning.push_back(80);
908     // ratio_binning.push_back(90);
909     ratio_binning.push_back(80);
910     // ratio_binning.push_back(110);
911 buchmann 1.1 ratio_binning.push_back(500);
912    
913     TH1F *observedtb = allsamples.Draw("observedtb", datajzb, ratio_binning, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets,data, luminosity);
914     TH1F *datapredictiontb = allsamples.Draw("datapredictiontb", "-"+datajzb, ratio_binning, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets,data, luminosity);
915 buchmann 1.6 TH1F *datapredictiontbo = allsamples.Draw("datapredictiontbo", "-"+datajzb, ratio_binning, "JZB (GeV)", "events", cutmass&&cutOSOF&&cutnJets,data, luminosity);
916     datapredictiontb->Add(datapredictiontbo,-1);
917     TH1F *analytical_background_predictionb = allsamples.Draw("analytical_background_predictionb",datajzb, ratio_binning, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets&&"mll<2",data, luminosity);
918 buchmann 1.1 for(int i=0;i<=ratio_binning.size();i++) {
919 buchmann 1.9 analytical_background_predictionb->SetBinContent(i+1,functions[1]->Eval(analytical_background_predictionb->GetBinCenter(i)));
920     analytical_background_predictionb->SetBinError(i+1,TMath::Sqrt(functions[1]->Eval(analytical_background_predictionb->GetBinCenter(i))));
921 buchmann 1.1 }
922    
923     TH1F *ratio = (TH1F*) observedtb->Clone();
924     ratio->Divide(datapredictiontb);
925 buchmann 1.6
926     for (int i=0;i<=ratio_binning.size();i++) {
927     cout << observedtb->GetBinLowEdge(i+1) << ";"<<observedtb->GetBinContent(i+1) << ";" << datapredictiontb->GetBinContent(i+1) << " --> " << ratio->GetBinContent(i+1) << "+/-" << ratio->GetBinError(i+1) << endl;
928     }
929    
930     // ratio->Divide(datapredictiontb);
931 buchmann 1.1 // ratio->Divide(analytical_background_predictionb);
932     // TGraphAsymmErrors *JZBratio= histRatio(observedtb,analytical_background_predictionb,data,ratio_binning);
933     // ratio->Divide(analytical_background_prediction);
934     // ratio->Divide(datapredictiont);
935     // ratio->GetYaxis()->SetTitle("obs/pred");
936     // JZBratio->Draw("AP");
937     ratio->GetYaxis()->SetRangeUser(0,10);
938     ratio->Draw();
939     //analytical_background_predictionb->Draw();
940     // JZBratio->SetTitle("");
941     TText *title4 = write_title("Ratio of observed to predicted");
942     title4->Draw();
943    
944     // CompleteSave(c4,"test/ttbar_discovery_dataprediction___analytical_function");
945     CompleteSave(c4,"test/ttbar_discovery_dataprediction__analytical__new_binning_one_huge_bin_from_80");
946    
947    
948    
949    
950    
951     }
952    
953     void calculate_upper_limits(string mcjzb, string datajzb) {
954 buchmann 1.11 write_warning("calculate_upper_limits","Upper limit calculation temporarily deactivated");
955     // write_warning("calculate_upper_limits","Calculation of SUSY upper limits has been temporarily suspended in favor of top discovery");
956     // rediscover_the_top(mcjzb,datajzb);
957 buchmann 1.6 /*
958 buchmann 1.1 TCanvas *c3 = new TCanvas("c3","c3");
959     c3->SetLogy(1);
960     vector<float> binning;
961     //binning=allsamples.get_optimal_binsize(mcjzb,cutmass&&cutOSSF&&cutnJets,20,50,800);
962     binning.push_back(50);
963     binning.push_back(100);
964     binning.push_back(150);
965     binning.push_back(200);
966     binning.push_back(500);
967     TH1F *datapredictiona = allsamples.Draw("datapredictiona", "-"+datajzb, binning, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets,mc, luminosity);
968     TH1F *datapredictionb = allsamples.Draw("datapredictionb", "-"+datajzb, binning, "JZB (GeV)", "events", cutmass&&cutOSOF&&cutnJets,mc, luminosity);
969     TH1F *datapredictionc = allsamples.Draw("datapredictionc", datajzb, binning, "JZB (GeV)", "events", cutmass&&cutOSOF&&cutnJets,mc, luminosity);
970 buchmann 1.5 TH1F *dataprediction = (TH1F*)datapredictiona->Clone();
971 buchmann 1.1 dataprediction->Add(datapredictionb,-1);
972     dataprediction->Add(datapredictionc);
973     TH1F *puresignal = allsamples.Draw("puresignal", mcjzb, binning, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets,mc, luminosity,allsamples.FindSample("LM4"));
974 buchmann 1.5 TH1F *signalpred = allsamples.Draw("signalpred", "-"+mcjzb, binning, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets,mc, luminosity,allsamples.FindSample("LM4"));
975     TH1F *signalpredlo = allsamples.Draw("signalpredlo", "-"+mcjzb, binning, "JZB (GeV)", "events", cutmass&&cutOSOF&&cutnJets,mc, luminosity,allsamples.FindSample("LM4"));
976     TH1F *signalpredro = allsamples.Draw("signalpredro", mcjzb, binning, "JZB (GeV)", "events", cutmass&&cutOSOF&&cutnJets,mc, luminosity,allsamples.FindSample("LM4"));
977 buchmann 1.1 TH1F *puredata = allsamples.Draw("puredata", datajzb,binning, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
978 buchmann 1.5 signalpred->Add(signalpredlo,-1);
979     signalpred->Add(signalpredro);
980     puresignal->Add(signalpred,-1);//subtracting signal contamination
981 buchmann 1.1 ofstream myfile;
982     myfile.open ("ShapeFit_log.txt");
983     establish_upper_limits(puredata,dataprediction,puresignal,"LM4",myfile);
984     myfile.close();
985 buchmann 1.6 */
986 buchmann 1.1 }
987    
988     void susy_scan_axis_labeling(TH2F *histo) {
989     histo->GetXaxis()->SetTitle("#Chi_{2}^{0}-LSP");
990     histo->GetXaxis()->CenterTitle();
991     histo->GetYaxis()->SetTitle("m_{#tilde{q}}");
992     histo->GetYaxis()->CenterTitle();
993     }
994    
995     void scan_susy_space(string mcjzb, string datajzb) {
996     TCanvas *c3 = new TCanvas("c3","c3");
997     vector<float> binning;
998     binning=allsamples.get_optimal_binsize(mcjzb,cutmass&&cutOSSF&&cutnJets,20,50,800);
999     /*
1000     binning.push_back(50);
1001     binning.push_back(75);
1002     binning.push_back(100);
1003     binning.push_back(150);
1004     binning.push_back(200);
1005     binning.push_back(500);
1006     */
1007     float arrbinning[binning.size()];
1008     for(int i=0;i<binning.size();i++) arrbinning[i]=binning[i];
1009     TH1F *puredata = allsamples.Draw("puredata", datajzb,binning, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
1010 buchmann 1.2 TH1F *allbgs = allsamples.Draw("allbgs", "-"+datajzb,binning, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
1011     TH1F *allbgsb = allsamples.Draw("allbgsb", "-"+datajzb,binning, "JZB (GeV)", "events", cutmass&&cutOSOF&&cutnJets,data,luminosity);
1012     TH1F *allbgsc = allsamples.Draw("allbgsc", datajzb,binning, "JZB (GeV)", "events", cutmass&&cutOSOF&&cutnJets,data,luminosity);
1013     allbgs->Add(allbgsb,-1);
1014     allbgs->Add(allbgsc);
1015 buchmann 1.1 int ndata=puredata->Integral();
1016     ofstream myfile;
1017     myfile.open ("susyscan_log.txt");
1018     TFile *susyscanfile = new TFile("/scratch/fronga/SMS/T5z_SqSqToQZQZ_38xFall10.root");
1019     TTree *suevents = (TTree*)susyscanfile->Get("events");
1020     TH2F *exclusionmap = new TH2F("exclusionmap","",20,0,500,20,0,1000);
1021     TH2F *exclusionmap1s = new TH2F("exclusionmap1s","",20,0,500,20,0,1000);
1022     TH2F *exclusionmap2s = new TH2F("exclusionmap2s","",20,0,500,20,0,1000);
1023     TH2F *exclusionmap3s = new TH2F("exclusionmap3s","",20,0,500,20,0,1000);
1024    
1025     susy_scan_axis_labeling(exclusionmap);
1026     susy_scan_axis_labeling(exclusionmap1s);
1027     susy_scan_axis_labeling(exclusionmap2s);
1028     susy_scan_axis_labeling(exclusionmap3s);
1029    
1030     Int_t MyPalette[100];
1031     Double_t r[] = {0., 0.0, 1.0, 1.0, 1.0};
1032     Double_t g[] = {0., 0.0, 0.0, 1.0, 1.0};
1033     Double_t b[] = {0., 1.0, 0.0, 0.0, 1.0};
1034     Double_t stop[] = {0., .25, .50, .75, 1.0};
1035     Int_t FI = TColor::CreateGradientColorTable(5, stop, r, g, b, 100);
1036     for (int i=0;i<100;i++) MyPalette[i] = FI+i;
1037    
1038     gStyle->SetPalette(100, MyPalette);
1039    
1040     for(int m23=50;m23<500;m23+=25) {
1041     for (int m0=(2*(m23-50)+150);m0<=1000;m0+=50)
1042     {
1043     c3->cd();
1044     stringstream drawcondition;
1045     drawcondition << "pfJetGoodNum>=3&&(TMath::Abs(masses[0]-"<<m0<<")<10&&TMath::Abs(masses[2]-masses[3]-"<<m23<<")<10)&&mll>5&&id1==id2";
1046     TH1F *puresignal = new TH1F("puresignal","puresignal",binning.size()-1,arrbinning);
1047 buchmann 1.4 TH1F *puresignall= new TH1F("puresignall","puresignal",binning.size()-1,arrbinning);
1048     stringstream drawvar,drawvar2;
1049 buchmann 1.1 drawvar<<mcjzb<<">>puresignal";
1050 buchmann 1.4 drawvar2<<"-"<<mcjzb<<">>puresignall";
1051 buchmann 1.1 suevents->Draw(drawvar.str().c_str(),drawcondition.str().c_str());
1052 buchmann 1.4 suevents->Draw(drawvar2.str().c_str(),drawcondition.str().c_str());
1053 buchmann 1.1 if(puresignal->Integral()<60) {
1054     delete puresignal;
1055     continue;
1056     }
1057 buchmann 1.4 puresignal->Add(puresignall,-1);//we need to correct for the signal contamination - we effectively only see (JZB>0)-(JZB<0) !!
1058 buchmann 1.1 puresignal->Scale(ndata/(20*puresignal->Integral()));//normalizing it to 5% of the data
1059     stringstream saveas;
1060     saveas<<"Model_Scan/m0_"<<m0<<"__m23_"<<m23;
1061 buchmann 1.5 cout << "PLEASE KEEP IN MIND THAT SIGNAL CONTAMINATION IS NOT REALLY TAKEN CARE OF YET DUE TO LOW STATISTICS! SHOULD BE SOMETHING LIKE THIS : "<< endl;
1062     // TH1F *signalpredlo = allsamples.Draw("signalpredlo", "-"+mcjzb, binning, "JZB (GeV)", "events", cutmass&&cutOSOF&&cutnJets,mc, luminosity,allsamples.FindSample("LM4"));
1063     // TH1F *signalpredro = allsamples.Draw("signalpredro", mcjzb, binning, "JZB (GeV)", "events", cutmass&&cutOSOF&&cutnJets,mc, luminosity,allsamples.FindSample("LM4"));
1064     // TH1F *puredata = allsamples.Draw("puredata", datajzb,binning, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
1065     // signalpred->Add(signalpredlo,-1);
1066     // signalpred->Add(signalpredro);
1067     // puresignal->Add(signalpred,-1);//subtracting signal contamination
1068     //---------------------
1069 buchmann 1.1 // cout << "(m0,m23)=("<<m0<<","<<m23<<") contains " << puresignal->Integral() << endl;
1070     // TH1F *puresignal = allsamples.Draw("puresignal",mcjzb, binning, "JZB (GeV)", "events", cutmass&&cutOSSF&&cutnJets,mc, luminosity,allsamples.FindSample("LM4"));
1071     vector<float> results=establish_upper_limits(puredata,allbgs,puresignal,saveas.str(),myfile);
1072     if(results.size()==0) {
1073     delete puresignal;
1074     continue;
1075     }
1076     exclusionmap->Fill(m23,m0,results[0]);
1077     exclusionmap1s->Fill(m23,m0,results[1]);
1078     exclusionmap2s->Fill(m23,m0,results[2]);
1079     exclusionmap3s->Fill(m23,m0,results[3]);
1080     delete puresignal;
1081     cout << "(m0,m23)=("<<m0<<","<<m23<<") : 3 sigma at " << results[3] << endl;
1082     }
1083     }//end of model scan for loop
1084    
1085     cout << "Exclusion Map contains" << exclusionmap->Integral() << " (integral) and entries: " << exclusionmap->GetEntries() << endl;
1086     c3->cd();
1087     exclusionmap->Draw("CONTZ");
1088     CompleteSave(c3,"Model_Scan/CONT/Model_Scan_Mean_values");
1089     exclusionmap->Draw("COLZ");
1090     CompleteSave(c3,"Model_Scan/COL/Model_Scan_Mean_values");
1091    
1092     exclusionmap1s->Draw("CONTZ");
1093     CompleteSave(c3,"Model_Scan/CONT/Model_Scan_1sigma_values");
1094     exclusionmap1s->Draw("COLZ");
1095     CompleteSave(c3,"Model_Scan/COL/Model_Scan_1sigma_values");
1096    
1097     exclusionmap2s->Draw("CONTZ");
1098     CompleteSave(c3,"Model_Scan/CONT/Model_Scan_2sigma_values");
1099     exclusionmap2s->Draw("COLZ");
1100     CompleteSave(c3,"Model_Scan/COL/Model_Scan_2sigma_values");
1101    
1102     exclusionmap3s->Draw("CONTZ");
1103     CompleteSave(c3,"Model_Scan/CONT/Model_Scan_3sigma_values");
1104     exclusionmap3s->Draw("COLZ");
1105     CompleteSave(c3,"Model_Scan/COL/Model_Scan_3sigma_values");
1106    
1107     TFile *exclusion_limits = new TFile("exclusion_limits.root","RECREATE");
1108     exclusionmap->Write();
1109     exclusionmap1s->Write();
1110     exclusionmap2s->Write();
1111     exclusionmap3s->Write();
1112     exclusion_limits->Close();
1113     susyscanfile->Close();
1114    
1115     myfile.close();
1116     }
1117 buchmann 1.2
1118 buchmann 1.4 void draw_ttbar_and_zjets_shape_for_one_configuration(string mcjzb, string datajzb, int leptontype=-1, int scenario=0,bool floating=false) {
1119 buchmann 1.10 //Step 1: Establishing cuts
1120 buchmann 1.3 stringstream jetcutstring;
1121 buchmann 1.10 string writescenario="";
1122    
1123 buchmann 1.3 if(scenario==0) jetcutstring << "(pfJetGoodNum>=3)&&"<<(const char*) basicqualitycut;
1124     if(scenario==1) jetcutstring << "(pfJetPt[0]>50&&pfJetPt[1]>50)&&"<<(const char*)basicqualitycut;
1125     TCut jetcut(jetcutstring.str().c_str());
1126     string leptoncut="mll>0";
1127     if(leptontype==0||leptontype==1) {
1128 buchmann 1.10 if(leptontype==0) {
1129     leptoncut="id1==0";
1130     writescenario="__ee";
1131     }
1132     else {
1133     leptoncut="id1==1";
1134     writescenario="__ee";
1135     }
1136 buchmann 1.3 }
1137     TCut lepcut(leptoncut.c_str());
1138    
1139 buchmann 1.6 TCanvas *c5 = new TCanvas("c5","c5",1500,500);
1140 buchmann 1.10 TCanvas *c6 = new TCanvas("c6","c6");
1141 buchmann 1.3 c5->Divide(3,1);
1142 buchmann 1.10
1143     //STEP 2: Extract Zjets shape in data
1144     c5->cd(1);
1145 buchmann 1.3 c5->cd(1)->SetLogy(1);
1146     TCut massat40("mll>40");
1147     TH1F *ossfleft = allsamples.Draw("ossfleft", "-"+datajzb,40,0,200, "JZB (GeV)", "events", massat40&&cutOSSF&&jetcut&&lepcut,data,luminosity);
1148     TH1F *osofleft = allsamples.Draw("osofleft", "-"+datajzb,40,0,200, "JZB (GeV)", "events", massat40&&cutOSOF&&jetcut&&lepcut,data,luminosity);
1149     ossfleft->SetLineColor(kRed);
1150     ossfleft->SetMarkerColor(kRed);
1151     ossfleft->Add(osofleft,-1);
1152 buchmann 1.9 vector<TF1*> functions = do_cb_fit_to_plot(ossfleft,10);
1153 buchmann 1.3 ossfleft->Draw();
1154 buchmann 1.9 functions[0]->Draw("same");functions[1]->Draw("same");functions[2]->Draw("same");
1155     TF1 *zjetsfunc = (TF1*) functions[1]->Clone();
1156     TF1 *zjetsfuncN = (TF1*) functions[0]->Clone();
1157     TF1 *zjetsfuncP = (TF1*) functions[2]->Clone();
1158 buchmann 1.6 zjetsfunc->Draw("same");zjetsfuncN->Draw("same");zjetsfuncP->Draw("same");
1159 buchmann 1.3 TLegend *leg1 = new TLegend(0.6,0.6,0.89,0.80);
1160     leg1->SetFillColor(kWhite);
1161     leg1->SetLineColor(kWhite);
1162     leg1->AddEntry(ossfleft,"OSSF (sub),JZB<peak","p");
1163 buchmann 1.6 leg1->AddEntry(zjetsfunc,"OSSF fit ('zjets')","l");
1164 buchmann 1.3 leg1->Draw("same");
1165     TText *titleleft = write_title("Extracting Z+Jets shape");
1166     titleleft->Draw();
1167    
1168 buchmann 1.10 //Step 3: Extract ttbar shape (in data or MC?)
1169     c5->cd(2);
1170 buchmann 1.3 c5->cd(2)->SetLogy(1);
1171 buchmann 1.4 TH1F *osof;
1172     TText *titlecenter;
1173 buchmann 1.10 bool frommc=false;
1174     if(frommc) {
1175 buchmann 1.4 osof = allsamples.Draw("osof",datajzb,40,-200,200, "JZB (GeV)", "events", massat40&&cutOSSF&&jetcut&&lepcut,mc,luminosity,allsamples.FindSample("TTJets"));
1176     titlecenter = write_title("Extracting ttbar shape (from ossf MC)");
1177     }
1178     else {
1179     osof = allsamples.Draw("osof",datajzb,40,-200,200, "JZB (GeV)", "events", massat40&&cutOSOF&&jetcut&&lepcut,data,luminosity);
1180     titlecenter = write_title("Extracting ttbar shape (from osof data)");
1181     }
1182 buchmann 1.3 osof->Draw();
1183 buchmann 1.10 vector<TF1*> ttbarfunctions = do_cb_fit_to_plot(osof,35,true);
1184     ttbarfunctions[0]->SetLineColor(kRed); ttbarfunctions[0]->SetLineStyle(2); ttbarfunctions[0]->Draw("same");
1185     ttbarfunctions[1]->SetLineColor(kRed); ttbarfunctions[1]->Draw("same");
1186     ttbarfunctions[2]->SetLineColor(kRed); ttbarfunctions[2]->SetLineStyle(2); ttbarfunctions[2]->Draw("same");
1187    
1188 buchmann 1.3 TLegend *leg2 = new TLegend(0.15,0.8,0.4,0.89);
1189     leg2->SetFillColor(kWhite);
1190     leg2->SetLineColor(kWhite);
1191 buchmann 1.10 if(frommc) {
1192 buchmann 1.4 leg2->AddEntry(osof,"t#bar{t} OSSF, MC","p");
1193 buchmann 1.10 leg2->AddEntry(ttbarfunctions[1],"Fit to t#bar{t} OSSF,MC","l");
1194 buchmann 1.4 } else {
1195     leg2->AddEntry(osof,"OSOF","p");
1196 buchmann 1.10 leg2->AddEntry(ttbarfunctions[1],"Fit to OSOF","l");
1197 buchmann 1.4 }
1198 buchmann 1.3 leg2->Draw("same");
1199     titlecenter->Draw();
1200 buchmann 1.10
1201     //--------------------------------------------------------------------------------------------------------------------------------
1202 buchmann 1.4 //STEP 4: Present it!
1203     // actually: if we wanna let it float we need to do that first :-)
1204 buchmann 1.3 c5->cd(3);
1205     c5->cd(3)->SetLogy(1);
1206     TH1F *observed = allsamples.Draw("observed",datajzb,100,0,500, "JZB (GeV)", "events", massat40&&cutOSSF&&jetcut&&lepcut,data,luminosity);
1207 buchmann 1.4
1208 buchmann 1.6 TF1 *logparc = new TF1("logparc",InvCrystalBall,0,1000,5); logparc->SetLineColor(kRed);
1209 buchmann 1.8 TF1 *logparcn = new TF1("logparcn",InvCrystalBallN,0,1000,5); logparcn->SetLineColor(kRed); logparcn->SetLineStyle(2);
1210     TF1 *logparcp = new TF1("logparcp",InvCrystalBallP,0,1000,5); logparcp->SetLineColor(kRed); logparcp->SetLineStyle(2);
1211 buchmann 1.6
1212     TF1 *zjetsc = new TF1("zjetsc",InvCrystalBall,0,1000,5); zjetsc->SetLineColor(kBlue);
1213     TF1 *zjetscn = new TF1("zjetscn",InvCrystalBallN,0,1000,5); zjetscn->SetLineColor(kBlue); zjetscn->SetLineStyle(2);
1214     TF1 *zjetscp = new TF1("zjetscp",InvCrystalBallP,0,1000,5); zjetscp->SetLineColor(kBlue); zjetscp->SetLineStyle(2);
1215    
1216 buchmann 1.10 TF1 *ZplusJetsplusTTbar = new TF1("ZplusJetsplusTTbar", DoubleInvCrystalBall,0,1000,10); ZplusJetsplusTTbar->SetLineColor(kBlue);
1217     TF1 *ZplusJetsplusTTbarP= new TF1("ZplusJetsplusTTbarP",DoubleInvCrystalBallP,0,1000,10); ZplusJetsplusTTbarP->SetLineColor(kBlue); ZplusJetsplusTTbarP->SetLineStyle(2);
1218     TF1 *ZplusJetsplusTTbarN= new TF1("ZplusJetsplusTTbarN",DoubleInvCrystalBallN,0,1000,10); ZplusJetsplusTTbarN->SetLineColor(kBlue); ZplusJetsplusTTbarN->SetLineStyle(2);
1219    
1220 buchmann 1.6 zjetsc->SetParameters(zjetsfunc->GetParameters());
1221     zjetscp->SetParameters(zjetsfunc->GetParameters());
1222     zjetscn->SetParameters(zjetsfunc->GetParameters());
1223 buchmann 1.10
1224 buchmann 1.8 TH1F *observeda = allsamples.Draw("observeda",datajzb,53,80,350, "JZB (GeV)", "events", massat40&&cutOSSF&&jetcut&&lepcut,data,luminosity);
1225 buchmann 1.7 //blublu
1226 buchmann 1.10 logparc->SetParameters(ttbarfunctions[1]->GetParameters());
1227     logparcn->SetParameters(ttbarfunctions[1]->GetParameters());
1228     logparcp->SetParameters(ttbarfunctions[1]->GetParameters());
1229 buchmann 1.4 if(floating) {
1230 buchmann 1.10 cout << "TTbar contribution assumed (before fitting) : " << logparc->GetParameter(0) << endl;
1231     logparc->SetParameters(ttbarfunctions[1]->GetParameters());
1232 buchmann 1.7 for(int i=0;i<10;i++) {
1233     if(i<5) ZplusJetsplusTTbar->FixParameter(i,zjetsfunc->GetParameter(i));
1234     if(i>=5) {
1235     if (i>5) ZplusJetsplusTTbar->FixParameter(i,logparc->GetParameter(i-5));
1236     if (i==5) ZplusJetsplusTTbar->SetParameter(i,logparc->GetParameter(i-5));
1237     }
1238     }//end of setting parameters
1239     observeda->Draw("same");
1240     ZplusJetsplusTTbar->Draw("same");
1241     observeda->Fit(ZplusJetsplusTTbar);
1242     cout << "--> Quality of Z+Jets / TTbar fit : chi2/ndf = " << ZplusJetsplusTTbar->GetChisquare() << "/" << ZplusJetsplusTTbar->GetNDF() << endl;
1243     ZplusJetsplusTTbar->Draw("same");
1244 buchmann 1.8 ZplusJetsplusTTbarP->SetParameters(ZplusJetsplusTTbar->GetParameters());
1245     ZplusJetsplusTTbarN->SetParameters(ZplusJetsplusTTbar->GetParameters());
1246 buchmann 1.10 cout << "TTbar contribution found (after fitting) : " << ZplusJetsplusTTbar->GetParameter(5) << endl;
1247 buchmann 1.7 float factor = ZplusJetsplusTTbar->GetParameter(5) / logparc->GetParameter(0);
1248     cout << "FACTOR: " << factor << endl;
1249 buchmann 1.10 logparc->SetParameter(0,factor*ttbarfunctions[1]->GetParameter(0));
1250     logparcn->SetParameter(0,factor*ttbarfunctions[1]->GetParameter(0));
1251     logparcp->SetParameter(0,factor*ttbarfunctions[1]->GetParameter(0));
1252 buchmann 1.4 }
1253 buchmann 1.6
1254 buchmann 1.4 c5->cd(3);
1255 buchmann 1.6 c5->cd(3)->SetLogy(1);
1256     observed->Draw();
1257     zjetsc->Draw("same");zjetscn->Draw("same");zjetscp->Draw("same");
1258 buchmann 1.3 logparc->Draw("same");
1259 buchmann 1.6 logparcn->Draw("same");
1260     logparcp->Draw("same");
1261    
1262 buchmann 1.3 TLegend *leg3 = new TLegend(0.6,0.6,0.89,0.80);
1263     leg3->SetFillColor(kWhite);
1264     leg3->SetLineColor(kWhite);
1265     leg3->AddEntry(observed,"OSSF,JZB>peak","p");
1266 buchmann 1.10 leg3->AddEntry(ttbarfunctions[1],"OSOF fit ('ttbar')","l");
1267 buchmann 1.6 leg3->AddEntry(zjetsfunc,"OSSF,JZB<0 fit ('zjets')","l");
1268 buchmann 1.3 leg3->Draw("same");
1269     TText *titleright = write_title("Summary of shapes and observed shape");
1270     titleright->Draw();
1271    
1272     c6->cd()->SetLogy(1);
1273 buchmann 1.6 observed->Draw();
1274     zjetsc->Draw("same");zjetscn->Draw("same");zjetscp->Draw("same");
1275 buchmann 1.3 logparc->Draw("same");
1276 buchmann 1.6 logparcn->Draw("same");
1277     logparcp->Draw("same");
1278 buchmann 1.3 leg3->Draw("same");
1279     titleright->Draw();
1280    
1281     if(scenario==0) {
1282     CompleteSave(c5,"Shapes2/Making_of___3jetsabove30"+writescenario);
1283 buchmann 1.6 CompleteSave(c5->cd(1),"Shapes2/Making_of___3jetsabove30"+writescenario+"__cd1");
1284     CompleteSave(c5->cd(2),"Shapes2/Making_of___3jetsabove30"+writescenario+"__cd2");
1285     CompleteSave(c5->cd(3),"Shapes2/Making_of___3jetsabove30"+writescenario+"__cd3");
1286 buchmann 1.3 CompleteSave(c6,"Shapes2/Background_Shapes___3jetsabove30"+writescenario);
1287     } else {
1288     CompleteSave(c5,"Shapes2/Making_of___2jetsabove50"+writescenario);
1289 buchmann 1.6 CompleteSave(c5->cd(1),"Shapes2/Making_of___2jetsabove50"+writescenario+"__cd1");
1290     CompleteSave(c5->cd(2),"Shapes2/Making_of___2jetsabove50"+writescenario+"__cd2");
1291     CompleteSave(c5->cd(3),"Shapes2/Making_of___2jetsabove50"+writescenario+"__cd3");
1292 buchmann 1.3 CompleteSave(c6,"Shapes2/Background_Shapes___2jetsabove50"+writescenario);
1293     }
1294 buchmann 1.6 cout << "Statistics about our fits: " << endl;
1295 buchmann 1.7 cout << "Z+Jets shape: Chi2/ndf = " << zjetsfunc->GetChisquare() << "/" << ossfleft->GetNbinsX() << endl;
1296 buchmann 1.10 cout << "ttbar shape: Chi2/ndf = " << ttbarfunctions[1]->GetChisquare() << "/" << osof->GetNbinsX() << endl;
1297 buchmann 1.4
1298 buchmann 1.7 c6->cd();
1299 buchmann 1.8 TLegend *additionallegend = new TLegend(0.6,0.6,0.89,0.89);
1300     additionallegend->SetFillColor(kWhite);
1301     additionallegend->SetLineColor(kWhite);
1302     additionallegend->AddEntry(observed,"Data","p");
1303     additionallegend->AddEntry(ZplusJetsplusTTbar,"Fitted Z+jets & TTbar","l");
1304     additionallegend->AddEntry(zjetsc,"Z+jets","l");
1305     additionallegend->AddEntry(logparc,"TTbar","l");
1306 buchmann 1.7 observed->Draw();
1307     ZplusJetsplusTTbar->SetLineColor(kGreen);
1308 buchmann 1.8 ZplusJetsplusTTbarP->SetLineColor(kGreen);
1309     ZplusJetsplusTTbarN->SetLineColor(kGreen);
1310     ZplusJetsplusTTbarP->SetLineStyle(2);
1311     ZplusJetsplusTTbarN->SetLineStyle(2);
1312     TF1 *ZplusJetsplusTTbar2 = new TF1("ZplusJetsplusTTbar2",DoubleInvCrystalBall,0,1000,10);
1313     ZplusJetsplusTTbar2->SetParameters(ZplusJetsplusTTbar->GetParameters());
1314     ZplusJetsplusTTbar2->SetLineColor(kGreen);
1315     ZplusJetsplusTTbarP->SetFillColor(TColor::GetColor("#81F781"));
1316     ZplusJetsplusTTbarN->SetFillColor(kWhite);
1317     ZplusJetsplusTTbarP->Draw("fcsame");
1318     ZplusJetsplusTTbarN->Draw("fcsame");
1319     TH1F *hZplusJetsplusTTbar = (TH1F*)ZplusJetsplusTTbar2->GetHistogram();
1320     TH1F *hZplusJetsplusTTbarN = (TH1F*)ZplusJetsplusTTbarN->GetHistogram();
1321     TH1F *hZplusJetsplusTTbarP = (TH1F*)ZplusJetsplusTTbarP->GetHistogram();
1322     hZplusJetsplusTTbar->SetMarkerSize(0);
1323     hZplusJetsplusTTbarP->SetMarkerSize(0);
1324     hZplusJetsplusTTbarN->SetMarkerSize(0);
1325     for (int i=1;i<=hZplusJetsplusTTbar->GetNbinsX();i++) {
1326     float newerror=hZplusJetsplusTTbarP->GetBinContent(i)-hZplusJetsplusTTbar->GetBinContent(i);
1327     hZplusJetsplusTTbar->SetBinError(i,newerror);
1328     if(hZplusJetsplusTTbar->GetBinContent(i)<0.05) hZplusJetsplusTTbar->SetBinContent(i,0); //avoiding a displaying probolem
1329     }
1330     hZplusJetsplusTTbarP->SetFillColor(kGreen);
1331     hZplusJetsplusTTbarN->SetFillColor(kWhite);
1332     hZplusJetsplusTTbarN->Draw("same");
1333    
1334     ZplusJetsplusTTbar2->SetMarkerSize(0);
1335     ZplusJetsplusTTbar2->Draw("same");
1336    
1337 buchmann 1.7 zjetsc->Draw("same");zjetscn->Draw("same");zjetscp->Draw("same");
1338     logparc->Draw("same");
1339     logparcn->Draw("same");
1340     logparcp->Draw("same");
1341 buchmann 1.8 additionallegend->Draw("same");
1342     if(scenario==0) {
1343     CompleteSave(c6,"Shapes2/Background_Shapes___3jetsabove30__allfits__"+writescenario);
1344     } else {
1345     CompleteSave(c6,"Shapes2/Background_Shapes___2jetsabove50__allfits__"+writescenario);
1346     }
1347 buchmann 1.10 //--------------------------------------------------------------------------------------------------------------------------------
1348 buchmann 1.3 }
1349    
1350 buchmann 1.9
1351    
1352 buchmann 1.2 void draw_ttbar_and_zjets_shape(string mcjzb, string datajzb) {
1353 buchmann 1.3 int all_leptons=-1;
1354     int electrons_only=0;
1355     int mu_only=1;
1356     int twojetswith50gev=1;
1357     int threejetswith30gev=0;
1358 buchmann 1.4 /*
1359 buchmann 1.3 draw_ttbar_and_zjets_shape_for_one_configuration(mcjzb,datajzb,all_leptons,twojetswith50gev);
1360     draw_ttbar_and_zjets_shape_for_one_configuration(mcjzb,datajzb,all_leptons,threejetswith30gev);
1361    
1362     draw_ttbar_and_zjets_shape_for_one_configuration(mcjzb,datajzb,electrons_only,twojetswith50gev);
1363     draw_ttbar_and_zjets_shape_for_one_configuration(mcjzb,datajzb,electrons_only,threejetswith30gev);
1364    
1365     draw_ttbar_and_zjets_shape_for_one_configuration(mcjzb,datajzb,mu_only,twojetswith50gev);
1366     draw_ttbar_and_zjets_shape_for_one_configuration(mcjzb,datajzb,mu_only,threejetswith30gev);
1367 buchmann 1.4 */
1368 buchmann 1.2
1369 buchmann 1.10 draw_ttbar_and_zjets_shape_for_one_configuration(mcjzb,datajzb,all_leptons,threejetswith30gev,true);
1370 buchmann 1.2 }
1371 buchmann 1.6
1372 buchmann 1.8 void find_correction_factors(string &jzbvardata,string &jzbvarmc) {
1373 buchmann 1.9 TCut zptcutforresponse("pt>30&&pt<300&&TMath::Abs(91.2-mll)<20&&id1==id2&&(ch1*ch2<0)");
1374 buchmann 1.8 TH1F *dresponse = allsamples.Draw("dresponse","sumJetPt[1]/pt",100,0,5, "JZB (GeV)", "events", zptcutforresponse,data,luminosity);
1375     float datacorrection=dresponse->GetMean();
1376     stringstream jzbvardatas;
1377 buchmann 1.9 if(datacorrection>1) jzbvardatas<<"(jzb[1]-"<<datacorrection-1<<"*pt)";
1378     if(datacorrection<1) jzbvardatas<<"(jzb[1]+"<<1-datacorrection<<"*pt)";
1379 buchmann 1.8 jzbvardata=jzbvardatas.str();
1380 buchmann 1.11 TH1F *mcresponse = allsamples.Draw("mcresponse","sumJetPt[1]/pt",100,0,5, "JZB (GeV)", "events", zptcutforresponse,mc,luminosity,allsamples.FindSample("DYJetsToLL_TuneZ2_M"));
1381 buchmann 1.8 float mccorrection=mcresponse->GetMean();
1382     stringstream jzbvarmcs;
1383 buchmann 1.9 if(mccorrection>1) jzbvarmcs<<"(jzb[1]-"<<mccorrection-1<<"*pt)";
1384     if(mccorrection<1) jzbvarmcs<<"(jzb[1]+"<<1-mccorrection<<"*pt)";
1385 buchmann 1.8 jzbvarmc=jzbvarmcs.str();
1386     cout << "JZB Z pt correction summary : " << endl;
1387     cout << " Data: The response is " << dresponse->GetMean() << " --> jzb variable is now : " << jzbvardata << endl;
1388     cout << " MC : The response is " << mcresponse->GetMean() << " --> jzb variable is now : " << jzbvarmc << endl;
1389     }
1390 buchmann 1.11
1391 buchmann 1.12 void do_experimental_pred_obs_calculation(float cut ,string mcjzb,string datajzb, int mcordata) {
1392 buchmann 1.11 cout << "Crunching the numbers for JZB>" << cut << endl;
1393 buchmann 1.12 string xlabel="JZB (GeV) -- for algoritm internal use only!";
1394     TH1F *ZOSSFP = allsamples.Draw("ZOSSFP",datajzb,1,cut,14000, xlabel, "events",cutmass&&cutOSSF&&cutnJets&&basiccut,mcordata,luminosity);
1395     TH1F *ZOSOFP = allsamples.Draw("ZOSOFP",datajzb,1,cut,14000, xlabel, "events",cutmass&&cutOSOF&&cutnJets&&basiccut,mcordata,luminosity);
1396     TH1F *ZOSSFN = allsamples.Draw("ZOSSFN","-"+datajzb,1,cut,14000, xlabel, "events",cutmass&&cutOSSF&&cutnJets&&basiccut,mcordata,luminosity);
1397     TH1F *ZOSOFN = allsamples.Draw("ZOSOFN","-"+datajzb,1,cut,14000, xlabel, "events",cutmass&&cutOSOF&&cutnJets&&basiccut,mcordata,luminosity);
1398    
1399     TH1F *SBOSSFP = allsamples.Draw("SBOSSFP",datajzb,1,cut,14000, xlabel, "events",cutOSSF&&cutnJets&&basiccut&&sidebandcut,mcordata,luminosity);
1400     TH1F *SBOSOFP = allsamples.Draw("SBOSOFP",datajzb,1,cut,14000, xlabel, "events",cutOSOF&&cutnJets&&basiccut&&sidebandcut,mcordata,luminosity);
1401     TH1F *SBOSSFN = allsamples.Draw("SBOSSFN","-"+datajzb,1,cut,14000, xlabel, "events",cutOSSF&&cutnJets&&basiccut&&sidebandcut,mcordata,luminosity);
1402     TH1F *SBOSOFN = allsamples.Draw("SBOSOFN","-"+datajzb,1,cut,14000, xlabel, "events",cutOSOF&&cutnJets&&basiccut&&sidebandcut,mcordata,luminosity);
1403 buchmann 1.11
1404     cout << " Observed : " << ZOSSFP->Integral() << endl;
1405     cout << " Predicted: " << ZOSSFN->Integral() << " + (1/3)*(" << ZOSOFP->Integral() << "-" << ZOSOFN->Integral()<<") + (1/3)*(" << SBOSSFP->Integral() << "-" << SBOSSFN->Integral()<<") + (1/3)*(" << SBOSOFP->Integral() << "-" << SBOSOFN->Integral()<<")" << endl;
1406 buchmann 1.12 cout << " P(ZJets ) \t " << ZOSSFN->Integral() << endl;
1407     cout << " P(e&mu;]) \t " << ZOSOFP->Integral() << "-" << ZOSOFN->Integral() << " = " << ZOSOFP->Integral()-ZOSOFN->Integral()<<endl;
1408     cout << " P(ossf,sb]) \t " << SBOSSFP->Integral() << "-" << SBOSSFN->Integral()<<" = "<<SBOSSFP->Integral()-SBOSSFN->Integral()<<endl;
1409     cout << " P(osof,sb]) \t " << SBOSOFP->Integral() << "-" << SBOSOFN->Integral()<<" = "<<SBOSOFP->Integral()-SBOSOFN->Integral()<<endl;
1410 buchmann 1.11
1411     delete ZOSSFP;
1412     delete ZOSOFP;
1413     delete ZOSSFN;
1414     delete ZOSOFN;
1415    
1416     delete SBOSSFP;
1417     delete SBOSOFP;
1418     delete SBOSSFN;
1419     delete SBOSOFN;
1420    
1421     }
1422    
1423     void look_at_sidebands(string mcjzb, string datajzb) {
1424    
1425     cout << "Looking at sidebands ... " << endl;
1426 buchmann 1.12 int mcordata=data;//data // you can perform this study for mc or data ...
1427    
1428 buchmann 1.11
1429     int nbins=75; float min=51;float max=201;string xlabel="mll (GeV)";
1430     TCanvas *c1 = new TCanvas("c1","c1");
1431     c1->SetLogy(1);
1432     TH1F *datahistoOSSF = allsamples.Draw("datahistoOSSF","mll",nbins,min,max, xlabel, "events",cutOSSF&&cutnJets&&basiccut,data,luminosity);
1433     THStack mcstackOSSF = allsamples.DrawStack("mcstackOSSF","mll",nbins,min,max, xlabel, "events",cutOSSF&&cutnJets&&basiccut,mc,luminosity);
1434    
1435     datahistoOSSF->SetMinimum(1);
1436     datahistoOSSF->Draw();
1437     mcstackOSSF.Draw("same");
1438     datahistoOSSF->Draw("same");
1439     TLegend *kinleg = allsamples.allbglegend();
1440     kinleg->AddEntry(datahistoOSSF,"OSSF (data)","p");
1441     kinleg->Draw();
1442     CompleteSave(c1,"test/OSSF");
1443    
1444     TH1F *datahistoOSOF = allsamples.Draw("datahistoOSOF","mll",nbins,min,max, xlabel, "events",cutOSOF&&cutnJets&&basiccut,data,luminosity);
1445     THStack mcstackOSOF = allsamples.DrawStack("mcstackOSOF","mll",nbins,min,max, xlabel, "events",cutOSOF&&cutnJets&&basiccut,mc,luminosity);
1446     datahistoOSOF->SetMinimum(0.4);
1447     datahistoOSOF->Draw();
1448     mcstackOSOF.Draw("same");
1449     datahistoOSOF->Draw("same");
1450     TLegend *kinleg2 = allsamples.allbglegend();
1451     kinleg2->AddEntry(datahistoOSOF,"OSOF (data)","p");
1452     kinleg2->Draw();
1453     CompleteSave(c1,"test/OSOF");
1454    
1455 buchmann 1.12 TH1F *rawmlleemmData = allsamples.Draw("rawmlleemmData","mll",200,0,200, "JZB (GeV)", "events", cutOSSF&&cutnJets&&sidebandcut,data, luminosity);
1456 buchmann 1.11 TH1F *rawmllemData = allsamples.Draw("rawmllemData" ,"mll",200,0,200, "JZB (GeV)", "events", cutmass&&cutOSOF&&cutnJets,data, luminosity);
1457     cout << "Number of events in peak for OSOF: " << rawmllemData->GetEntries() << endl;
1458     cout << "Number of events in SB for OSSF: " << rawmlleemmData->GetEntries() << endl;
1459    
1460     TH1F *SFttbarZpeak = allsamples.Draw("SFttbarZpeak",mcjzb,100,-200,400, "JZB (GeV)", "events",cutmass&&cutOSSF&&cutnJets,mc,luminosity,allsamples.FindSample("TTJets"));
1461     TH1F *OFttbarZpeak = allsamples.Draw("OFttbarZpeak",mcjzb,100,-200,400, "JZB (GeV)", "events",cutmass&&cutOSOF&&cutnJets,mc,luminosity,allsamples.FindSample("TTJets"));
1462 buchmann 1.12 TH1F *SFttbarsideb = allsamples.Draw("SFttbarsideb",mcjzb,100,-200,400, "JZB (GeV)", "events",cutOSSF&&cutnJets&&sidebandcut,mc,luminosity,allsamples.FindSample("TTJets"));
1463 buchmann 1.11
1464     SFttbarZpeak->SetLineColor(kBlack);
1465     OFttbarZpeak->SetLineColor(kBlue);
1466     SFttbarsideb->SetLineColor(kPink);
1467    
1468     SFttbarZpeak->Draw("histo");
1469     OFttbarZpeak->Draw("histo,same");
1470     SFttbarsideb->Draw("histo,same");
1471    
1472     TLegend *leg3 = new TLegend(0.6,0.8,0.89,0.89);
1473     leg3->AddEntry(SFttbarZpeak,"SF ttbar Z peak","l");
1474     leg3->AddEntry(OFttbarZpeak,"OF ttbar Z peak","l");
1475     leg3->AddEntry(SFttbarsideb,"SF ttbar SB","l");
1476     leg3->SetFillColor(kWhite);
1477     leg3->SetLineColor(kWhite);
1478     leg3->Draw();
1479     CompleteSave(c1,"test/ttbar_comparison");
1480    
1481     cout << "Moving on to predicted / observed yields! " << endl;
1482 buchmann 1.12 cout << "Sideband definition: " << (const char*) sidebandcut << endl;
1483     do_experimental_pred_obs_calculation(50,mcjzb,datajzb,mcordata);
1484     do_experimental_pred_obs_calculation(75,mcjzb,datajzb,mcordata);
1485     do_experimental_pred_obs_calculation(100,mcjzb,datajzb,mcordata);
1486     do_experimental_pred_obs_calculation(125,mcjzb,datajzb,mcordata);
1487     do_experimental_pred_obs_calculation(150,mcjzb,datajzb,mcordata);
1488     }
1489    
1490     void pick_up_events(string cut) {
1491     cout << "Picking up events with cut " << cut << endl;
1492     allsamples.PickUpEvents(cut);
1493 buchmann 1.11 }
1494    
1495     void test() {
1496     TCanvas *testcanv = new TCanvas("testcanv","testcanv");
1497     testcanv->cd();
1498     switch_overunderflow(true);
1499     TH1F *ptdistr = allsamples.Draw("ptdistr","pt1",100,30,200, "p_{T} (GeV)", "events", cutOSSF,data,luminosity);
1500     switch_overunderflow(false);
1501     ptdistr->Draw();
1502     testcanv->SaveAs("test.png");
1503     cout << "HELLO there!" << endl;
1504     }