ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/Plotting_Functions.C
Revision: 1.24
Committed: Mon Jun 11 20:42:13 2012 UTC (12 years, 10 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.23: +16 -12 lines
Log Message:
Adapted peak finding for btag case; removed unnecessary legends on bpred plot (referring to fit that isn't drawn anymore); adapted bpred ratio range (otherwise lm4 is outside the drawn range);

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     #include <TProfile.h>
17     #include <TLegendEntry.h>
18    
19     using namespace std;
20    
21     using namespace PlottingSetup;
22    
23     void todo() {
24     dout << "My to do list: " << endl;
25     dout << " - ExperimentalModule::Poisson_ratio_plot : Get the second part to work!" << endl;
26 buchmann 1.5 dout << " - compare_onpeak_offpeak_signal_distributions: Implement this function ... " << endl;
27 buchmann 1.1 dout << "Info : The trigger requirement is currently set to " << (const char*) passtrig << endl;
28 buchmann 1.4 dout << "Info : The pt requirement is currently set to " << (const char*) passtrig << endl;
29     dout << "Info : The mll requirement is currently set to " << (const char*) cutmass << endl;
30     dout << "Info : The lepton requirement is currently set to " << (const char*) leptoncut << endl;
31 buchmann 1.1 }
32    
33     void find_peaks(float &MCPeak,float &MCPeakError, float &DataPeak, float &DataPeakError, float &MCSigma, float &DataSigma, stringstream &result, bool doPUreweighting = true )
34     {
35     // Temporarily switch off PU reweighting, if asked
36     TCut weightbackup=cutWeight;
37     if ( !doPUreweighting ) cutWeight="1.0";
38 buchmann 1.21
39 buchmann 1.24 int nbins=100;
40     if(PlottingSetup::DoBTag) nbins=25;
41    
42 buchmann 1.1 TCanvas *tempcan = new TCanvas("tempcan","Temporary canvas for peak finding preparations");
43 buchmann 1.24 TH1F *rawJZBeemmMC = allsamples.Draw("rawJZBeemmMC",jzbvariablemc,nbins,-50,50, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc, luminosity);
44     TH1F *rawJZBeemmData = allsamples.Draw("rawJZBeemmData",jzbvariabledata,nbins, -50,50, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data, luminosity);
45     TH1F *rawJZBemMC = allsamples.Draw("rawJZBemMC",jzbvariablemc,nbins,-50,50, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc, luminosity);
46     TH1F *rawJZBemData = allsamples.Draw("rawJZBemData",jzbvariabledata,nbins, -50,50, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data, luminosity);
47 buchmann 1.1 TH1F *rawttbarjzbeemmMC;
48    
49     if(method==doKM) {
50     //we only need this histo for the KM fitting...
51 buchmann 1.24 rawttbarjzbeemmMC = allsamples.Draw("rawttbarjzbeemmMC",jzbvariablemc,nbins, -50,50, "JZB [GeV]", "events",cutmass&&cutOSSF&&cutnJets,mc,luminosity,allsamples.FindSample("TTJet"));
52 buchmann 1.1 MCPeak=find_peak(rawJZBeemmMC, rawttbarjzbeemmMC, -40, 40, mc, MCPeakError,MCSigma,method);
53     DataPeak=find_peak(rawJZBeemmData, rawJZBeemmData, -40, 40, data, DataPeakError,DataSigma,method);
54     }
55     else {
56     TH1F *reducedMC = (TH1F*)rawJZBeemmMC->Clone();
57     TH1F *reducedData = (TH1F*)rawJZBeemmData->Clone();
58     reducedMC->Add(rawJZBemMC,-1);
59     reducedData->Add(rawJZBemData,-1);
60     //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)
61     MCPeak=find_peak(reducedMC, rawttbarjzbeemmMC, -40, 40, mc, MCPeakError,MCSigma,method);
62     DataPeak=find_peak(reducedData, rawJZBeemmData, -40, 40, data, DataPeakError,DataSigma,method);
63    
64     }
65    
66     // Revert to original PU reweighting
67     if ( !doPUreweighting ) cutWeight = weightbackup;
68    
69     // MCPeak=find_peak(rawJZBeemmMC, rawttbarjzbeemmMC, -40, 40, mc, MCPeakError,MCSigma,method);
70     // DataPeak=find_peak(rawJZBeemmData, rawJZBeemmData, -40, 40, data, DataPeakError,DataSigma,method);
71     dout << "We have found the peak in Data at " << DataPeak << " +/- " << DataPeakError << " with sigma=" << DataSigma << " +/- ?? (not impl.)" << endl;
72     result << "We have found the peak in Data at " << DataPeak << " +/- " << DataPeakError << " with sigma=" << DataSigma << " +/- ?? (not impl.)" << endl;
73     dout << "We have found the peak in MC at " << MCPeak << " +/- " << MCPeakError << " with sigma=" << MCSigma << " +/- ?? (not impl.)" << endl;
74     result << "We have found the peak in MC at " << MCPeak << " +/- " << MCPeakError << " with sigma=" << MCSigma << " +/- ?? (not impl.)" << endl;
75     delete tempcan;
76     }
77    
78 buchmann 1.11 void make_special_obs_pred_mll_plot(string mcjzb, float jzbthreshold) {
79     float min=70.0;
80     float max=115.0;
81     if(!PlottingSetup::RestrictToMassPeak) {
82     min=10;
83     max=150;
84     }
85     int nbins=int((max-min)/5);
86    
87     TCanvas *ckin = new TCanvas("ckin","Kinematic Plots (in the making)",600,600);
88    
89     stringstream largerzeroS;
90     largerzeroS << "(" << mcjzb << ">" << jzbthreshold << ")";
91     TCut largerzero(largerzeroS.str().c_str());
92    
93     stringstream smallerzeroS;
94     smallerzeroS << "(" << mcjzb << "<-" << jzbthreshold << ")";
95     TCut smallerzero(smallerzeroS.str().c_str());
96    
97     TH1F *RcorrJZBeemm = allsamples.Draw("RcorrJZBeemm", "mll",nbins,min,max, "m_{ll} [GeV]", "events", cutmass&&cutOSSF&&cutnJets&&largerzero,data,luminosity);
98     THStack mcRcorrJZBeemm = allsamples.DrawStack("mcRcorrJZBeemm","mll",nbins,min,max, "m_{ll} [GeV}", "events", cutmass&&cutOSSF&&cutnJets&&largerzero,mc,luminosity);
99     TH1F *LcorrJZBeemm = allsamples.Draw("LcorrJZBeemm", "mll",nbins,min,max, "m_{ll} [GeV]", "events", cutmass&&cutOSSF&&cutnJets&&smallerzero,data,luminosity);
100     TH1F *RcorrJZBem = allsamples.Draw("RcorrJZBem", "mll",nbins,min,max, "m_{ll} [GeV]", "events", cutmass&&cutOSOF&&cutnJets&&largerzero,data,luminosity);
101     TH1F *LcorrJZBem = allsamples.Draw("LcorrJZBem", "mll",nbins,min,max, "m_{ll} [GeV]", "events", cutmass&&cutOSOF&&cutnJets&&smallerzero,data,luminosity);
102    
103     TH1F *RcorrJZBSBem;
104     TH1F *LcorrJZBSBem;
105     TH1F *RcorrJZBSBeemm;
106     TH1F *LcorrJZBSBeemm;
107    
108 buchmann 1.16 // TH1F *RcorrJZBeemmNoS;
109 buchmann 1.11
110     if(PlottingSetup::RestrictToMassPeak) {
111     RcorrJZBSBem = allsamples.Draw("RcorrJZBSBem", "mll",nbins,min,max, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets&&largerzero,data, luminosity);
112     LcorrJZBSBem = allsamples.Draw("LcorrJZBSBem", "mll",nbins,min,max, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets&&smallerzero,data, luminosity);
113     RcorrJZBSBeemm = allsamples.Draw("RcorrJZBSBeemm","mll",nbins,min,max, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets&&largerzero,data, luminosity);
114     LcorrJZBSBeemm = allsamples.Draw("LcorrJZBSBeemm","mll",nbins,min,max, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets&&smallerzero,data, luminosity);
115     }
116    
117     TH1F *Bpred = (TH1F*)LcorrJZBeemm->Clone("Bpred");
118     if(PlottingSetup::RestrictToMassPeak) {
119     Bpred->Add(RcorrJZBem,1.0/3.);
120     Bpred->Add(LcorrJZBem,-1.0/3.);
121     Bpred->Add(RcorrJZBSBem,1.0/3.);
122     Bpred->Add(LcorrJZBSBem,-1.0/3.);
123     Bpred->Add(RcorrJZBSBeemm,1.0/3.);
124     Bpred->Add(LcorrJZBSBeemm,-1.0/3.);
125     } else {
126     Bpred->Add(RcorrJZBem,1.0);
127     Bpred->Add(LcorrJZBem,-1.0);
128     }
129    
130     Bpred->SetLineColor(kRed);
131    
132     RcorrJZBeemm->Draw();
133     mcRcorrJZBeemm.Draw("same");
134     Bpred->Draw("histo,same");
135     RcorrJZBeemm->Draw("same");
136    
137     TLegend *leg = allsamples.allbglegend();
138     leg->SetX1(0.58);
139     leg->AddEntry(RcorrJZBeemm,"observed (data)","l");
140     leg->AddEntry(Bpred,"predicted (data)","l");
141     leg->Draw("same");
142    
143     stringstream saveas;
144     saveas << "kin/Mll_After_Cut/Cut_At" << jzbthreshold;
145     CompleteSave(ckin,saveas.str());
146    
147     delete RcorrJZBeemm;
148     delete LcorrJZBeemm;
149     delete RcorrJZBem;
150     delete LcorrJZBem;
151     if(PlottingSetup::RestrictToMassPeak) {
152     delete RcorrJZBSBeemm;
153     delete LcorrJZBSBeemm;
154     delete RcorrJZBSBem;
155     delete LcorrJZBSBem;
156     }
157     delete Bpred;
158     delete ckin;
159     }
160    
161 buchmann 1.1 void make_special_mll_plot(int nbins, float min, float max, bool logscale,string xlabel) {
162    
163     TCanvas *ckin = new TCanvas("ckin","Kinematic Plots (in the making)",600,600);
164    
165     TH1F *datahistoOSSF = allsamples.Draw("datahistoOSSF","mll",nbins,min,max, xlabel, "events",cutOSSF&&cutnJets&&basiccut,data,luminosity);
166     THStack mcstackOSSF = allsamples.DrawStack("mcstackOSSF","mll",nbins,min,max, xlabel, "events",cutOSSF&&cutnJets&&basiccut,mc,luminosity);
167     TH1F *datahistoOSOF = allsamples.Draw("datahistoOSOF","mll",nbins,min,max, xlabel, "events",cutOSOF&&cutnJets&&basiccut,data,luminosity);
168    
169     if(logscale) ckin->SetLogy(1);
170     datahistoOSSF->SetMarkerSize(DataMarkerSize);
171     datahistoOSSF->GetXaxis()->SetTitle(xlabel.c_str());
172     datahistoOSSF->GetXaxis()->CenterTitle();
173     datahistoOSSF->GetYaxis()->SetTitle("events");
174     datahistoOSSF->GetYaxis()->CenterTitle();
175     datahistoOSOF->SetMarkerSize(DataMarkerSize);
176     datahistoOSSF->SetMarkerSize(DataMarkerSize);
177     datahistoOSSF->Draw();
178    
179     mcstackOSSF.Draw("same");
180     datahistoOSSF->Draw("same");
181    
182     datahistoOSOF->SetMarkerColor(TColor::GetColor("#FE642E"));
183     datahistoOSOF->SetLineColor(kRed);
184     datahistoOSOF->SetMarkerStyle(21);
185     datahistoOSOF->Draw("same");
186    
187     // Try to re-arrange legend...
188     TLegend *bgleg = allsamples.allbglegend("",datahistoOSSF);
189     TLegend *kinleg = make_legend();
190     kinleg->AddEntry(datahistoOSSF,"SF (data)","p");
191     kinleg->AddEntry(datahistoOSOF,"OF (data)","p");
192     TIter next(bgleg->GetListOfPrimitives());
193     TObject* obj;
194     // Copy the nice bkgd legend skipping the "data"
195 buchmann 1.16 while ( (obj = next()) )
196 buchmann 1.1 if ( strcmp(((TLegendEntry*)obj)->GetObject()->GetName(),"datahistoOSSF") )
197     kinleg->GetListOfPrimitives()->Add(obj);
198    
199     kinleg->Draw();
200     CompleteSave(ckin,"kin/mll_ossf_osof_distribution");
201    
202     delete datahistoOSOF;
203     delete datahistoOSSF;
204     delete ckin;
205     }
206    
207    
208     void draw_ratio_plot(TH1* hdata, THStack& hmc, float ymin=0.5, float ymax=1.5) {
209    
210     // Make a histogram from stack
211     TIter next(hmc.GetHists());
212     TObject* obj;
213     TH1* hratio = NULL;
214 buchmann 1.16 while ( (obj = next()) ) {
215 buchmann 1.1 if ( !hratio ) {
216     hratio = (TH1*)obj->Clone();
217     hratio->SetName("hratio");
218     } else hratio->Add( (TH1*)obj );
219     }
220     hratio->Divide(hdata);
221     hratio->SetMaximum(ymax);
222     hratio->SetMinimum(ymin);
223     hratio->SetMarkerStyle(2);
224     hratio->SetLineWidth(1);
225     hratio->GetYaxis()->SetLabelSize(0.08);
226     hratio->GetXaxis()->SetLabelSize(0.0);
227    
228     TPad* rpad = new TPad("rpad","",0.15,0.73,0.4,0.88);
229     rpad->SetTopMargin(0.0);
230     rpad->SetBottomMargin(0.0);
231     rpad->SetRightMargin(0.0);
232     rpad->Draw();
233     rpad->cd();
234     // hratio->GetXaxis()->SetNdivisions(0);
235     hratio->GetYaxis()->SetNdivisions(502,false);
236     hratio->Draw("e1x0");
237    
238     TF1* oneline = new TF1("","1.0",0,1000);
239     oneline->SetLineColor(kBlue);
240     oneline->SetLineStyle(1);
241     oneline->SetLineWidth(1);
242     oneline->Draw("same");
243     }
244    
245     float lastrange_min=0;
246     float lastrange_max=0;
247    
248     void make_kin_plot(string variable, string addcut, int nbins, float min, float max, bool logscale,
249     string xlabel, string filename, bool isPF=true, bool plotratio=true, bool loadlastminmax=false ) {
250     // TCut basiccut("(pfJetGoodNum>=2&&pfJetGoodID[0])&&(pfJetGoodNum>=2&&pfJetGoodID[1])&&(passed_triggers||!is_data)");
251     TCut ibasiccut=basiccut;
252     bool draw_separation_lines=false;
253    
254     if(isPF) ibasiccut=basiccut&&"pfjzb[0]>-998";
255    
256     if(addcut != "") ibasiccut = ibasiccut && addcut.c_str();
257    
258     //Step 1: Adapt the variable (if we're dealing with PF we need to adapt the variable!)
259     if(isPF) {
260     if(variable=="mll") variable="pfmll[0]";
261     if(variable=="jetpt[0]") variable="pfJetGoodPt[0]";
262     if(variable=="jeteta[0]") variable="pfJetGoodEta[0]";
263     if(variable=="pt") variable="pfpt[0]";
264     if(variable=="pt1") variable="pfpt1[0]";
265     if(variable=="pt2") variable="pfpt2[0]";
266     if(variable=="eta1") variable="pfeta1[0]";
267     if(variable=="jzb[1]") variable="pfjzb[0]";
268     //if(variable=="pfJetGoodNum") variable="pfJetGoodNum"; // pointless.
269     }
270    
271     //Step 2: Refine the cut
272     TCut cut;
273     cut=cutmass&&cutOSSF&&cutnJets&&ibasiccut;
274     if(filename=="nJets") cut=cutmass&&cutOSSF&&ibasiccut;
275 buchmann 1.22 if(filename=="nJets_osof") cut=cutmass&&cutOSOF&&ibasiccut;
276 buchmann 1.1 if(filename=="nJets_nocuts_except_mll_ossf") cut=cutmass&&cutOSSF;
277     if(filename=="mll") {
278     cut=cutOSSF&&cutnJets&&ibasiccut;
279     draw_separation_lines=true;
280     }
281     if(filename=="mll_ee") cut=cutOSSF&&cutnJets&&ibasiccut&&"id1==0";
282     if(filename=="mll_osof") {
283     cut=cutOSOF&&cutnJets&&ibasiccut;
284     draw_separation_lines=true;
285     }
286     if(filename=="mll_mm") cut=cutOSSF&&cutnJets&&ibasiccut&&"id1==1";
287 buchmann 1.11 if(Contains(filename,"mll_aboveJZB100")) cut=cutOSSF&&cutnJets&&ibasiccut;
288     if(Contains(filename,"mll_osof_aboveJZB")) cut=cutOSOF&&cutnJets&&ibasiccut;
289 buchmann 1.1 if(filename=="mll_inclusive"||filename=="mll_inclusive_highrange") cut=cutOSSF;
290     if(filename=="mll_inclusive_osof") cut=cutOSOF;
291     if(filename=="mll_inclusive_ee") cut=cutOSSF&&"id1==0";
292     if(filename=="mll_inclusive_mm") cut=cutOSSF&&"id1==1";
293     if(filename=="pfJetGoodEta_0") cut=cutOSSF&&cutmass&&ibasiccut&&cutnJets;
294     if(filename=="pfJetGoodPt_0") cut=cutOSSF&&cutmass&&ibasiccut&&cutnJets;
295    
296     TCanvas *ckin = new TCanvas("ckin","Kinematic Plots (in the making)",600,600);
297     ckin->SetLogy(logscale);
298     TH1F *datahisto = allsamples.Draw("datahisto",variable,nbins,min,max, xlabel, "events",cut,data,luminosity);
299     datahisto->SetMarkerSize(DataMarkerSize);
300     THStack mcstack = allsamples.DrawStack("mcstack",variable,nbins,min,max, xlabel, "events",cut,mc,luminosity);
301     if(variable=="pfJetGoodPt[0]") datahisto->SetMaximum(10*datahisto->GetMaximum());
302     if(variable=="pt") datahisto->SetMaximum(10*datahisto->GetMaximum());
303 buchmann 1.22 if(filename=="mll_inclusive"||filename=="mll_inclusive_mm"||filename=="mll_inclusive_ee") datahisto->SetMinimum(1);
304 buchmann 1.1 if(filename=="mll_osof") datahisto->SetMaximum(10*datahisto->GetMaximum());
305     if(filename=="mll_osof") datahisto->SetMinimum(9);
306     datahisto->SetMaximum(5.3*datahisto->GetMaximum());
307     datahisto->Draw("e1");
308     ckin->Update();
309     if(loadlastminmax) {
310     datahisto->SetMinimum(lastrange_min);
311     datahisto->SetMaximum(lastrange_max);
312     if(logscale) {
313     datahisto->SetMinimum(pow(10,lastrange_min));
314     datahisto->SetMaximum(pow(10,lastrange_max));
315     }
316     }
317    
318     mcstack.Draw("same");
319     datahisto->Draw("same,e1");
320     TLegend *kinleg = allsamples.allbglegend();
321     kinleg->Draw();
322     if(filename=="mll_osof") {
323     kinleg->SetHeader("Opposite flavor");
324     kinleg->SetX1(0.58);
325     }
326     if(filename=="mll") {
327     kinleg->SetHeader("Same flavor");
328     kinleg->SetX1(0.58);
329     }
330     TText* write_cut = write_cut_on_canvas(decipher_cut(cut,basicqualitycut));
331     write_cut->Draw();
332     TText* write_variable = write_text(0.99,0.01,variable);
333     write_variable->SetTextAlign(31);
334     write_variable->SetTextSize(0.02);
335    
336     TLine *lowerboundary;
337     TLine *upperboundary;
338    
339     if(RestrictToMassPeak&&draw_separation_lines) {
340     Color_t linecolor=kBlue;
341     float linemin=pow(10,ckin->GetUymin());
342     if(filename=="mll_osof") linemin=pow(10,lastrange_min);
343     lowerboundary = new TLine(71,linemin,71,datahisto->GetMaximum());
344     upperboundary = new TLine(111,linemin,111,datahisto->GetMaximum());
345     lowerboundary->SetLineColor(linecolor);
346     lowerboundary->SetLineStyle(2);
347     upperboundary->SetLineColor(linecolor);
348     upperboundary->SetLineStyle(2);
349     }
350    
351     lastrange_min=ckin->GetUymin();
352     lastrange_max=ckin->GetUymax();
353    
354    
355     if ( plotratio ) {
356     TPad *kinpad = new TPad("kinpad","kinpad",0,0,1,1);
357     kinpad->cd();
358     kinpad->SetLogy(logscale);
359     datahisto->Draw("e1");
360     mcstack.Draw("same");
361     datahisto->Draw("same,e1");
362     datahisto->Draw("same,axis");
363     if(RestrictToMassPeak&&draw_separation_lines) {
364     lowerboundary->Draw("same");
365     upperboundary->Draw("same");
366     }
367    
368     kinleg->Draw();
369     write_cut->Draw();
370     DrawPrelim();
371     string saveas="kin/"+filename;
372     if(isPF) saveas="kin/"+filename+"__PF";
373     save_with_ratio(datahisto,mcstack,kinpad->cd(),saveas);
374     // if(isPF) CompleteSave(with_ratio,"kin/"+filename+"__PF_withratio");
375     // else CompleteSave(with_ratio,"kin/"+filename+"_withratio");
376     // delete with_ratio;
377     } else {
378     if(isPF) CompleteSave(ckin,"kin/"+filename+"__PF");
379     else CompleteSave(ckin,"kin/"+filename);
380     }
381     datahisto->Delete();
382     delete ckin;
383     }
384    
385     void make_JES_plot() {
386    
387     int nbins=10;
388     float min=-0.5;
389     float max = 9.5;
390     bool logscale=true;
391     string xlabel="nJets";
392     TCut cut=cutmass&&cutOSSF&&basiccut;
393    
394     TCanvas *ckin = new TCanvas("ckin","Kinematic Plots (in the making)",600,600);
395     ckin->SetLogy(logscale);
396     TH1F *datahisto = allsamples.Draw("datahisto","pfJetGoodNum",nbins,min,max, xlabel, "events",cut,data,luminosity);
397     datahisto->SetMarkerSize(DataMarkerSize);
398     THStack mcstack = allsamples.DrawStack("mcstack","pfJetGoodNum",nbins,min,max, xlabel, "events",cut,mc,luminosity);
399     TH1F *JESup = allsamples.Draw("JESup","pfJetGoodNum315",nbins,min,max, xlabel, "events",cut,mc,luminosity);
400     TH1F *JESdn = allsamples.Draw("JESdn","pfJetGoodNum285",nbins,min,max, xlabel, "events",cut,mc,luminosity);
401    
402     datahisto->SetMinimum(1);
403     datahisto->SetMaximum(5.3*datahisto->GetMaximum()); // in line with kinematic plots style
404    
405     float xs[nbins],ys[nbins],exs[nbins],eys[nbins];
406     for(int i=1;i<JESup->GetNbinsX();i++) {
407     float up=JESup->GetBinContent(i);
408     float dn=JESdn->GetBinContent(i);
409     xs[i]=JESup->GetBinCenter(i);
410     ys[i]=0.5*(up+dn);
411     exs[i]=0.5*JESup->GetBinWidth(i);
412     eys[i]=0.5*TMath::Abs(up-dn);
413     }
414    
415     TGraphAsymmErrors *JESunc = new TGraphAsymmErrors(nbins, xs,ys,exs,exs,eys,eys);
416     JESunc->SetFillColor(TColor::GetColor("#00ADE1"));
417     JESunc->SetFillStyle(3002);
418     datahisto->Draw("e1");
419     mcstack.Draw("same");
420     JESunc->Draw("2");
421     datahisto->Draw("same,e1");
422     TLegend *kinleg = allsamples.allbglegend();
423     kinleg->AddEntry(JESunc,"JES uncertainty","f");
424     kinleg->Draw();
425     CompleteSave(ckin,"Systematics/JES");
426     datahisto->Delete();
427     delete ckin;
428    
429     }
430    
431     void do_kinematic_plots(string mcjzb, string datajzb, bool doPF=false)
432     {
433     bool dolog=true;
434     bool nolog=false;
435     if(doPF) write_warning(__FUNCTION__,"Please use caution when trying to produce PF plots; not all versions of the JZB trees have these variables!");
436 buchmann 1.17 float mll_low=50;
437 buchmann 1.2 float mll_hi=160;
438     if(!PlottingSetup::RestrictToMassPeak) {
439 buchmann 1.17 mll_low=20;
440 buchmann 1.11 mll_hi=210;
441 buchmann 1.2 }
442 buchmann 1.16
443 buchmann 1.19 make_kin_plot("met[4]","",40,0,200,dolog,"MET [GeV]","met",doPF,true);
444 buchmann 1.2 make_kin_plot("mll","",(int)((mll_hi-mll_low)/5),mll_low,mll_hi,dolog,"m_{ll} [GeV]","mll",doPF,true);
445     make_kin_plot("mll","",(int)((mll_hi-mll_low)/5),mll_low,mll_hi,dolog,"m_{ll} [GeV]","mll_osof",doPF,true,true);
446     make_kin_plot("mll","",(int)((mll_hi-mll_low)/5),mll_low,mll_hi,dolog,"m_{ll} [GeV]","mll_ee",doPF,true);
447     make_kin_plot("mll","",(int)((mll_hi-mll_low)/5),mll_low,mll_hi,dolog,"m_{ll} [GeV]","mll_mm",doPF,true);
448     make_kin_plot("mll","",(int)((mll_hi-mll_low)),mll_low,mll_hi,dolog,"m_{ll} [GeV]","mll_inclusive",doPF,true);
449     make_kin_plot("mll","",(int)((mll_hi-mll_low)),mll_low,mll_hi,dolog,"m_{ll} [GeV]","mll_inclusive_ee",doPF,true);
450     make_kin_plot("mll","",(int)((mll_hi-mll_low)),mll_low,mll_hi,dolog,"m_{ll} [GeV]","mll_inclusive_mm",doPF,true);
451 buchmann 1.22 make_kin_plot("mll","",(int)((mll_hi-mll_low))/2,mll_low,mll_hi,dolog,"m_{ll} [GeV]","mll_inclusive_osof",doPF,true);
452 buchmann 1.16 make_kin_plot("mll","",(int)((350-mll_low)),mll_low,350,dolog,"m_{ll} [GeV]","mll_inclusive_highrange",doPF);
453 buchmann 1.11 make_kin_plot("numVtx","",(int)(30.5-(-0.5)),-0.5,30.5,nolog,"N(Vtx)","numVtx",doPF);
454 buchmann 1.22 // make_kin_plot("jetpt[0]","",40,0,200,dolog,"leading jet p_{T} [GeV]","pfJetGoodPt_0",doPF);
455     // make_kin_plot("jeteta[0]","",40,-5,5,nolog,"leading jet #eta","pfJetGoodEta_0",doPF);
456 buchmann 1.1 make_kin_plot("pt","",50,0,400,dolog,"Z p_{T} [GeV]","Zpt",doPF);
457     make_kin_plot("pt1","",50,0,100,nolog,"p_{T} [GeV]","pt1",doPF);
458     make_kin_plot("pt2","",50,0,100,nolog,"p_{T} [GeV]","pt2",doPF);
459     make_kin_plot("eta1","",40,-5,5,nolog,"#eta_{l}","eta",doPF);
460     make_kin_plot("jzb[1]","",100,-150,150,dolog,"JZB [GeV]","jzb_ossf",doPF);
461 buchmann 1.22 make_kin_plot("pfJetGoodNum40","",9,-0.5,8.5,dolog,"nJets","nJets",doPF);
462     make_kin_plot("pfJetGoodNum40","",9,-0.5,8.5,dolog,"nJets","nJets_osof",doPF);
463     make_kin_plot("pfJetGoodNum40","",9,-0.5,8.5,dolog,"nJets","nJets_nocuts_except_mll_ossf",doPF);
464 buchmann 1.2 if(!doPF) make_special_mll_plot((int)((mll_hi-mll_low)/5),mll_low,mll_hi,dolog,"m_{ll} [GeV]");
465 buchmann 1.1 stringstream jzbcut;
466     jzbcut << "((is_data&&("<<datajzb<<")>100)||(!is_data&&("<<mcjzb<<")>100))";
467 buchmann 1.2 make_kin_plot("mll",jzbcut.str(),(int)((mll_hi-mll_low)/5),mll_low,mll_hi,dolog,"m_{ll} [GeV]","mll_aboveJZB100",doPF,true);
468 buchmann 1.11 make_kin_plot("mll",jzbcut.str(),(int)((mll_hi-mll_low)/5),mll_low,mll_hi,dolog,"m_{ll} [GeV]","mll_osof_aboveJZB100",doPF,true);
469     make_kin_plot("mll",jzbcut.str(),(int)((mll_hi-mll_low)/5),mll_low,mll_hi,dolog,"m_{ll} [GeV]","mll_aboveJZB100",doPF,true);
470     make_kin_plot("mll",jzbcut.str(),(int)((mll_hi-mll_low)/5),mll_low,mll_hi,dolog,"m_{ll} [GeV]","mll_osof_aboveJZB100",doPF,true);
471     stringstream jzbcut2;
472     jzbcut2 << "((is_data&&("<<datajzb<<")>150)||(!is_data&&("<<mcjzb<<")>150))";
473     make_kin_plot("mll",jzbcut2.str(),(int)((mll_hi-mll_low)/5),mll_low,mll_hi,dolog,"m_{ll} [GeV]","mll_aboveJZB150",doPF,true);
474     make_kin_plot("mll",jzbcut2.str(),(int)((mll_hi-mll_low)/5),mll_low,mll_hi,dolog,"m_{ll} [GeV]","mll_osof_aboveJZB150",doPF,true);
475     stringstream jzbcut3;
476     jzbcut3 << "((is_data&&("<<datajzb<<")>50)||(!is_data&&("<<mcjzb<<")>50))";
477 buchmann 1.24 make_kin_plot("mll",jzbcut3.str(),(int)((mll_hi-mll_low)/5),mll_low,mll_hi,dolog,"m_{ll} [GeV]","mll_aboveJZB50",doPF,true);
478     make_kin_plot("mll",jzbcut3.str(),(int)((mll_hi-mll_low)/5),mll_low,mll_hi,dolog,"m_{ll} [GeV]","mll_osof_aboveJZB50",doPF,true,true);
479 buchmann 1.16
480     make_special_obs_pred_mll_plot(mcjzb,0);
481 buchmann 1.11 make_special_obs_pred_mll_plot(mcjzb,50);
482     make_special_obs_pred_mll_plot(mcjzb,100);
483     make_special_obs_pred_mll_plot(mcjzb,150);
484     make_special_obs_pred_mll_plot(mcjzb,200);
485     make_special_obs_pred_mll_plot(mcjzb,250);
486 buchmann 1.1 }
487    
488     void make_comp_plot( string var, string xlabel, string filename, float jzbcut, string mcjzb, string datajzb,
489     int nbins, float xmin, float xmax, bool log,
490     float ymin=0, float ymax=0, bool leftJustified=false ) {
491     flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak ---- WATCH OUT: the argument in the function changed!
492    
493     TCut weightbackup=cutWeight;//backing up the correct weight (restoring below!)
494     if(weightbackup==TCut("1.0")||weightbackup==TCut("1")) write_warning(__FUNCTION__,"WATCH OUT THE WEIGHT HAS POSSIBLY NOT BEEN RESET!!!! PLEASE CHANGE LINE "+any2string(__LINE__));
495     if(var=="numVtx") cutWeight=TCut("1.0");
496     TCut jzbData[]= { TCut(TString(datajzb+">"+any2string(jzbcut))),TCut(TString(datajzb+"<-"+any2string(jzbcut))) };
497     TCut jzbMC[] = { TCut(TString(mcjzb+">"+any2string(jzbcut))),TCut(TString(mcjzb+"<-"+any2string(jzbcut))) };
498    
499     flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak ---- below: the next ~20 lines changed!
500     int nRegions=4;
501     if(!PlottingSetup::RestrictToMassPeak) {
502     nRegions=2;
503     }
504    
505     string sRegions[] = { "SFZP","OFZP","SFSB","OFSB" };
506     TCut kRegions[] = { cutOSSF&&cutnJets&&cutmass, cutOSOF&&cutnJets&&cutmass,
507     cutOSSF&&cutnJets&&sidebandcut, cutOSOF&&cutnJets&&sidebandcut };
508    
509     //find ymax
510     TH1F *Refdatahisto = allsamples.Draw("datahisto", var,nbins,xmin,xmax,xlabel,"events",kRegions[0]&&jzbData[0],data,luminosity);
511     ymax=int((Refdatahisto->GetMaximum()+2*TMath::Sqrt(Refdatahisto->GetMaximum()))*1.2+0.5);
512     delete Refdatahisto;
513    
514     for ( int iregion=0; iregion<nRegions; ++iregion )
515     for ( int ijzb=0; ijzb<2; ++ijzb ) {
516     TCanvas *ccomp = new TCanvas("ccomp","Comparison plot",600,400);
517     ccomp->SetLogy(log);
518     TH1F *datahisto = allsamples.Draw("datahisto", var,nbins,xmin,xmax,xlabel,"events",kRegions[iregion]&&jzbData[ijzb],data,luminosity);
519     TH1F *lm4histo = signalsamples.Draw("lm4histo", var,nbins,xmin,xmax,xlabel,"events",kRegions[iregion]&&jzbData[ijzb],data,luminosity,signalsamples.FindSample("LM4"));
520     THStack mcstack = allsamples.DrawStack("mcstack",var,nbins,xmin,xmax,xlabel,"events",kRegions[iregion]&&jzbMC[ijzb], mc, luminosity);
521     datahisto->SetMarkerSize(DataMarkerSize);
522     if (ymax>ymin) datahisto->SetMaximum(ymax);
523     lm4histo->SetLineStyle(2);
524     datahisto->Draw("e1");
525     mcstack.Draw("same");
526     datahisto->Draw("same,e1");
527     // lm4histo->Draw("hist,same");
528     TLegend *kinleg = allsamples.allbglegend((sRegions[iregion]+(ijzb?"neg":"pos")).c_str());
529     if ( leftJustified ) {
530     Float_t w = kinleg->GetX2()-kinleg->GetX1();
531     kinleg->SetX1(0.2);
532     kinleg->SetX2(0.2+w);
533     }
534     // kinleg->AddEntry(lm4histo,"LM4","l");
535     kinleg->Draw();
536     TText* write_variable = write_text(0.99,0.01,var);
537     write_variable->SetTextAlign(31);
538     write_variable->SetTextSize(0.02);
539     ccomp->RedrawAxis();
540     CompleteSave(ccomp,"compare/JZBcut_at_"+any2string(jzbcut)+"/"+filename+"/"+filename+sRegions[iregion]+(ijzb?"neg":"pos"));
541     delete datahisto;
542     delete ccomp;
543     delete lm4histo;
544     }
545     cutWeight=weightbackup;
546     }
547    
548    
549     void region_comparison_plots(string mcjzb, string datajzb, vector<float> jzb_cuts) {
550     dout << "Creating comparison plots for signal and control regions" << endl;
551     // Compare a few quantities in the signal region and all 7 control regions
552    
553     switch_overunderflow(true); // switching overflow/underflow bins on
554    
555    
556     flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak ---- the arguments changed
557 buchmann 1.16 for(int ijzb=0;ijzb<(int)jzb_cuts.size();ijzb++) {
558 buchmann 1.1 float jzbcut=jzb_cuts[ijzb]; // Comparison plots are done for this JZB cut
559     float mll_low=50;float mll_high=170;
560     if(!PlottingSetup::RestrictToMassPeak) {
561     mll_high=200;
562     mll_low=50;
563     }
564     make_comp_plot("mll","m_{ll} [GeV]","mll",jzbcut,mcjzb,datajzb,30,mll_low,mll_high,false,0,16.);
565     make_comp_plot("met[4]","pfMET [GeV]","pfmet",jzbcut,mcjzb,datajzb,18,0,360,false,0,16.);
566     make_comp_plot("pfJetGoodNum","#(jets)","njets",jzbcut,mcjzb,datajzb,10,0,10, false,0,35.);
567     make_comp_plot("pt","Z p_{T} [GeV]","Zpt",jzbcut,mcjzb,datajzb,26,0,525,false,0.,21.);
568     make_comp_plot("numVtx","#(prim. vertices)","nvtx",jzbcut,mcjzb,datajzb,20,0.,20.,false,0,16.);
569     make_comp_plot("TMath::Abs(dphi)","#Delta#phi(leptons)","dphilep",jzbcut,mcjzb,datajzb,10,0.,3.1415,false,0,16.,true);
570     make_comp_plot("TMath::Abs(dphi_sumJetVSZ[1])","#Delta#phi(Z,jets)","dphiZjets",jzbcut,mcjzb,datajzb,10,0.,3.1415,false,0,16.,true);
571     }
572    
573     switch_overunderflow(false); // switching overflow/underflow bins off
574     }
575    
576    
577    
578     void do_kinematic_PF_plots(string mcjzb, string datajzb)
579     {
580     do_kinematic_plots(mcjzb,datajzb,true);
581     }
582    
583     void signal_bg_comparison()
584     {
585     TCanvas *can = new TCanvas("can","Signal Background Comparison Canvas");
586     can->SetLogy(1);
587    
588     int sbg_nbins=130;
589     float sbg_min=-500; //-110;
590     float sbg_max=800; //jzbHigh;
591    
592     float simulatedlumi=luminosity;//in pb please - adjust to your likings
593    
594 buchmann 1.3 TH1F *JZBplotZJETs = allsamples.Draw("JZBplotZJETs",jzbvariablemc,sbg_nbins,sbg_min,sbg_max, "JZB [GeV]", "events",cutmass&&cutOSSF&&cutnJets,mc,simulatedlumi,allsamples.FindSample("/DY"));
595 buchmann 1.7 TH1F *JZBplotLM4;
596     if(PlottingSetup::RestrictToMassPeak) JZBplotLM4 = allsamples.Draw("JZBplotLM4",jzbvariablemc,sbg_nbins,sbg_min,sbg_max, "JZB [GeV]", "events",cutmass&&cutOSSF&&cutnJets,mc,simulatedlumi,allsamples.FindSample("LM4"));
597     else JZBplotLM4 = allsamples.Draw("JZBplotLM4",jzbvariablemc,sbg_nbins,sbg_min,sbg_max, "JZB [GeV]", "events",cutmass&&cutOSSF&&cutnJets,mc,simulatedlumi,allsamples.FindSample("LM3"));
598 buchmann 1.1 TH1F *JZBplotTtbar = allsamples.Draw("JZBplotTtbar",jzbvariablemc,sbg_nbins,sbg_min,sbg_max, "JZB [GeV]", "events",cutmass&&cutOSSF&&cutnJets,mc,simulatedlumi,allsamples.FindSample("TTJets"));
599    
600     JZBplotTtbar->SetLineColor(allsamples.GetColor("TTJet"));
601     JZBplotZJETs->SetFillColor(allsamples.GetColor("DY"));
602     JZBplotZJETs->SetLineColor(kBlack);
603     JZBplotLM4->SetLineStyle(2);
604     JZBplotZJETs->SetMaximum(JZBplotZJETs->GetMaximum()*5);
605     JZBplotZJETs->SetMinimum(1);
606    
607     JZBplotTtbar->SetMaximum(JZBplotZJETs->GetMaximum());
608     JZBplotTtbar->SetMinimum(0.01);
609     JZBplotTtbar->SetFillColor(allsamples.GetColor("TTJets"));
610     JZBplotTtbar->DrawClone("histo");
611     JZBplotZJETs->Draw("histo,same");
612     JZBplotTtbar->SetFillColor(0);
613     JZBplotTtbar->DrawClone("histo,same");
614     JZBplotTtbar->SetFillColor(allsamples.GetColor("TTJets"));
615     JZBplotLM4->Draw("histo,same");
616    
617    
618     TLegend *signal_bg_comparison_leg2 = make_legend("",0.55,0.75,false);
619     signal_bg_comparison_leg2->AddEntry(JZBplotZJETs,"Z+Jets","f");
620     signal_bg_comparison_leg2->AddEntry(JZBplotTtbar,"t#bar{t}","f");
621 buchmann 1.7 if(PlottingSetup::RestrictToMassPeak) signal_bg_comparison_leg2->AddEntry(JZBplotLM4,"LM4","f");
622     else signal_bg_comparison_leg2->AddEntry(JZBplotLM4,"LM3","f");
623 buchmann 1.1 signal_bg_comparison_leg2->Draw();
624     DrawMCPrelim(simulatedlumi);
625     CompleteSave(can,"jzb_bg_vs_signal_distribution");
626    
627     // Define illustrative set of SMS points
628     TCut kSMS1("MassGlu==250&&MassLSP==75");
629     TCut kSMS2("MassGlu==800&&MassLSP==200");
630     TCut kSMS3("MassGlu==1050&&MassLSP==850");
631     TCut kSMS4("MassGlu==1200&&MassLSP==100");
632    
633     //If the scan samples haven't been loaded yet, this is a good point to load them (all of them!)
634     if((scansample.collection).size()<2) define_SMS_sample(false, allsamples, signalsamples, scansample, true); // loading ALL zones for the scans, not only the basic one.
635    
636    
637     TH1F *JZBplotSMS1 = scansample.Draw("JZBplotSMS1",jzbvariablemc,sbg_nbins,sbg_min,sbg_max, "JZB [GeV]", "events",cutmass&&cutOSSF&&cutnJets&&kSMS1,mc,simulatedlumi,scansample.FindSample("t"));
638     JZBplotSMS1->Scale(JZBplotLM4->Integral()/JZBplotSMS1->Integral());
639    
640     TH1F *JZBplotSMS2 = scansample.Draw("JZBplotSMS2",jzbvariablemc,sbg_nbins,sbg_min,sbg_max, "JZB [GeV]", "events",cutmass&&cutOSSF&&cutnJets&&kSMS2,mc,simulatedlumi,scansample.FindSample("t"));
641     JZBplotSMS2->Scale(JZBplotLM4->Integral()/JZBplotSMS2->Integral());
642    
643     TH1F *JZBplotSMS3 = scansample.Draw("JZBplotSMS3",jzbvariablemc,sbg_nbins,sbg_min,sbg_max, "JZB [GeV]", "events",cutmass&&cutOSSF&&cutnJets&&kSMS3,mc,simulatedlumi,scansample.FindSample("t"));
644     JZBplotSMS3->Scale(JZBplotLM4->Integral()/JZBplotSMS3->Integral());
645    
646     TH1F *JZBplotSMS4 = scansample.Draw("JZBplotSMS4",jzbvariablemc,sbg_nbins,sbg_min,sbg_max, "JZB [GeV]", "events",cutmass&&cutOSSF&&cutnJets&&kSMS4,mc,simulatedlumi,scansample.FindSample("t"));
647     JZBplotSMS4->Scale(JZBplotLM4->Integral()/JZBplotSMS4->Integral());
648    
649     // Draw all plots overlaid
650     JZBplotTtbar->SetMaximum(JZBplotZJETs->GetMaximum());
651     JZBplotTtbar->SetMinimum(0.01);
652     JZBplotTtbar->SetFillColor(allsamples.GetColor("TTJets"));
653     JZBplotTtbar->DrawClone("histo");
654     JZBplotZJETs->Draw("histo,same");
655     JZBplotTtbar->SetFillColor(0);
656     JZBplotTtbar->DrawClone("histo,same");
657     JZBplotTtbar->SetFillColor(allsamples.GetColor("TTJets"));
658    
659     JZBplotSMS1->SetLineColor(kRed+1);
660     JZBplotSMS2->SetLineColor(kBlue+1);
661     JZBplotSMS3->SetLineColor(kRed+1);
662     JZBplotSMS4->SetLineColor(kBlue+1);
663     JZBplotSMS3->SetLineStyle(2);
664     JZBplotSMS4->SetLineStyle(2);
665    
666     JZBplotSMS1->Draw("histo,same");
667     JZBplotSMS2->Draw("histo,same");
668     JZBplotSMS3->Draw("histo,same");
669     JZBplotSMS4->Draw("histo,same");
670     JZBplotLM4->SetLineColor(kGreen);JZBplotLM4->Draw("histo,same");
671     TLegend *signal_bg_comparison_leg6 = make_legend("",0.55,0.55,false);
672     signal_bg_comparison_leg6->AddEntry(JZBplotZJETs,"Z+Jets","f");
673     signal_bg_comparison_leg6->AddEntry(JZBplotTtbar,"t#bar{t}","f");
674     signal_bg_comparison_leg6->AddEntry(JZBplotSMS1,"","");
675     signal_bg_comparison_leg6->AddEntry(JZBplotSMS1,"SMS parameters","");
676     signal_bg_comparison_leg6->AddEntry(JZBplotSMS1,"(250,75) [GeV]","f");
677     signal_bg_comparison_leg6->AddEntry(JZBplotSMS2,"(800,200) [GeV]","f");
678     signal_bg_comparison_leg6->AddEntry(JZBplotSMS3,"(1050,850) [GeV]","f");
679     signal_bg_comparison_leg6->AddEntry(JZBplotSMS4,"(1200,100) [GeV]","f");
680     signal_bg_comparison_leg6->AddEntry(JZBplotLM4,"LM4","f");
681     signal_bg_comparison_leg6->Draw();
682     DrawMCPrelim(simulatedlumi);
683     CompleteSave(can,"jzb_bg_vs_signal_distribution_SMS__summary");
684    
685     while((scansample.collection).size() > 1) scansample.RemoveLastSample();
686    
687     }
688    
689     vector<TF1*> do_cb_fit_to_plot(TH1F *histo, float Sigma, float doingfitacrosstheboard=false) {
690     TF1 *BpredFunc = new TF1("BpredFunc",InvCrystalBall,0,1000,5);
691     BpredFunc->SetParameter(0,histo->GetBinContent(1));
692     if(doingfitacrosstheboard) BpredFunc->SetParameter(0,histo->GetMaximum());
693     BpredFunc->SetParameter(1,0.);
694     if(method==1) BpredFunc->SetParameter(2,10*Sigma);//KM
695     else BpredFunc->SetParameter(2,Sigma);//Gaussian based methods
696     if(method==-99) BpredFunc->SetParameter(2,2.0*Sigma);//Kostas
697     BpredFunc->SetParameter(3,1.8);
698     BpredFunc->SetParameter(4,2.5);
699     histo->Fit(BpredFunc,"QN0");
700     BpredFunc->SetLineColor(kBlue);
701    
702     TF1 *BpredFuncP = new TF1("BpredFuncP",InvCrystalBallP,-1000,histo->GetBinLowEdge(histo->GetNbinsX())+histo->GetBinWidth(histo->GetNbinsX()),5);
703     TF1 *BpredFuncN = new TF1("BpredFuncN",InvCrystalBallN,-1000,histo->GetBinLowEdge(histo->GetNbinsX())+histo->GetBinWidth(histo->GetNbinsX()),5);
704    
705     BpredFuncP->SetParameters(BpredFunc->GetParameters());
706     BpredFuncP->SetLineColor(kBlue);
707     BpredFuncP->SetLineStyle(2);
708    
709     BpredFuncN->SetParameters(BpredFunc->GetParameters());
710     BpredFuncN->SetLineColor(kBlue);
711     BpredFuncN->SetLineStyle(2);
712    
713     vector<TF1*> functions;
714     functions.push_back(BpredFuncN);
715     functions.push_back(BpredFunc);
716     functions.push_back(BpredFuncP);
717     return functions;
718     }
719    
720    
721     TF1* do_logpar_fit_to_plot(TH1F *osof) {
722     TCanvas *logpar_fit_can = new TCanvas("logpar_fit_can","Fit canvas for LogPar");
723     TF1 *logparfunc = new TF1("logparfunc",LogParabola,0,300,3);
724     TF1 *logparfunc2 = new TF1("logparfunc2",LogParabola,0,(osof->GetXaxis()->GetBinLowEdge(osof->GetNbinsX())+osof->GetXaxis()->GetBinWidth(osof->GetNbinsX())),3);
725     TF1 *logparfuncN = new TF1("logparfuncN",LogParabolaN,0,(osof->GetXaxis()->GetBinLowEdge(osof->GetNbinsX())+osof->GetXaxis()->GetBinWidth(osof->GetNbinsX())),3);
726     TF1 *logparfuncP = new TF1("logparfuncP",LogParabolaP,0,(osof->GetXaxis()->GetBinLowEdge(osof->GetNbinsX())+osof->GetXaxis()->GetBinWidth(osof->GetNbinsX())),3);
727     osof->SetMinimum(0);
728     osof->Fit(logparfunc,"QR");
729     osof->Draw();
730     logparfunc->SetLineWidth(2);
731     logparfunc2->SetParameters(logparfunc->GetParameters());
732     logparfuncN->SetParameters(logparfunc->GetParameters());
733     logparfuncP->SetParameters(logparfunc->GetParameters());
734     stringstream fitinfo;
735     fitinfo << "#Chi^{2} / NDF : " << logparfunc->GetChisquare() << " / " << logparfunc->GetNDF();
736     TText *writefitinfo = write_text(0.8,0.8,fitinfo.str());
737     writefitinfo->SetTextSize(0.03);
738     DrawPrelim();
739     writefitinfo->Draw();
740     logparfunc->Draw("same");
741     logparfunc2->Draw("same");
742     logparfuncN->SetLineStyle(2);
743     logparfuncP->SetLineStyle(2);
744     logparfuncN->Draw("same");
745     logparfuncP->Draw("same");
746     CompleteSave(logpar_fit_can,"MakingOfBpredFunction/Bpred_Data_LogPar_Fit_To_TTbarPred");
747     delete logpar_fit_can;
748     return logparfunc2;
749     }
750    
751     vector<TF1*> do_extended_fit_to_plot(TH1F *prediction, TH1F *Tprediction, TH1F *ossf, TH1F *osof,int isdata) {
752     /* there are mainly two background contributions: Z+Jets (a) and ttbar (b). So:
753     a) The case is clear - we take the OSSF prediction - OSOF prediction, and fit a crystal ball function to it. We then extract the CB parameters.
754     b) For ttbar, we use the OSOF distribution and look at the [10,100] GeV JZB range, and fit our log parabola. We then extract the LP parameters.
755     Once we have these two components, we use the combined parameters to get the final function and we're done.
756     */
757     //Step 1: take the OSSF prediction - OSOF prediction, and fit a crystal ball function to it
758     TH1F *step1cb = (TH1F*)ossf->Clone("step1cb");
759     step1cb->Add(osof,-1);
760     vector<TF1*> functions = do_cb_fit_to_plot(step1cb,PlottingSetup::JZBPeakWidthData);
761     TF1 *zjetscrystalball = functions[1];
762    
763     //Step 2: use the OSOF distribution and look at the [10,100] GeV JZB range, and fit our log parabola
764     // TH1F *ttbarprediction=(TH1F*)prediction->Clone("ttbarprediction");
765     // ttbarprediction->Add(ossf,-1);//without the Z+Jets estimate, this is really just the ttbar estimate!
766     // the line above is not necessary anymore as we're now looking at a prediction without Z+Jets, and not multiplied with (1.0/3)
767     TF1 *ttbarlogpar = do_logpar_fit_to_plot(Tprediction);
768     flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
769     if(PlottingSetup::RestrictToMassPeak) ttbarlogpar->SetParameter(0,1.0/3*ttbarlogpar->GetParameter(0));//correcting for the fact that we didn't multiply with (1.0/3);
770    
771    
772     TF1 *ttbarlogparP = new TF1("ttbarlogparP",LogParabolaP,0,(prediction->GetXaxis()->GetBinLowEdge(prediction->GetNbinsX())+prediction->GetXaxis()->GetBinWidth(prediction->GetNbinsX())),8);
773     TF1 *ttbarlogparN = new TF1("ttbarlogparN",LogParabolaN,0,(prediction->GetXaxis()->GetBinLowEdge(prediction->GetNbinsX())+prediction->GetXaxis()->GetBinWidth(prediction->GetNbinsX())),8);
774    
775     //and now fuse the two!
776     TF1 *kmlp = new TF1("kmlp", CrystalBallPlusLogParabola, 0,(prediction->GetXaxis()->GetBinLowEdge(prediction->GetNbinsX())+prediction->GetXaxis()->GetBinWidth(prediction->GetNbinsX())),8);
777     TF1 *kmlpP= new TF1("kmlpP",CrystalBallPlusLogParabolaP,0,(prediction->GetXaxis()->GetBinLowEdge(prediction->GetNbinsX())+prediction->GetXaxis()->GetBinWidth(prediction->GetNbinsX())),8);
778     TF1 *kmlpN= new TF1("kmlpN",CrystalBallPlusLogParabolaN,0,(prediction->GetXaxis()->GetBinLowEdge(prediction->GetNbinsX())+prediction->GetXaxis()->GetBinWidth(prediction->GetNbinsX())),8);
779     double kmlp_pars[10];
780     for(int i=0;i<5;i++) kmlp_pars[i]=zjetscrystalball->GetParameter(i);
781     for(int i=0;i<3;i++) kmlp_pars[5+i]=ttbarlogpar->GetParameter(i);
782     ttbarlogparP->SetParameters(ttbarlogpar->GetParameters());
783     ttbarlogparN->SetParameters(ttbarlogpar->GetParameters());
784     kmlp->SetParameters(kmlp_pars);
785     prediction->Fit(kmlp,"Q");//fitting the final result (done this in the past but kicked it)
786     /*
787     if you want to start from scratch (without the partial fitting and only fitting the whole thing, some good start values could be :
788     */
789     kmlp_pars[0]=kmlp->GetParameter(0);
790     kmlp_pars[1]=3.6198;
791     kmlp_pars[2]=16.4664;
792     kmlp_pars[3]=1.92253;
793     kmlp_pars[4]=3.56099;
794     kmlp_pars[5]=5.83;
795     kmlp_pars[6]=0.000757479;
796     kmlp_pars[7]=95.6157;
797     kmlp_pars[8]=0;
798     kmlp_pars[9]=0;
799     kmlp->SetParameters(kmlp_pars);
800     /**/
801     prediction->Fit(kmlp,"Q");//fitting the final result (done this in the past but kicked it)
802    
803     kmlpP->SetParameters(kmlp->GetParameters());
804     kmlpN->SetParameters(kmlp->GetParameters());
805    
806     // now that we're done, let's save all of this so we can have a look at it afterwards.
807     TCanvas *can = new TCanvas("can","Prediction Fit Canvas");
808     can->SetLogy(1);
809     prediction->SetMarkerColor(kRed);
810     prediction->Draw();
811    
812     kmlp->SetLineColor(TColor::GetColor("#04B404"));
813     kmlpP->SetLineColor(TColor::GetColor("#04B404"));
814     kmlpN->SetLineColor(TColor::GetColor("#04B404"));
815     kmlp->Draw("same");
816     kmlpN->SetLineStyle(2);
817     kmlpP->SetLineStyle(2);
818     kmlpN->Draw("same");
819     kmlpP->Draw("same");
820    
821     ttbarlogpar->SetLineColor(TColor::GetColor("#CC2EFA"));//purple
822     ttbarlogpar->Draw("same");
823     ttbarlogparP->SetLineColor(TColor::GetColor("#CC2EFA"));//purple
824     ttbarlogparN->SetLineColor(TColor::GetColor("#CC2EFA"));//purple
825     ttbarlogparP->SetLineStyle(2);
826     ttbarlogparN->SetLineStyle(2);
827     ttbarlogparP->Draw("same");
828     ttbarlogparN->Draw("same");
829    
830     functions[0]->Draw("same");functions[1]->Draw("same");functions[2]->Draw("same");
831    
832     TLegend *analyticalBpredLEG = make_legend("",0.5,0.55);
833     analyticalBpredLEG->AddEntry(prediction,"predicted","p");
834     analyticalBpredLEG->AddEntry(functions[1],"Crystal Ball fit","l");
835     analyticalBpredLEG->AddEntry(functions[0],"1#sigma Crystal Ball fit","l");
836     analyticalBpredLEG->AddEntry(ttbarlogparN,"TTbar fit","l");
837     analyticalBpredLEG->AddEntry(ttbarlogpar,"1#sigma TTbar fit","l");
838     analyticalBpredLEG->AddEntry(kmlp,"Combined function","l");
839     analyticalBpredLEG->AddEntry(kmlpN,"1#sigma combined function","l");
840     analyticalBpredLEG->Draw("same");
841    
842     if(isdata==0) CompleteSave(can,"MakingOfBpredFunction/Bpred_MC_Analytical_Function_Composition");
843     if(isdata==1) CompleteSave(can,"MakingOfBpredFunction/Bpred_data_Analytical_Function_Composition");
844     if(isdata==2) CompleteSave(can,"MakingOfBpredFunction/Bpred_MCBnS_Analytical_Function_Composition");
845     delete can;
846    
847     //and finally: prep return functions
848     vector<TF1*> return_functions;
849     return_functions.push_back(kmlpN);
850     return_functions.push_back(kmlp);
851     return_functions.push_back(kmlpP);
852    
853     return_functions.push_back(ttbarlogparN);
854     return_functions.push_back(ttbarlogpar);
855     return_functions.push_back(ttbarlogparP);
856    
857     return_functions.push_back(functions[0]);
858     return_functions.push_back(functions[1]);
859     return_functions.push_back(functions[2]);
860    
861     return return_functions;
862     }
863    
864 buchmann 1.12 void do_prediction_plot(string jzb, TCanvas *globalcanvas, float sigma, float high, int use_data, bool overlay_signal = false,string subdir="" )
865 buchmann 1.1 {
866 buchmann 1.12 // note: sigma is not used ATM
867 buchmann 1.10 switch_overunderflow(true);
868 buchmann 1.1 bool is_data=false;
869     bool use_signal=false;
870     if(use_data==1) is_data=true;
871     if(use_data==2) use_signal=true;
872     int nbins=50;//100;
873     if(is_data) nbins=50;
874     float low=0;
875     float hi=500;
876 buchmann 1.10
877     stringstream largerzeroS;
878     largerzeroS << "("<<jzb<<">0)";
879     TCut largerzero(largerzeroS.str().c_str());
880     stringstream smallerzeroS;
881     smallerzeroS << "("<<jzb<<"<0)";
882     TCut smallerzero(smallerzeroS.str().c_str());
883    
884 buchmann 1.1 TH1F *blankback = new TH1F("blankback","blankback",int(high/10),0,high);
885 buchmann 1.10 TH1F *RcorrJZBeemm = allsamples.Draw("RcorrJZBeemm",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets&&largerzero,is_data, luminosity,use_signal);
886     TH1F *LcorrJZBeemm = allsamples.Draw("LcorrJZBeemm",("-"+jzb).c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets&&smallerzero,is_data, luminosity,use_signal);
887     TH1F *RcorrJZBem = allsamples.Draw("RcorrJZBem",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets&&largerzero,is_data, luminosity,use_signal);
888     TH1F *LcorrJZBem = allsamples.Draw("LcorrJZBem",("-"+jzb).c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets&&smallerzero,is_data, luminosity,use_signal);
889    
890 buchmann 1.1 blankback->GetXaxis()->SetTitle(RcorrJZBeemm->GetXaxis()->GetTitle());
891     blankback->GetYaxis()->SetTitle(RcorrJZBeemm->GetYaxis()->GetTitle());
892     blankback->GetXaxis()->CenterTitle();
893     blankback->GetYaxis()->CenterTitle();
894    
895     flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak ---- prediction changed.
896     TH1F *RcorrJZBSBem;
897     TH1F *LcorrJZBSBem;
898     TH1F *RcorrJZBSBeemm;
899     TH1F *LcorrJZBSBeemm;
900    
901     TH1F *RcorrJZBeemmNoS;
902    
903     //these are for the ratio
904 buchmann 1.10
905     TH1F *JRcorrJZBeemm = allsamples.Draw("JRcorrJZBeemm",jzb.c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets&&largerzero,is_data, luminosity,use_signal);
906     TH1F *JLcorrJZBeemm = allsamples.Draw("JLcorrJZBeemm",("-"+jzb).c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets&&smallerzero,is_data, luminosity,use_signal);
907     TH1F *JRcorrJZBem = allsamples.Draw("JRcorrJZBem",jzb.c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets&&largerzero,is_data, luminosity,use_signal);
908     TH1F *JLcorrJZBem = allsamples.Draw("JLcorrJZBem",("-"+jzb).c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets&&smallerzero,is_data, luminosity,use_signal);
909 buchmann 1.1
910     TH1F *JRcorrJZBSBem;
911     TH1F *JLcorrJZBSBem;
912     TH1F *JRcorrJZBSBeemm;
913     TH1F *JLcorrJZBSBeemm;
914    
915 buchmann 1.10 if(use_data==2 || overlay_signal) RcorrJZBeemmNoS = allsamples.Draw("RcorrJZBeemmNoS",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets&&largerzero,is_data, luminosity,false);
916 buchmann 1.1
917    
918     if(PlottingSetup::RestrictToMassPeak) {
919 buchmann 1.10 RcorrJZBSBem = allsamples.Draw("RcorrJZBSBem",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets&&largerzero,is_data, luminosity,use_signal);
920     LcorrJZBSBem = allsamples.Draw("LcorrJZBSBem",("-"+jzb).c_str(),nbins,low,hi, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets&&smallerzero,is_data, luminosity,use_signal);
921     RcorrJZBSBeemm = allsamples.Draw("RcorrJZBSBeemm",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets&&largerzero,is_data, luminosity,use_signal);
922     LcorrJZBSBeemm = allsamples.Draw("LcorrJZBSBeemm",("-"+jzb).c_str(),nbins,low,hi, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets&&smallerzero,is_data, luminosity,use_signal);
923    
924 buchmann 1.1 //these are for the ratio
925 buchmann 1.10 JRcorrJZBSBem = allsamples.Draw("JRcorrJZBSBem",jzb.c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets&&largerzero,is_data, luminosity,use_signal);
926     JLcorrJZBSBem = allsamples.Draw("JLcorrJZBSBem",("-"+jzb).c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets&&smallerzero,is_data, luminosity,use_signal);
927     JRcorrJZBSBeemm = allsamples.Draw("JRcorrJZBSBeemm",jzb.c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets&&largerzero,is_data, luminosity,use_signal);
928     JLcorrJZBSBeemm = allsamples.Draw("JLcorrJZBSBeemm",("-"+jzb).c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets&&smallerzero,is_data, luminosity,use_signal);
929 buchmann 1.1 }
930    
931     TH1F *lm4RcorrJZBeemm;
932 buchmann 1.10 if(overlay_signal || use_data == 2 || use_data == 1) lm4RcorrJZBeemm = allsamples.Draw("lm4RcorrJZBSBeemm",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets&&largerzero,is_data, luminosity,allsamples.FindSample("LM"));
933 buchmann 1.1
934     flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak ---- prediction changed.
935    
936     TH1F *Zjetspred = (TH1F*)LcorrJZBeemm->Clone("Zjetspred");
937     TH1F *TTbarpred = (TH1F*)RcorrJZBem->Clone("TTbarpred");
938    
939     TH1F *Bpred = (TH1F*)LcorrJZBeemm->Clone("Bpred");
940     TH1F *JBpred = (TH1F*)JLcorrJZBeemm->Clone("Bpred");
941 buchmann 1.23
942     TH1F *BpredSys = (TH1F*)LcorrJZBeemm->Clone("Bpredsys");
943     ClearHisto(BpredSys);
944     write_error(__FUNCTION__,"Implemented systematic error but NOT drawing it yet. Need to still implement drawing, particularly in ratio.");
945    
946 buchmann 1.1 if(PlottingSetup::RestrictToMassPeak) {
947     Bpred->Add(RcorrJZBem,1.0/3.);
948     Bpred->Add(LcorrJZBem,-1.0/3.);
949     Bpred->Add(RcorrJZBSBem,1.0/3.);
950     Bpred->Add(LcorrJZBSBem,-1.0/3.);
951     Bpred->Add(RcorrJZBSBeemm,1.0/3.);
952     Bpred->Add(LcorrJZBSBeemm,-1.0/3.);
953    
954     TTbarpred->Scale(1.0/3);
955     Zjetspred->Add(LcorrJZBem,-1.0/3.);
956     Zjetspred->Add(LcorrJZBSBem,-1.0/3.);
957     TTbarpred->Add(RcorrJZBSBem,1.0/3.);
958     Zjetspred->Add(LcorrJZBSBeemm,-1.0/3.);
959     TTbarpred->Add(RcorrJZBSBeemm,1.0/3.);
960    
961     //these are for the ratio
962     JBpred->Add(JRcorrJZBem,1.0/3.);
963     JBpred->Add(JLcorrJZBem,-1.0/3.);
964     JBpred->Add(JRcorrJZBSBem,1.0/3.);
965     JBpred->Add(JLcorrJZBSBem,-1.0/3.);
966     JBpred->Add(JRcorrJZBSBeemm,1.0/3.);
967     JBpred->Add(JLcorrJZBSBeemm,-1.0/3.);
968 buchmann 1.23
969     //Systematics:
970     AddSquared(BpredSys,LcorrJZBeemm,zjetsestimateuncertONPEAK*zjetsestimateuncertONPEAK);
971     AddSquared(BpredSys,RcorrJZBem,emuncertONPEAK*emuncertONPEAK*(1.0/9));
972     AddSquared(BpredSys,LcorrJZBem,emuncertONPEAK*emuncertONPEAK*(1.0/9));
973     AddSquared(BpredSys,RcorrJZBSBem,emsidebanduncertONPEAK*emsidebanduncertONPEAK*(1.0/9));
974     AddSquared(BpredSys,LcorrJZBSBem,emsidebanduncertONPEAK*emsidebanduncertONPEAK*(1.0/9));
975     AddSquared(BpredSys,RcorrJZBSBeemm,eemmsidebanduncertONPEAK*eemmsidebanduncertONPEAK*(1.0/9));
976     AddSquared(BpredSys,LcorrJZBSBeemm,eemmsidebanduncertONPEAK*eemmsidebanduncertONPEAK*(1.0/9));
977 buchmann 1.1 } else {
978     Bpred->Add(RcorrJZBem,1.0);
979     Bpred->Add(LcorrJZBem,-1.0);
980    
981     Zjetspred->Add(LcorrJZBem,-1.0);
982    
983     //these are for the ratio
984     JBpred->Add(JRcorrJZBem,1.0);
985     JBpred->Add(JLcorrJZBem,-1.0);
986 buchmann 1.23
987     //Systematics
988     AddSquared(BpredSys,LcorrJZBeemm,zjetsestimateuncertOFFPEAK*zjetsestimateuncertOFFPEAK);
989     AddSquared(BpredSys,RcorrJZBem,emuncertOFFPEAK*emuncertOFFPEAK);
990     AddSquared(BpredSys,LcorrJZBem,emuncertOFFPEAK*emuncertOFFPEAK);
991    
992 buchmann 1.1 }
993    
994 buchmann 1.23 SQRT(BpredSys);
995    
996 buchmann 1.1
997     flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak ---- prediction changed
998     TH1F *Tpred = (TH1F*)RcorrJZBem->Clone("Bpred");
999     Tpred->Add(LcorrJZBem,-1.0);
1000     if(PlottingSetup::RestrictToMassPeak) {
1001     Tpred->Add(RcorrJZBSBem,1.0);
1002     Tpred->Add(LcorrJZBSBem,-1.0);
1003     Tpred->Add(RcorrJZBSBeemm,1.0);
1004     Tpred->Add(LcorrJZBSBeemm,-1.0);
1005     }
1006    
1007     globalcanvas->cd();
1008     globalcanvas->SetLogy(1);
1009    
1010     RcorrJZBeemm->SetMarkerStyle(20);
1011     RcorrJZBeemm->GetXaxis()->SetRangeUser(0,high);
1012     RcorrJZBeemm->SetMinimum(0.1);
1013    
1014     Bpred->SetLineColor(kRed);
1015     Bpred->SetStats(0);
1016    
1017     int versok=false;
1018     if(gROOT->GetVersionInt()>=53000) versok=true;
1019    
1020    
1021     if ( overlay_signal || use_data==2 ) lm4RcorrJZBeemm->SetLineColor(TColor::GetColor("#088A08"));
1022     RcorrJZBeemm->SetMarkerSize(DataMarkerSize);
1023    
1024     TLegend *legBpred = make_legend("",0.6,0.55);
1025     TLegend *legBpred2 = make_legend("",0.6,0.55);
1026    
1027    
1028     vector<TF1*> analytical_function;
1029     TPad *kinpad = new TPad("kinpad","kinpad",0,0,1,1);
1030     kinpad->cd();
1031     kinpad->SetLogy(1);
1032    
1033     string Bpredsaveas="Bpred_Data";
1034     blankback->SetMaximum(5*RcorrJZBeemm->GetMaximum());
1035     blankback->SetMinimum(0.1);
1036     if(use_data!=1) blankback->SetMinimum(0.1);
1037     blankback->Draw();
1038     if(use_data==1)
1039     {
1040 buchmann 1.14 //Bpred->SetLineWidth(3); //paper style.overruled.
1041     //lm4RcorrJZBeemm->SetLineWidth(3); //paper style.overruled.
1042 buchmann 1.1 analytical_function = do_extended_fit_to_plot(Bpred,Tpred,LcorrJZBeemm,LcorrJZBem,is_data);
1043     kinpad->cd();//necessary because the extended fit function creates its own canvas
1044     RcorrJZBeemm->Draw("e1x0,same");
1045    
1046     Bpred->Draw("hist,same");
1047 buchmann 1.14 //analytical_function[0]->Draw("same"); analytical_function[1]->Draw("same");analytical_function[2]->Draw("same");
1048 buchmann 1.1 RcorrJZBeemm->Draw("e1x0,same");//HAVE IT ON TOP!
1049 buchmann 1.10 lm4RcorrJZBeemm->Draw("hist,same");
1050 buchmann 1.1 legBpred->AddEntry(RcorrJZBeemm,"observed","p");
1051     legBpred->AddEntry(Bpred,"predicted","l");
1052 buchmann 1.24 // legBpred->AddEntry(analytical_function[1],"predicted fit","l");
1053     // legBpred->AddEntry(analytical_function[2],"stat. uncert.","l");
1054     legBpred->AddEntry(lm4RcorrJZBeemm,"LM4","l");
1055 buchmann 1.1 legBpred->Draw();
1056     DrawPrelim();
1057    
1058     //this plot shows what the prediction is composed of
1059 buchmann 1.10 TPad *predcomppad = new TPad("predcomppad","predcomppad",0,0,1,1);
1060     float CurrentBpredLineWidth=Bpred->GetLineWidth();
1061     Bpred->SetLineWidth(2);
1062     predcomppad->cd();
1063     predcomppad->SetLogy(1);
1064    
1065     TH1F *jzbnegative = (TH1F*)LcorrJZBeemm->Clone("jzbnegative");
1066     TH1F *sidebandsemu = (TH1F*)Bpred->Clone("sidebandsemu");
1067     sidebandsemu->Add(jzbnegative,-1);
1068    
1069     jzbnegative->SetFillColor(allsamples.GetColor((allsamples.FindSample("DY"))[0]));
1070     jzbnegative->SetLineColor(allsamples.GetColor((allsamples.FindSample("DY"))[0]));
1071     sidebandsemu->SetLineColor(allsamples.GetColor((allsamples.FindSample("TTJets"))[0]));
1072     sidebandsemu->SetFillColor(allsamples.GetColor((allsamples.FindSample("TTJets"))[0]));
1073    
1074 buchmann 1.1 THStack predcomposition("predcomposition","prediction composition");
1075 buchmann 1.10 predcomposition.Add(sidebandsemu);
1076     predcomposition.Add(jzbnegative);
1077 buchmann 1.1 blankback->Draw();
1078     RcorrJZBeemm->Draw("e1x0,same");
1079     predcomposition.Draw("histo,same");//
1080     Bpred->Draw("hist,same");
1081 buchmann 1.10 // analytical_function[0]->Draw("same"); analytical_function[1]->Draw("same");analytical_function[2]->Draw("same");
1082 buchmann 1.1 RcorrJZBeemm->Draw("e1x0,same");//HAVE IT ON TOP!
1083 buchmann 1.14 // lm4RcorrJZBeemm->SetLineColor(kOrange+1);
1084 buchmann 1.10 lm4RcorrJZBeemm->SetLineWidth(2);
1085 buchmann 1.14 //lm4RcorrJZBeemm->SetLineWidth(2); // paper style. overruled.
1086 buchmann 1.10 lm4RcorrJZBeemm->Draw("histo,same");
1087 buchmann 1.1 DrawPrelim();
1088 buchmann 1.10 TLegend *speciallegBpred = make_legend("",0.45,0.55);
1089 buchmann 1.14 //TLegend *speciallegBpred = make_legend("",0.35,0.55); // paper style. overruled.
1090 buchmann 1.10 speciallegBpred->AddEntry(RcorrJZBeemm,"Data","pl");
1091     speciallegBpred->AddEntry(Bpred,"Total background","l");
1092     speciallegBpred->AddEntry(jzbnegative,"JZB<0 (data)","f");
1093     speciallegBpred->AddEntry(sidebandsemu,"Sidebands/e#mu (data)","f");
1094     // speciallegBpred->AddEntry(lm4RcorrJZBeemmC,"LM4","l");
1095     speciallegBpred->AddEntry(lm4RcorrJZBeemm,"LM4","l");
1096     speciallegBpred->Draw();
1097 buchmann 1.24 save_with_ratio(JRcorrJZBeemm,JBpred,predcomppad,subdir+"Bpred_Data_____PredictionComposition",true,true,"data/pred",BpredSys);
1098 buchmann 1.10
1099     TCanvas *specialcanv = new TCanvas("specialcanv","specialcanv");
1100 buchmann 1.1 THStack kostack = allsamples.DrawStack("RcorrJZBeemm",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,!is_data, luminosity,use_signal);
1101     blankback->Draw();
1102     kostack.Draw("same");
1103     Bpred->Draw("hist,same");
1104     analytical_function[0]->Draw("same"); analytical_function[1]->Draw("same");analytical_function[2]->Draw("same");
1105     RcorrJZBeemm->Draw("e1x0,same");//HAVE IT ON TOP!
1106     legBpred->Draw();
1107     DrawPrelim();
1108 buchmann 1.12 CompleteSave(specialcanv,subdir+"Bpred_Data_____PredictionCompositioninMC");
1109 buchmann 1.10 Bpred->SetLineWidth((int)CurrentBpredLineWidth);
1110 buchmann 1.1
1111 buchmann 1.10 delete speciallegBpred;
1112 buchmann 1.1 delete Zjetspred;
1113     delete TTbarpred;
1114    
1115     kinpad->cd();
1116     }
1117     if(use_data==0) {
1118     RcorrJZBeemm->Draw("e1x0,same");
1119 buchmann 1.14 //Bpred->SetLineWidth(3); // paper style. overruled.
1120 buchmann 1.1 Bpred->Draw("hist,same");
1121     RcorrJZBeemm->Draw("e1x0,same");//HAVE IT ON TOP!
1122 buchmann 1.10 legBpred->AddEntry(RcorrJZBeemm,"MC true","p");
1123 buchmann 1.18 legBpred->AddEntry(Bpred,"MC predicted","l");
1124 buchmann 1.1 if(versok) legBpred->AddEntry((TObject*)0,"",""); // Just for alignment // causes seg fault on root v5.18
1125     if(versok) legBpred->AddEntry((TObject*)0,"",""); // causes seg fault on root v5.18
1126     if ( overlay_signal ) legBpred->AddEntry(lm4RcorrJZBeemm,"LM4","l");
1127     legBpred->Draw();
1128     DrawMCPrelim();
1129     Bpredsaveas="Bpred_MC";
1130     // CompleteSave(globalcanvas,"Bpred_MC"); // done below in save_with_ratio
1131     }
1132     if(use_data==2) {
1133     RcorrJZBeemm->Draw("e1x0,same");
1134 buchmann 1.14 //Bpred->SetLineWidth(3); // paper style. overruled.
1135 buchmann 1.1 Bpred->Draw("hist,same");
1136     RcorrJZBeemm->Draw("e1x0,same");//HAVE IT ON TOP!
1137 buchmann 1.10 legBpred->AddEntry(RcorrJZBeemm,"MC true","p");
1138 buchmann 1.1 legBpred->AddEntry(Bpred,"MC predicted","l");
1139 buchmann 1.10 legBpred2->AddEntry(RcorrJZBeemm,"MC true","p");
1140 buchmann 1.1 legBpred2->AddEntry(Bpred,"MC predicted","l");
1141     {
1142     if(versok) legBpred->AddEntry((TObject*)0,"",""); // Just for alignment // causes seg fault on root v5.18 --> now only allowed for root >=v5.30
1143     if(versok) legBpred->AddEntry((TObject*)0,"",""); // causes seg fault on root v5.18 --> now only allowed for root >=v5.30
1144     legBpred->Draw();
1145     DrawMCPrelim();
1146     Bpredsaveas="Bpred_MCwithS";
1147     // CompleteSave(globalcanvas,"Bpred_MCwithS"); // done below in save_with_ratio
1148     }
1149     {
1150 buchmann 1.14 //lm4RcorrJZBeemm->SetLineWidth(3); //paper style. overruled.
1151     //RcorrJZBeemmNoS->SetLineWidth(3); //paper style. overruled.
1152     //lm4RcorrJZBeemm->SetLineStyle(2); //paper style. overruled.
1153     //RcorrJZBeemmNoS->SetLineStyle(3); //paper style. overruled.
1154     //lm4RcorrJZBeemm->SetLineColor(kOrange+1); //paper style. overruled.
1155    
1156 buchmann 1.1 RcorrJZBeemmNoS->SetLineStyle(2);
1157     legBpred2->AddEntry(RcorrJZBeemmNoS,"MC B","l");
1158     legBpred2->AddEntry(lm4RcorrJZBeemm,"MC S","l");
1159     legBpred2->Draw();
1160     RcorrJZBeemmNoS->SetLineColor(TColor::GetColor("#61210B"));
1161     RcorrJZBeemmNoS->Draw("histo,same");
1162     RcorrJZBeemm->Draw("e1x0,same");//HAVE IT ON TOP!
1163     lm4RcorrJZBeemm->Draw("histo,same");
1164     DrawMCPrelim();
1165     Bpredsaveas="Bpred_MCwithS__plus";
1166     // CompleteSave(globalcanvas,"Bpred_MCwithS__plus"); // done below in save_with_ratio
1167     }
1168     }
1169    
1170    
1171     //3rd last argument: do special bpred ratio, 2nd last argument: extended range!, last: y-axis title
1172     string ytitle("ratio");
1173     if ( use_data==1 ) ytitle = "data/pred";
1174 buchmann 1.8 //save_with_ratio(JRcorrJZBeemm,JBpred,kinpad,Bpredsaveas,true,use_data!=1,ytitle);
1175 buchmann 1.24 save_with_ratio(JRcorrJZBeemm,JBpred,kinpad,subdir+Bpredsaveas,true,true,ytitle);//not extending the y range anymore up to 4
1176 buchmann 1.1
1177    
1178     flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
1179 buchmann 1.24 // The part below is meaningless for the offpeak analysis (it's a comparison of the different estimates but there is but one estimate!)
1180 buchmann 1.11 if(PlottingSetup::RestrictToMassPeak) {
1181     TH1F *Bpredem = (TH1F*)LcorrJZBeemm->Clone("Bpredem");
1182     Bpredem->Add(RcorrJZBem);
1183     Bpredem->Add(LcorrJZBem,-1);
1184     TH1F *BpredSBem = (TH1F*)LcorrJZBeemm->Clone("BpredSBem");
1185     BpredSBem->Add(RcorrJZBSBem);
1186     Bpred->Add(LcorrJZBSBem,-1);
1187     TH1F *BpredSBeemm = (TH1F*)LcorrJZBeemm->Clone("BpredSBeemm");
1188     BpredSBeemm->Add(RcorrJZBSBeemm);
1189     BpredSBeemm->Add(LcorrJZBSBeemm,-1.0);
1190     globalcanvas->cd();
1191     globalcanvas->SetLogy(1);
1192    
1193     RcorrJZBeemm->SetMarkerStyle(20);
1194     RcorrJZBeemm->GetXaxis()->SetRangeUser(0,high);
1195     blankback->Draw();
1196     RcorrJZBeemm->Draw("e1x0,same");
1197     RcorrJZBeemm->SetMarkerSize(DataMarkerSize);
1198    
1199     Bpredem->SetLineColor(kRed+1);
1200     Bpredem->SetStats(0);
1201     Bpredem->Draw("hist,same");
1202    
1203     BpredSBem->SetLineColor(kGreen+2);//TColor::GetColor("#0B6138"));
1204     BpredSBem->SetLineStyle(2);
1205     BpredSBem->Draw("hist,same");
1206    
1207     BpredSBeemm->SetLineColor(kBlue+1);
1208     BpredSBeemm->SetLineStyle(3);
1209     BpredSBeemm->Draw("hist,same");
1210     RcorrJZBeemm->Draw("e1x0,same");//HAVE IT ON TOP!
1211    
1212     TLegend *legBpredc = make_legend("",0.6,0.55);
1213     if(use_data==1)
1214     {
1215     legBpredc->AddEntry(RcorrJZBeemm,"observed","p");
1216     legBpredc->AddEntry(Bpredem,"OFZP","l");
1217     legBpredc->AddEntry(BpredSBem,"OFSB","l");
1218     legBpredc->AddEntry(BpredSBeemm,"SFSB","l");
1219     legBpredc->Draw();
1220 buchmann 1.12 CompleteSave(globalcanvas,subdir+"Bpred_Data_comparison");
1221 buchmann 1.11 }
1222     if(use_data==0) {
1223     legBpredc->AddEntry(RcorrJZBeemm,"MC true","p");
1224     legBpredc->AddEntry(Bpredem,"MC OFZP","l");
1225     legBpredc->AddEntry(BpredSBem,"MC OFSB","l");
1226     legBpredc->AddEntry(BpredSBeemm,"MC SFSB","l");
1227     legBpredc->Draw();
1228 buchmann 1.12 CompleteSave(globalcanvas,subdir+"Bpred_MC_comparison");
1229 buchmann 1.11 }
1230     if(use_data==2) {
1231     legBpredc->AddEntry(RcorrJZBeemm,"MC true","p");
1232     legBpredc->AddEntry(Bpredem,"MC OFZP","l");
1233     legBpredc->AddEntry(BpredSBem,"MC OFSB","l");
1234     legBpredc->AddEntry(BpredSBeemm,"MC SFSB","l");
1235     if ( overlay_signal ) legBpred->AddEntry(lm4RcorrJZBeemm,"LM4","l");
1236     legBpredc->Draw();
1237 buchmann 1.12 CompleteSave(globalcanvas,subdir+"Bpred_MCwithS_comparison");
1238 buchmann 1.11 }
1239     }
1240 buchmann 1.1
1241     delete RcorrJZBeemm;
1242     delete LcorrJZBeemm;
1243     delete RcorrJZBem;
1244     delete LcorrJZBem;
1245    
1246     delete JRcorrJZBeemm;
1247     delete JLcorrJZBeemm;
1248     delete JRcorrJZBem;
1249     delete JLcorrJZBem;
1250    
1251     delete blankback;
1252    
1253     if(PlottingSetup::RestrictToMassPeak) {
1254     delete RcorrJZBSBem;
1255     delete LcorrJZBSBem;
1256     delete RcorrJZBSBeemm;
1257     delete LcorrJZBSBeemm;
1258    
1259     delete JRcorrJZBSBem;
1260     delete JLcorrJZBSBem;
1261     delete JRcorrJZBSBeemm;
1262     delete JLcorrJZBSBeemm;
1263     }
1264     if(overlay_signal || use_data==2) delete lm4RcorrJZBeemm;
1265 buchmann 1.10 switch_overunderflow(false);
1266 buchmann 1.1 }
1267    
1268     void do_prediction_plots(string mcjzb, string datajzb, float DataSigma, float MCSigma, bool overlay_signal ) {
1269     TCanvas *globalcanvas = new TCanvas("globalcanvas","Prediction Canvas");
1270     do_prediction_plot(datajzb,globalcanvas,DataSigma,jzbHigh ,data,overlay_signal);
1271 buchmann 1.14 if ( !PlottingSetup::Approved ) {
1272     do_prediction_plot(mcjzb,globalcanvas,MCSigma,jzbHigh ,mc,overlay_signal);
1273     do_prediction_plot(mcjzb,globalcanvas,MCSigma,jzbHigh ,mcwithsignal,overlay_signal);
1274     } else {
1275     write_info(__FUNCTION__,"You set approved to true, therefore not producing prediction/observation plots for MC with and without signal.");
1276 buchmann 1.15 }
1277 buchmann 1.1 }
1278    
1279     void do_ratio_plot(int is_data,vector<float> binning, string jzb, TCanvas *can, float high=-9999) {
1280     bool do_data=0;
1281     bool dosignal=0;
1282     if(is_data==1) do_data=1;
1283     if(is_data==2) dosignal=1;
1284     TH1F *RcorrJZBeemm = allsamples.Draw("RcorrJZBeemm",jzb.c_str(),binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,do_data, luminosity,dosignal);
1285     TH1F *LcorrJZBeemm = allsamples.Draw("LcorrJZBeemm",("-"+jzb).c_str(),binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,do_data, luminosity,dosignal);
1286     TH1F *RcorrJZBem = allsamples.Draw("RcorrJZBem",jzb.c_str(),binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,do_data, luminosity,dosignal);
1287     TH1F *LcorrJZBem = allsamples.Draw("LcorrJZBem",("-"+jzb).c_str(),binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,do_data, luminosity,dosignal);
1288    
1289     TH1F *RcorrJZBSBem;
1290     TH1F *LcorrJZBSBem;
1291     TH1F *RcorrJZBSBeemm;
1292     TH1F *LcorrJZBSbeemm;
1293    
1294     if(PlottingSetup::RestrictToMassPeak) {
1295     RcorrJZBSBem = allsamples.Draw("RcorrJZBSbem",jzb.c_str(),binning, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets,do_data, luminosity,dosignal);
1296     LcorrJZBSBem = allsamples.Draw("LcorrJZBSbem",("-"+jzb).c_str(),binning, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets,do_data, luminosity,dosignal);
1297     RcorrJZBSBeemm = allsamples.Draw("RcorrJZBSbeemm",jzb.c_str(),binning, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets,do_data, luminosity,dosignal);
1298     LcorrJZBSbeemm = allsamples.Draw("LcorrJZBSbeemm",("-"+jzb).c_str(),binning, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets,do_data, luminosity,dosignal);
1299     }
1300    
1301    
1302    
1303    
1304     flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
1305     TH1F *Bpred = (TH1F*)LcorrJZBeemm->Clone("Bpred");
1306     if(PlottingSetup::RestrictToMassPeak) {
1307     Bpred->Add(RcorrJZBem,1.0/3);
1308     Bpred->Add(LcorrJZBem,-1.0/3);
1309     Bpred->Add(RcorrJZBSBem,1.0/3);
1310     Bpred->Add(LcorrJZBSBem,-1.0/3);
1311     Bpred->Add(RcorrJZBSBeemm,1.0/3);
1312     Bpred->Add(LcorrJZBSbeemm,-1.0/3);
1313     } else {
1314     Bpred->Add(RcorrJZBem,1.0);
1315     Bpred->Add(LcorrJZBem,-1.0);
1316     }
1317    
1318     can->cd();
1319     can->SetLogy(0);
1320     Bpred->SetLineColor(kRed);
1321     Bpred->SetStats(0);
1322     if(high>0) Bpred->GetXaxis()->SetRangeUser(0,high);
1323     TH1F *JZBratioforfitting=(TH1F*)RcorrJZBeemm->Clone("JZBratioforfitting");
1324     JZBratioforfitting->Divide(Bpred);
1325     TGraphAsymmErrors *JZBratio = histRatio(RcorrJZBeemm,Bpred,is_data,binning,false);
1326    
1327    
1328     JZBratio->SetTitle("");
1329     JZBratio->GetYaxis()->SetRangeUser(0.0,9.0);
1330     // if(is_data==1) JZBratio->GetXaxis()->SetRangeUser(0,jzbHigh );
1331    
1332     TF1 *pol0 = new TF1("pol0","[0]",0,1000);
1333     TF1 *pol0d = new TF1("pol0","[0]",0,1000);
1334     // straightline_fit->SetParameter(0,1);
1335     JZBratioforfitting->Fit(pol0,"Q0R","",0,30);
1336     pol0d->SetParameter(0,pol0->GetParameter(0));
1337    
1338     JZBratio->GetYaxis()->SetTitle("Observed/Predicted");
1339     JZBratio->GetXaxis()->SetTitle("JZB [GeV]");
1340     if ( high>0 ) JZBratio->GetXaxis()->SetRangeUser(0.0,high);
1341     JZBratio->GetYaxis()->SetNdivisions(519);
1342     JZBratio->GetYaxis()->SetRangeUser(0.0,4.0);
1343     JZBratio->GetYaxis()->CenterTitle();
1344     JZBratio->GetXaxis()->CenterTitle();
1345     JZBratio->SetMarkerSize(DataMarkerSize);
1346     JZBratio->Draw("AP");
1347     /////----------------------------
1348     TPaveText *writeresult = new TPaveText(0.15,0.78,0.49,0.91,"blNDC");
1349     writeresult->SetFillStyle(4000);
1350     writeresult->SetFillColor(kWhite);
1351     writeresult->SetTextFont(42);
1352     writeresult->SetTextSize(0.03);
1353     writeresult->SetTextAlign(12);
1354     ostringstream jzb_agreement_data_text;
1355     jzb_agreement_data_text<< setprecision(2) << "mean =" << pol0->GetParameter(0) << " #pm " << setprecision(1) << pol0->GetParError(0);
1356     if(is_data==1) fitresultconstdata=pol0->GetParameter(0);// data
1357     if(is_data==0) fitresultconstmc=pol0->GetParameter(0); // monte carlo, no signal
1358     /* if(is_data) writeresult->AddText("Data closure test");
1359     else writeresult->AddText("MC closure test");
1360     */
1361     writeresult->AddText(jzb_agreement_data_text.str().c_str());
1362     // writeresult->Draw("same");
1363     // pol0d->Draw("same");
1364     TF1 *topline = new TF1("","1.5",0,1000);
1365     TF1 *bottomline = new TF1("","0.5",0,1000);
1366     topline->SetLineColor(kBlue);
1367     topline->SetLineStyle(2);
1368     bottomline->SetLineColor(kBlue);
1369     bottomline->SetLineStyle(2);
1370     // topline->Draw("same");
1371     // bottomline->Draw("same");
1372     TF1 *oneline = new TF1("","1.0",0,1000);
1373     oneline->SetLineColor(kBlue);
1374     oneline->SetLineStyle(1);
1375     oneline->Draw("same");
1376     TLegend *phony_leg = make_legend("ratio",0.6,0.55,false);//this line is just to have the default CMS Preliminary (...) on the canvas as well.
1377     if(is_data==1) DrawPrelim();
1378     else DrawMCPrelim();
1379     TLegend *leg = new TLegend(0.55,0.75,0.89,0.89);
1380     leg->SetTextFont(42);
1381     leg->SetTextSize(0.04);
1382     // if(is_data==1) leg->SetHeader("Ratio (data)");
1383     // else leg->SetHeader("Ratio (MC)");
1384    
1385     TString MCtitle("MC ");
1386     if (is_data==1) MCtitle = "";
1387    
1388     leg->SetFillStyle(4000);
1389     leg->SetFillColor(kWhite);
1390     leg->SetTextFont(42);
1391     // leg->AddEntry(topline,"+20\% sys envelope","l");
1392     leg->AddEntry(JZBratio,MCtitle+"obs / "+MCtitle+"pred","p");
1393     leg->AddEntry(oneline,"ratio = 1","l");
1394     // leg->AddEntry(pol0d,"fit in [0,30] GeV","l");
1395     // leg->AddEntry(bottomline,"#pm50% envelope","l");
1396    
1397    
1398     //leg->Draw("same"); // no longer drawing legend
1399    
1400     if(is_data==1) CompleteSave(can, "jzb_ratio_data");
1401     if(is_data==0) CompleteSave(can, "jzb_ratio_mc");
1402     if(is_data==2) CompleteSave(can, "jzb_ratio_mc_BandS");//special case, MC with signal!
1403    
1404     delete RcorrJZBeemm;
1405     delete LcorrJZBeemm;
1406     delete RcorrJZBem;
1407     delete LcorrJZBem;
1408    
1409     delete RcorrJZBSBem;
1410     delete LcorrJZBSBem;
1411     delete RcorrJZBSBeemm;
1412     delete LcorrJZBSbeemm;
1413     }
1414    
1415     void do_ratio_plots(string mcjzb,string datajzb,vector<float> ratio_binning) {
1416     TCanvas *globalc = new TCanvas("globalc","Ratio Plot Canvas");
1417     globalc->SetLogy(0);
1418    
1419     do_ratio_plot(mc,ratio_binning,mcjzb,globalc, jzbHigh );
1420     do_ratio_plot(data,ratio_binning,datajzb,globalc, jzbHigh );
1421     do_ratio_plot(mcwithsignal,ratio_binning,mcjzb,globalc, jzbHigh );
1422     }
1423    
1424     string give_jzb_expression(float peak, int type) {
1425     stringstream val;
1426     if(type==data) {
1427     if(peak<0) val << jzbvariabledata << "+" << TMath::Abs(peak);
1428     if(peak>0) val << jzbvariabledata << "-" << TMath::Abs(peak);
1429     if(peak==0) val << jzbvariabledata;
1430     }
1431     if(type==mc) {
1432     if(peak<0) val << jzbvariablemc << "+" << TMath::Abs(peak);
1433     if(peak>0) val << jzbvariablemc << "-" << TMath::Abs(peak);
1434     if(peak==0) val << jzbvariablemc;
1435     }
1436     return val.str();
1437     }
1438    
1439    
1440     void lepton_comparison_plots() {
1441     Float_t ymin = 1.e-5, ymax = 0.25;
1442     TCanvas *can = new TCanvas("can","Lepton Comparison Canvas");
1443     can->SetLogy(1);
1444 buchmann 1.3 TH1F *eemc = allsamples.Draw("eemc","mll",50,50,150, "mll [GeV]", "events", cutOSSF&&cutnJets&&"(id1==0)",mc, luminosity,allsamples.FindSample("/DY"));
1445     TH1F *mmmc = allsamples.Draw("mmmc","mll",50,50,150, "mll [GeV]", "events", cutOSSF&&cutnJets&&"(id1==1)",mc, luminosity,allsamples.FindSample("/DY"));
1446 buchmann 1.1 eemc->SetLineColor(kBlue);
1447     mmmc->SetLineColor(kRed);
1448     eemc->SetMinimum(0.1);
1449     eemc->SetMaximum(10*eemc->GetMaximum());
1450     eemc->Draw("histo");
1451     mmmc->Draw("histo,same");
1452     TLegend *leg = make_legend();
1453     leg->AddEntry(eemc,"ZJets->ee (MC)","l");
1454     leg->AddEntry(mmmc,"ZJets->#mu#mu (MC)","l");
1455     leg->Draw("same");
1456     CompleteSave(can, "lepton_comparison/mll_effratio_mc");
1457    
1458     TH1F *eed = allsamples.Draw("eed","mll",50,50,150, "mll [GeV]", "events", cutOSSF&&cutnJets&&"(id1==0)",data, luminosity);
1459     TH1F *mmd = allsamples.Draw("mmd","mll",50,50,150, "mll [GeV]", "events", cutOSSF&&cutnJets&&"(id1==1)",data, luminosity);
1460     eed->SetLineColor(kBlue);
1461     mmd->SetLineColor(kRed);
1462     eed->SetMinimum(0.1);
1463     eed->SetMaximum(10*eed->GetMaximum());
1464     eed->Draw("histo");
1465     mmd->Draw("histo,same");
1466     TLegend *leg2 = make_legend();
1467     leg2->AddEntry(eed,"ZJets->ee (data)","l");
1468     leg2->AddEntry(mmd,"ZJets->#mu#mu (data)","l");
1469     leg2->Draw();
1470     CompleteSave(can, "lepton_comparison/mll_effratio_data");
1471    
1472     TH1F *jeed = allsamples.Draw("jeed",jzbvariabledata, int((jzbHigh+150)/5),-150,jzbHigh , "JZB [GeV]", "events", cutOSSF&&cutnJets&&"(id1==0)",data, luminosity);
1473     TH1F *jmmd = allsamples.Draw("jmmd",jzbvariabledata, int((jzbHigh+150)/5),-150,jzbHigh , "JZB [GeV]", "events", cutOSSF&&cutnJets&&"(id1==1)",data, luminosity);
1474     TH1F *jeemmd = allsamples.Draw("jeemmd",jzbvariabledata,int((jzbHigh+150)/5),-150,jzbHigh , "JZB [GeV]", "events", cutOSSF&&cutnJets,data, luminosity);
1475     dout << "ee : " << jeed->GetMean() << "+/-" << jeed->GetMeanError() << endl;
1476     dout << "ee : " << jmmd->GetMean() << "+/-" << jmmd->GetMeanError() << endl;
1477     dout << "eemd : " << jeemmd->GetMean() << "+/-" << jeemmd->GetMeanError() << endl;
1478     jeemmd->SetLineColor(kBlack);
1479     jeemmd->SetMarkerStyle(25);
1480     jeed->SetLineColor(kBlue);
1481     jmmd->SetLineColor(kRed);
1482     jeed->SetMinimum(0.1);
1483     jeed->SetMaximum(10*eed->GetMaximum());
1484     TH1* njeemmd = jeemmd->DrawNormalized();
1485     njeemmd->SetMinimum(ymin);
1486     njeemmd->SetMaximum(ymax);
1487    
1488     jeed->DrawNormalized("histo,same");
1489     jmmd->DrawNormalized("histo,same");
1490     jeemmd->DrawNormalized("same");
1491     TLegend *jleg2 = make_legend(" ");
1492     jleg2->AddEntry(jeemmd,"ee and #mu#mu","p");
1493     jleg2->AddEntry(jeed,"ee","l");
1494     jleg2->AddEntry(jmmd,"#mu#mu","l");
1495     jleg2->Draw();
1496     CompleteSave(can,"lepton_comparison/jzb_effratio_data");
1497    
1498     TPad *eemmpad = new TPad("eemmpad","eemmpad",0,0,1,1);
1499     eemmpad->cd();
1500     eemmpad->SetLogy(1);
1501     jeed->Draw("histo");
1502     jmmd->Draw("histo,same");
1503     TLegend *eemmlegend = make_legend(" ");
1504     eemmlegend->AddEntry(jeed,"ee","l");
1505     eemmlegend->AddEntry(jmmd,"#mu#mu","l");
1506     eemmlegend->Draw();
1507     DrawPrelim();
1508     save_with_ratio(jeed,jmmd,eemmpad->cd(),"lepton_comparison/jzb_Comparing_ee_mm_data");
1509    
1510 buchmann 1.3 TH1F *zjeed = allsamples.Draw("zjeed",jzbvariablemc, int((jzbHigh+150)/5),-150,jzbHigh , "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets&&"(id1==0)",mc, luminosity,allsamples.FindSample("/DY"));
1511     TH1F *zjmmd = allsamples.Draw("zjmmd",jzbvariablemc, int((jzbHigh+150)/5),-150,jzbHigh , "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets&&"(id1==1)",mc, luminosity,allsamples.FindSample("/DY"));
1512     TH1F *zjeemmd = allsamples.Draw("zjeemmd",jzbvariablemc,int((jzbHigh+150)/5),-150,jzbHigh , "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets, mc, luminosity,allsamples.FindSample("/DY"));
1513 buchmann 1.1 dout << "Z+Jets ee : " << zjeed->GetMean() << "+/-" << zjeed->GetMeanError() << endl;
1514     dout << "Z+Jets ee : " << zjmmd->GetMean() << "+/-" << zjmmd->GetMeanError() << endl;
1515     dout << "Z+Jets eemd : " << zjeemmd->GetMean() << "+/-" << zjeemmd->GetMeanError() << endl;
1516     zjeemmd->SetLineColor(kBlack);
1517     zjeemmd->SetMarkerStyle(25);
1518     zjeed->SetLineColor(kBlue);
1519     zjmmd->SetLineColor(kRed);
1520     zjeed->SetMinimum(0.1);
1521     zjeed->SetMaximum(10*eed->GetMaximum());
1522    
1523     TH1* nzjeemmd = zjeemmd->DrawNormalized();
1524     nzjeemmd->SetMinimum(ymin);
1525     nzjeemmd->SetMaximum(ymax);
1526     zjeed->DrawNormalized("histo,same");
1527     zjmmd->DrawNormalized("histo,same");
1528     zjeemmd->DrawNormalized("same");
1529     TLegend *zjleg2 = make_legend("Z+jets MC");
1530     zjleg2->AddEntry(jeemmd,"ee and #mu#mu","p");
1531     zjleg2->AddEntry(jeed,"ee","l");
1532     zjleg2->AddEntry(jmmd,"#mu#mu","l");
1533     zjleg2->Draw();
1534     CompleteSave(can,"lepton_comparison/jzb_effratio_ZJets");
1535    
1536     TH1F *ld = allsamples.Draw("ld","mll",50,50,150, "mll [GeV]", "events", cutOSSF&&cutnJets,data, luminosity);
1537     ld->DrawNormalized("e1");
1538     eed->DrawNormalized("histo,same");
1539     mmd->DrawNormalized("histo,same");
1540     TLegend *leg3 = make_legend();
1541     leg3->AddEntry(ld,"ZJets->ll (data)","p");
1542     leg3->AddEntry(eed,"ZJets->ee (data)","l");
1543     leg3->AddEntry(mmd,"ZJets->#mu#mu (data)","l");
1544     leg3->Draw();
1545     CompleteSave(can,"lepton_comparison/mll_effratio_data__all_compared");
1546     /*
1547     TH1F *jzbld = allsamples.Draw("jzbld",jzbvariable,75,-150,150, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data, luminosity);
1548     TH1F *jzbemd = allsamples.Draw("jzbemd",jzbvariable,75,-150,150, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data, luminosity);
1549     */
1550     TH1F *jzbld = allsamples.Draw("jzbld",jzbvariabledata,int((jzbHigh+110)/5),-110,jzbHigh , "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data, luminosity);
1551     TH1F *jzbemd = allsamples.Draw("jzbemd",jzbvariabledata,int((jzbHigh+110)/5),-110,jzbHigh , "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data, luminosity);
1552     jzbld->SetMarkerColor(kBlack);
1553     jzbld->SetMarkerStyle(26);
1554     jzbemd->SetMarkerStyle(25);
1555     jzbemd->SetMarkerColor(kRed);
1556     jzbemd->SetLineColor(kRed);
1557     jzbld->SetMinimum(0.35);
1558     jzbld->Draw("e1");
1559     jzbemd->Draw("e1,same");
1560     TLegend *leg4 = make_legend();
1561     leg4->AddEntry(jzbld,"SFZP","p");
1562     leg4->AddEntry(jzbemd,"OFZP","p");
1563     leg4->AddEntry((TObject*)0,"",""); //causes segmentation violation
1564     leg4->AddEntry((TObject*)0,"",""); //causes segmentation violation
1565     leg4->Draw();
1566     CompleteSave(can,"lepton_comparison/jzb_eemumu_emu_data");
1567    
1568     TH1F *ttbarjzbld = allsamples.Draw("ttbarjzbld",jzbvariablemc,int((jzbHigh+150)/5),-150,jzbHigh, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data, luminosity,allsamples.FindSample("TTJet"));
1569     TH1F *ttbarjzbemd = allsamples.Draw("ttbarjzbemd",jzbvariablemc,int((jzbHigh+150)/5),-150,jzbHigh, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data, luminosity,allsamples.FindSample("TTJet"));
1570     ttbarjzbld->SetLineColor(allsamples.GetColor("TTJet"));
1571     ttbarjzbemd->SetLineColor(allsamples.GetColor("TTJet"));
1572     ttbarjzbld->Draw("histo");
1573     ttbarjzbemd->SetLineStyle(2);
1574     ttbarjzbemd->Draw("histo,same");
1575     TLegend *leg5 = make_legend();
1576     leg5->AddEntry(ttbarjzbld,"t#bar{t}->(ee or #mu#mu)","l");
1577     leg5->AddEntry(ttbarjzbemd,"t#bar{t}->e#mu","l");
1578     leg5->Draw();
1579     CompleteSave(can,"lepton_comparison/ttbar_emu_mc");
1580    
1581     }
1582    
1583     bool is_OF(TCut cut) {
1584     string scut = (const char*) cut;
1585     if((int)scut.find("id1!=id2")>-1) return true;
1586     if((int)scut.find("id1==id2")>-1) return false;
1587     return false;
1588     }
1589    
1590     bool is_ZP(TCut cut) {
1591     string scut = (const char*) cut;
1592     if((int)scut.find("91")>-1) return true;
1593     return false;
1594     }
1595    
1596    
1597     void draw_pure_jzb_histo(TCut cut,string datavariable, string mcvariable, string savename, TCanvas *can,vector<float> binning) {
1598     TPad *jzbpad = new TPad("jzbpad","jzbpad",0,0,1,1);
1599     jzbpad->cd();
1600     jzbpad->SetLogy(1);
1601     string xlabel="JZB [GeV]";
1602    
1603     TH1F *datahisto = allsamples.Draw("datahisto",datavariable,binning, xlabel, "events",cut,data,luminosity);
1604     THStack mcstack = allsamples.DrawStack("mcstack",mcvariable,binning, xlabel, "events",cut,mc,luminosity);
1605    
1606     datahisto->SetMinimum(0.1);
1607     datahisto->SetMarkerSize(DataMarkerSize);
1608     datahisto->Draw("e1");
1609     mcstack.Draw("same");
1610     datahisto->Draw("same,e1");
1611    
1612     TLegend *leg;
1613     if(is_OF(cut) && is_ZP(cut)) leg = allsamples.allbglegend("OFZP",datahisto);
1614     else if(!is_OF(cut) && is_ZP(cut)) leg = allsamples.allbglegend("SFZP",datahisto);
1615     else if( is_OF(cut) && !is_ZP(cut)) leg = allsamples.allbglegend("OFSB",datahisto);
1616     else if(!is_OF(cut) && !is_ZP(cut)) leg = allsamples.allbglegend("SFSB",datahisto);
1617     else {
1618     std::cerr << "Unable to decode cut: " << cut.GetTitle() << std::endl;
1619     exit(-1);
1620     }
1621     leg->Draw();
1622     string write_cut = decipher_cut(cut,"");
1623     TText *writeline1 = write_cut_on_canvas(write_cut.c_str());
1624     writeline1->SetTextSize(0.035);
1625     writeline1->Draw();
1626 buchmann 1.12 if(!Contains(savename,"Dibosons")) save_with_ratio(datahisto,mcstack,jzbpad->cd(),("jzb/"+savename));
1627     else save_with_ratio(datahisto,mcstack,jzbpad->cd(),savename);
1628 buchmann 1.1 TPad *jzbpad2 = new TPad("jzbpad2","jzbpad2",0,0,1,1);
1629     jzbpad2->cd();
1630     jzbpad2->SetLogy(1);
1631     datahisto->GetXaxis()->SetRangeUser(0,binning[binning.size()-1]);
1632     datahisto->SetMinimum(0.1);
1633     datahisto->SetMarkerSize(DataMarkerSize);
1634     datahisto->Draw("e1");
1635     mcstack.Draw("same");
1636     datahisto->Draw("same,e1");
1637     leg->SetHeader("");
1638     leg->Draw();
1639     writeline1->SetTextSize(0.035);
1640     writeline1->Draw();
1641     DrawPrelim();
1642 buchmann 1.12 if(!Contains(savename,"Dibosons")) save_with_ratio(datahisto,mcstack,jzbpad2->cd(),("jzb/PositiveSideOnly/"+savename+""));
1643     else save_with_ratio(datahisto,mcstack,jzbpad2->cd(),(savename+"__PosOnly"));
1644 buchmann 1.1 datahisto->Delete();
1645     mcstack.Delete();
1646     }
1647    
1648     Double_t GausR(Double_t *x, Double_t *par) {
1649     return gRandom->Gaus(x[0],par[0]);
1650     }
1651 buchmann 1.12
1652     void produce_stretched_jzb_plots(string mcjzb, string datajzb,vector<float> ratio_binning) {
1653     TCanvas *dican = new TCanvas("dican","JZB Plots Canvas");
1654     float max=jzbHigh ;
1655     float min=-120;
1656     int nbins=(int)((max-min)/5.0); // we want 5 GeV/bin
1657     int coarserbins=int(nbins/2.0);
1658     int rebinnedbins=int(nbins/4.0);
1659 buchmann 1.1
1660 buchmann 1.12 vector<float>binning;vector<float>coarse_binning;vector<float>coarsest_binning;
1661     for(int i=0;i<=nbins;i++)binning.push_back(min+i*(max-min)/((float)nbins));
1662     for(int i=0;i<=coarserbins;i++)coarse_binning.push_back(min+i*(max-min)/((float)coarserbins));
1663     for(int i=0;i<=rebinnedbins;i++)coarsest_binning.push_back(min+i*(max-min)/((float)rebinnedbins));
1664    
1665     draw_pure_jzb_histo(cutOSSF&&cutnJets&&cutmass,datajzb,mcjzb,"Dibosons/jzb_OS_SFZP",dican,binning);
1666     draw_pure_jzb_histo(cutOSOF&&cutnJets&&cutmass,datajzb,mcjzb,"Dibosons/jzb_OS_OFZP",dican,binning);
1667     draw_pure_jzb_histo(cutOSSF&&cutnJets&&cutmass&&"id1==0",datajzb,mcjzb,"Dibosons/ee/jzb_OS_SFZP",dican,binning);
1668     draw_pure_jzb_histo(cutOSSF&&cutnJets&&cutmass&&"id1==1",datajzb,mcjzb,"Dibosons/mm/jzb_OS_SFZP",dican,binning);
1669     draw_pure_jzb_histo(cutOSOF&&cutnJets&&cutmass&&"id1==0",datajzb,mcjzb,"Dibosons/ee/jzb_OS_OFZP",dican,binning);
1670     draw_pure_jzb_histo(cutOSOF&&cutnJets&&cutmass&&"id1==1",datajzb,mcjzb,"Dibosons/mm/jzb_OS_OFZP",dican,binning);
1671     if(PlottingSetup::RestrictToMassPeak) draw_pure_jzb_histo(cutOSSF&&cutnJets&&sidebandcut,datajzb,mcjzb,"Dibosons/jzb_OS_SFSB",dican,binning);
1672     if(PlottingSetup::RestrictToMassPeak) draw_pure_jzb_histo(cutOSOF&&cutnJets&&sidebandcut,datajzb,mcjzb,"Dibosons/jzb_OS_OFSB",dican,binning);
1673    
1674     draw_pure_jzb_histo(cutOSSF&&cutnJets&&cutmass,datajzb,mcjzb,"Dibosons/jzb_OS_SFZP_coarse",dican,coarse_binning);
1675     draw_pure_jzb_histo(cutOSOF&&cutnJets&&cutmass,datajzb,mcjzb,"Dibosons/jzb_OS_OFZP_coarse",dican,coarse_binning);
1676     if(PlottingSetup::RestrictToMassPeak) draw_pure_jzb_histo(cutOSSF&&cutnJets&&sidebandcut,datajzb,mcjzb,"Dibosons/jzb_OS_SFSB_coarse",dican,coarse_binning);
1677     if(PlottingSetup::RestrictToMassPeak) draw_pure_jzb_histo(cutOSOF&&cutnJets&&sidebandcut,datajzb,mcjzb,"Dibosons/jzb_OS_OFSB_coarse",dican,coarse_binning);
1678    
1679     delete dican;
1680     }
1681    
1682    
1683     void diboson_plots(string mcjzb, string datajzb,vector<float> ratio_binning) {
1684     vector<int> SamplesToBeModified = allsamples.FindSampleBySampleName("WW/WZ/ZZ");
1685    
1686     if(SamplesToBeModified.size()==0 || SamplesToBeModified[0]==-1) {
1687     write_error(__FUNCTION__,"Could not find any diboson samples - aborting diboson plots");
1688     return;
1689     }
1690    
1691     float stretchfactor = 100.0;
1692     vector<string> labels;
1693    
1694    
1695     dout << "Going to increase the cross section for diboson samples ... " << endl;
1696 buchmann 1.16 for(int i=0;i<(int)SamplesToBeModified.size();i++) {
1697 buchmann 1.12 float origxs=(allsamples.collection)[SamplesToBeModified[i]].xs;
1698     (allsamples.collection)[SamplesToBeModified[i]].xs=origxs*stretchfactor;
1699     dout << " Increased xs for sample " << (allsamples.collection)[SamplesToBeModified[i]].filename << " from " << origxs << " to " << (allsamples.collection)[SamplesToBeModified[i]].xs << " (by a factor of " << stretchfactor << ")" << endl;
1700     labels.push_back((allsamples.collection)[SamplesToBeModified[i]].samplename);
1701     (allsamples.collection)[SamplesToBeModified[i]].samplename=any2string(int(stretchfactor))+" x "+(allsamples.collection)[SamplesToBeModified[i]].samplename;
1702     dout << " (also renamed it to " << (allsamples.collection)[SamplesToBeModified[i]].samplename << " )" << endl;
1703     }
1704    
1705     dout << "Going to produce JZB plots" << endl;
1706 buchmann 1.13 produce_stretched_jzb_plots(mcjzb,datajzb,ratio_binning);
1707 buchmann 1.12 TCanvas *gloca = new TCanvas("gloca","gloca");
1708     float sigma=123456;
1709    
1710     dout << "Going to produce prediction plots" << endl;
1711     do_prediction_plot(mcjzb, gloca, sigma, PlottingSetup::jzbHigh, 0, false,"Dibosons/Bpred/" ); // do only MC plots, no signal
1712     do_prediction_plot(mcjzb, gloca, sigma, PlottingSetup::jzbHigh, 0, false,"Dibosons/Bpred/" ); // do MC plots with signal
1713     delete gloca;
1714    
1715     dout << "Going to reset the cross section for diboson samples ... " << endl;
1716 buchmann 1.16 for(int i=0;i<(int)SamplesToBeModified.size();i++) {
1717 buchmann 1.12 float Upxs=(allsamples.collection)[SamplesToBeModified[i]].xs;
1718     (allsamples.collection)[SamplesToBeModified[i]].xs=(allsamples.collection)[SamplesToBeModified[i]].xs*(1.0/stretchfactor);
1719     string Upname=(allsamples.collection)[SamplesToBeModified[i]].samplename;
1720     (allsamples.collection)[SamplesToBeModified[i]].samplename=labels[i];
1721     dout << " Reset xs for sample " << (allsamples.collection)[SamplesToBeModified[i]].samplename << " from " << Upxs << " to " << (allsamples.collection)[SamplesToBeModified[i]].xs << " (by a factor of " << stretchfactor << ") and reset the correct name (from " << Upname << ")" << endl;
1722    
1723     }
1724    
1725     }
1726 buchmann 1.1 void jzb_plots(string mcjzb, string datajzb,vector<float> ratio_binning) {
1727     TCanvas *can = new TCanvas("can","JZB Plots Canvas");
1728     float max=jzbHigh ;
1729     float min=-120;
1730     int nbins=(int)((max-min)/5.0); // we want 5 GeV/bin
1731     int coarserbins=int(nbins/2.0);
1732     int rebinnedbins=int(nbins/4.0);
1733    
1734     vector<float>binning;vector<float>coarse_binning;vector<float>coarsest_binning;
1735     for(int i=0;i<=nbins;i++)binning.push_back(min+i*(max-min)/((float)nbins));
1736     for(int i=0;i<=coarserbins;i++)coarse_binning.push_back(min+i*(max-min)/((float)coarserbins));
1737     for(int i=0;i<=rebinnedbins;i++)coarsest_binning.push_back(min+i*(max-min)/((float)rebinnedbins));
1738    
1739 buchmann 1.14 if ( !PlottingSetup::Approved ) {
1740     draw_pure_jzb_histo(cutOSSF&&cutnJets&&cutmass,datajzb,mcjzb,"jzb_OS_SFZP",can,binning);
1741     draw_pure_jzb_histo(cutOSOF&&cutnJets&&cutmass,datajzb,mcjzb,"jzb_OS_OFZP",can,binning);
1742     draw_pure_jzb_histo(cutOSSF&&cutnJets&&cutmass&&"id1==0",datajzb,mcjzb,"ee/jzb_OS_SFZP",can,binning);
1743     draw_pure_jzb_histo(cutOSSF&&cutnJets&&cutmass&&"id1==1",datajzb,mcjzb,"mm/jzb_OS_SFZP",can,binning);
1744     draw_pure_jzb_histo(cutOSOF&&cutnJets&&cutmass&&"id1==0",datajzb,mcjzb,"ee/jzb_OS_OFZP",can,binning);
1745     draw_pure_jzb_histo(cutOSOF&&cutnJets&&cutmass&&"id1==1",datajzb,mcjzb,"mm/jzb_OS_OFZP",can,binning);
1746     flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
1747     if(PlottingSetup::RestrictToMassPeak) draw_pure_jzb_histo(cutOSSF&&cutnJets&&sidebandcut,datajzb,mcjzb,"jzb_OS_SFSB",can,binning);
1748     if(PlottingSetup::RestrictToMassPeak) draw_pure_jzb_histo(cutOSOF&&cutnJets&&sidebandcut,datajzb,mcjzb,"jzb_OS_OFSB",can,binning);
1749     }
1750 buchmann 1.1
1751     draw_pure_jzb_histo(cutOSSF&&cutnJets&&cutmass,datajzb,mcjzb,"jzb_OS_SFZP_coarse",can,coarse_binning);
1752 buchmann 1.14 if ( !PlottingSetup::Approved ) {
1753     draw_pure_jzb_histo(cutOSOF&&cutnJets&&cutmass,datajzb,mcjzb,"jzb_OS_OFZP_coarse",can,coarse_binning);
1754     flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
1755     if(PlottingSetup::RestrictToMassPeak) draw_pure_jzb_histo(cutOSSF&&cutnJets&&sidebandcut,datajzb,mcjzb,"jzb_OS_SFSB_coarse",can,coarse_binning);
1756     if(PlottingSetup::RestrictToMassPeak) draw_pure_jzb_histo(cutOSOF&&cutnJets&&sidebandcut,datajzb,mcjzb,"jzb_OS_OFSB_coarse",can,coarse_binning);
1757     }
1758 buchmann 1.12 delete can;
1759 buchmann 1.1 }
1760    
1761    
1762     void calculate_all_yields(string mcdrawcommand,vector<float> jzb_cuts) {
1763     dout << "Calculating background yields in MC:" << endl;
1764     jzb_cuts.push_back(14000);
1765     TH1F *allbgs = allsamples.Draw("allbgs",jzbvariablemc,jzb_cuts, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc,luminosity);
1766     float cumulative=0;
1767     for(int i=allbgs->GetNbinsX();i>=1;i--) {
1768     cumulative+=allbgs->GetBinContent(i);
1769     dout << "Above " << allbgs->GetBinLowEdge(i) << " GeV/c : " << cumulative << endl;
1770     }
1771     }
1772    
1773    
1774     TCut give_jzb_expression(string mcjzb, float jzbcut, string posneg="pos") {
1775     stringstream res;
1776     res << "(" << mcjzb;
1777     if(posneg=="pos") res << ">";
1778     else res << "<-";
1779     res << jzbcut << ")";
1780     return TCut(res.str().c_str());
1781     }
1782    
1783     string sigdig(float number, int nsigdig=3, float min=0) {
1784     //this function tries to extract n significant digits, and if the number is below (min), it returns "<min"
1785     if(number<min) return "< "+any2string(min);
1786     stringstream sol;
1787     sol << setprecision(nsigdig) << number;
1788    
1789     return sol.str();
1790     }
1791    
1792     string jzb_tex_command(string region, string posneg) {
1793     if(posneg=="pos") posneg="POS";
1794     else posneg="NEG";
1795     stringstream texcommand;
1796     texcommand<<"\\"<<region <<"JZB"<<posneg;
1797     return texcommand.str();
1798     }
1799     // \SFZPJZBPOS
1800     // Sample & \OFZPJZBPOS & \OFZPJZBNEG & \SFZPJZBPOS & \SFZPJZBNEG & \OFSBJZBPOS & \OFSBJZBNEG & \SFSBJZBPOS & \SFSBJZBNEG \\\hline
1801    
1802     void compute_MC_yields(string mcjzb,vector<float> jzb_cuts) {
1803     dout << "Calculating background yields in MC:" << endl;
1804    
1805     TCanvas *yica = new TCanvas("yica","yield canvas");
1806    
1807     int nRegions=4;
1808     if(!PlottingSetup::RestrictToMassPeak) nRegions=2;
1809     string tsRegions[] = {"SFZP","OFZP","SFSB","OFSB"};
1810     string posneg[] = {"pos","neg"};
1811     TCut tkRegions[] = {cutOSSF&&cutnJets&&cutmass,cutOSOF&&cutnJets&&cutmass,cutOSSF&&cutnJets&&sidebandcut,cutOSOF&&cutnJets&&sidebandcut};
1812    
1813 buchmann 1.16 for(int ijzb=0;ijzb<(int)jzb_cuts.size();ijzb++) {
1814 buchmann 1.1 TCut jzbMC[] = { give_jzb_expression(mcjzb,jzb_cuts[ijzb],"pos"), give_jzb_expression(mcjzb,jzb_cuts[ijzb],"neg") };
1815     dout << "_________________________________________________________" << endl;
1816     dout << "Table for JZB> " << jzb_cuts[ijzb] << endl;
1817 buchmann 1.16 for(int isample=0;isample<(int)(allsamples.collection).size();isample++) {
1818 buchmann 1.1 if(!(allsamples.collection)[isample].is_data) dout << (allsamples.collection)[isample].samplename << " & ";
1819     else dout << "Sample & ";
1820     for(int iregion=0;iregion<nRegions;iregion++) {
1821     for(int ipos=0;ipos<2;ipos++) {
1822     if((allsamples.collection)[isample].is_data) dout << jzb_tex_command(tsRegions[iregion],posneg[ipos]) << " & ";
1823     else {
1824     vector<int> specific;specific.push_back(isample);
1825     TH1F *shisto = allsamples.Draw("shisto","mll",1,0,500,"tester","events",tkRegions[iregion]&&jzbMC[ipos],mc,luminosity,specific);
1826     dout << sigdig(shisto->Integral(),3,0.05) <<" & ";
1827     delete shisto;
1828     }
1829     }//end of ipos
1830     }//end of iregion
1831     dout << " \\\\" << endl;
1832     }//end of isample
1833     }//end of ijzb
1834     dout << " \\hline" << endl;
1835    
1836     delete yica;
1837     }
1838    
1839     void draw_ttbar_and_zjets_shape_for_one_configuration(string mcjzb, string datajzb, int leptontype=-1, int scenario=0,bool floating=false) {
1840     //Step 1: Establishing cuts
1841     stringstream jetcutstring;
1842     string writescenario="";
1843    
1844     if(scenario==0) jetcutstring << "(pfJetGoodNum>=3)&&"<<(const char*) basicqualitycut;
1845     if(scenario==1) jetcutstring << "(pfJetPt[0]>50&&pfJetPt[1]>50)&&"<<(const char*)basicqualitycut;
1846     TCut jetcut(jetcutstring.str().c_str());
1847     string leptoncut="mll>0";
1848     if(leptontype==0||leptontype==1) {
1849     if(leptontype==0) {
1850     leptoncut="id1==0";
1851     writescenario="__ee";
1852     }
1853     else {
1854     leptoncut="id1==1";
1855     writescenario="__ee";
1856     }
1857     }
1858     TCut lepcut(leptoncut.c_str());
1859    
1860     TCanvas *c5 = new TCanvas("c5","c5",1500,500);
1861     TCanvas *c6 = new TCanvas("c6","c6");
1862     c5->Divide(3,1);
1863    
1864     //STEP 2: Extract Zjets shape in data
1865     c5->cd(1);
1866     c5->cd(1)->SetLogy(1);
1867     TCut massat40("mll>40");
1868     TH1F *ossfleft = allsamples.Draw("ossfleft", "-"+datajzb,40,0,200, "JZB [GeV]", "events", massat40&&cutOSSF&&jetcut&&lepcut,data,luminosity);
1869     TH1F *osofleft = allsamples.Draw("osofleft", "-"+datajzb,40,0,200, "JZB [GeV]", "events", massat40&&cutOSOF&&jetcut&&lepcut,data,luminosity);
1870     ossfleft->SetLineColor(kRed);
1871     ossfleft->SetMarkerColor(kRed);
1872     ossfleft->Add(osofleft,-1);
1873     vector<TF1*> functions = do_cb_fit_to_plot(ossfleft,10);
1874     ossfleft->SetMarkerSize(DataMarkerSize);
1875     ossfleft->Draw();
1876     functions[0]->Draw("same");functions[1]->Draw("same");functions[2]->Draw("same");
1877     TF1 *zjetsfunc = (TF1*) functions[1]->Clone();
1878     TF1 *zjetsfuncN = (TF1*) functions[0]->Clone();
1879     TF1 *zjetsfuncP = (TF1*) functions[2]->Clone();
1880     zjetsfunc->Draw("same");zjetsfuncN->Draw("same");zjetsfuncP->Draw("same");
1881     TLegend *leg1 = new TLegend(0.6,0.6,0.89,0.80);
1882     leg1->SetFillColor(kWhite);
1883     leg1->SetLineColor(kWhite);
1884     leg1->AddEntry(ossfleft,"OSSF (sub),JZB<peak","p");
1885     leg1->AddEntry(zjetsfunc,"OSSF fit ('zjets')","l");
1886     leg1->Draw("same");
1887     TText *titleleft = write_title("Extracting Z+Jets shape");
1888     titleleft->Draw();
1889    
1890     //Step 3: Extract ttbar shape (in data or MC?)
1891     c5->cd(2);
1892     c5->cd(2)->SetLogy(1);
1893     TH1F *osof;
1894     TText *titlecenter;
1895     bool frommc=false;
1896     if(frommc) {
1897     osof = allsamples.Draw("osof",datajzb,40,-200,200, "JZB [GeV]", "events", massat40&&cutOSSF&&jetcut&&lepcut,mc,luminosity,allsamples.FindSample("TTJets"));
1898     titlecenter = write_title("Extracting ttbar shape (from ossf MC)");
1899     }
1900     else {
1901     osof = allsamples.Draw("osof",datajzb,40,-200,200, "JZB [GeV]", "events", massat40&&cutOSOF&&jetcut&&lepcut,data,luminosity);
1902     titlecenter = write_title("Extracting ttbar shape (from osof data)");
1903     }
1904     osof->SetMarkerSize(DataMarkerSize);
1905     osof->Draw();
1906     vector<TF1*> ttbarfunctions = do_cb_fit_to_plot(osof,35,true);
1907     ttbarfunctions[0]->SetLineColor(kRed); ttbarfunctions[0]->SetLineStyle(2); ttbarfunctions[0]->Draw("same");
1908     ttbarfunctions[1]->SetLineColor(kRed); ttbarfunctions[1]->Draw("same");
1909     ttbarfunctions[2]->SetLineColor(kRed); ttbarfunctions[2]->SetLineStyle(2); ttbarfunctions[2]->Draw("same");
1910    
1911     TLegend *leg2 = new TLegend(0.15,0.8,0.4,0.89);
1912     leg2->SetFillColor(kWhite);
1913     leg2->SetLineColor(kWhite);
1914     if(frommc) {
1915     leg2->AddEntry(osof,"t#bar{t} OSSF, MC","p");
1916     leg2->AddEntry(ttbarfunctions[1],"Fit to t#bar{t} OSSF,MC","l");
1917     } else {
1918     leg2->AddEntry(osof,"OSOF","p");
1919     leg2->AddEntry(ttbarfunctions[1],"Fit to OSOF","l");
1920     }
1921     leg2->Draw("same");
1922     titlecenter->Draw();
1923    
1924     //--------------------------------------------------------------------------------------------------------------------------------
1925     //STEP 4: Present it!
1926     // actually: if we wanna let it float we need to do that first :-)
1927     c5->cd(3);
1928     c5->cd(3)->SetLogy(1);
1929     TH1F *observed = allsamples.Draw("observed",datajzb,100,0,500, "JZB [GeV]", "events", massat40&&cutOSSF&&jetcut&&lepcut,data,luminosity);
1930     observed->SetMarkerSize(DataMarkerSize);
1931    
1932     TF1 *logparc = new TF1("logparc",InvCrystalBall,0,1000,5); logparc->SetLineColor(kRed);
1933     TF1 *logparcn = new TF1("logparcn",InvCrystalBallN,0,1000,5); logparcn->SetLineColor(kRed); logparcn->SetLineStyle(2);
1934     TF1 *logparcp = new TF1("logparcp",InvCrystalBallP,0,1000,5); logparcp->SetLineColor(kRed); logparcp->SetLineStyle(2);
1935    
1936     TF1 *zjetsc = new TF1("zjetsc",InvCrystalBall,0,1000,5); zjetsc->SetLineColor(kBlue);
1937     TF1 *zjetscn = new TF1("zjetscn",InvCrystalBallN,0,1000,5); zjetscn->SetLineColor(kBlue); zjetscn->SetLineStyle(2);
1938     TF1 *zjetscp = new TF1("zjetscp",InvCrystalBallP,0,1000,5); zjetscp->SetLineColor(kBlue); zjetscp->SetLineStyle(2);
1939    
1940     TF1 *ZplusJetsplusTTbar = new TF1("ZplusJetsplusTTbar", DoubleInvCrystalBall,0,1000,10); ZplusJetsplusTTbar->SetLineColor(kBlue);
1941     TF1 *ZplusJetsplusTTbarP= new TF1("ZplusJetsplusTTbarP",DoubleInvCrystalBallP,0,1000,10); ZplusJetsplusTTbarP->SetLineColor(kBlue); ZplusJetsplusTTbarP->SetLineStyle(2);
1942     TF1 *ZplusJetsplusTTbarN= new TF1("ZplusJetsplusTTbarN",DoubleInvCrystalBallN,0,1000,10); ZplusJetsplusTTbarN->SetLineColor(kBlue); ZplusJetsplusTTbarN->SetLineStyle(2);
1943    
1944     zjetsc->SetParameters(zjetsfunc->GetParameters());
1945     zjetscp->SetParameters(zjetsfunc->GetParameters());
1946     zjetscn->SetParameters(zjetsfunc->GetParameters());
1947    
1948     TH1F *observeda = allsamples.Draw("observeda",datajzb,53,80,jzbHigh, "JZB [GeV]", "events", massat40&&cutOSSF&&jetcut&&lepcut,data,luminosity);
1949     //blublu
1950     logparc->SetParameters(ttbarfunctions[1]->GetParameters());
1951     logparcn->SetParameters(ttbarfunctions[1]->GetParameters());
1952     logparcp->SetParameters(ttbarfunctions[1]->GetParameters());
1953     if(floating) {
1954     dout << "TTbar contribution assumed (before fitting) : " << logparc->GetParameter(0) << endl;
1955     logparc->SetParameters(ttbarfunctions[1]->GetParameters());
1956     for(int i=0;i<10;i++) {
1957     if(i<5) ZplusJetsplusTTbar->FixParameter(i,zjetsfunc->GetParameter(i));
1958     if(i>=5) {
1959     if (i>5) ZplusJetsplusTTbar->FixParameter(i,logparc->GetParameter(i-5));
1960     if (i==5) ZplusJetsplusTTbar->SetParameter(i,logparc->GetParameter(i-5));
1961     }
1962     }//end of setting parameters
1963     observeda->Draw("same");
1964     ZplusJetsplusTTbar->Draw("same");
1965     observeda->Fit(ZplusJetsplusTTbar);
1966     dout << "--> Quality of Z+Jets / TTbar fit : chi2/ndf = " << ZplusJetsplusTTbar->GetChisquare() << "/" << ZplusJetsplusTTbar->GetNDF() << endl;
1967     ZplusJetsplusTTbar->Draw("same");
1968     ZplusJetsplusTTbarP->SetParameters(ZplusJetsplusTTbar->GetParameters());
1969     ZplusJetsplusTTbarN->SetParameters(ZplusJetsplusTTbar->GetParameters());
1970     dout << "TTbar contribution found (after fitting) : " << ZplusJetsplusTTbar->GetParameter(5) << endl;
1971     float factor = ZplusJetsplusTTbar->GetParameter(5) / logparc->GetParameter(0);
1972     dout << "FACTOR: " << factor << endl;
1973     logparc->SetParameter(0,factor*ttbarfunctions[1]->GetParameter(0));
1974     logparcn->SetParameter(0,factor*ttbarfunctions[1]->GetParameter(0));
1975     logparcp->SetParameter(0,factor*ttbarfunctions[1]->GetParameter(0));
1976     }
1977    
1978     c5->cd(3);
1979     c5->cd(3)->SetLogy(1);
1980     observed->Draw();
1981     zjetsc->Draw("same");zjetscn->Draw("same");zjetscp->Draw("same");
1982     logparc->Draw("same");
1983     logparcn->Draw("same");
1984     logparcp->Draw("same");
1985    
1986     TLegend *leg3 = new TLegend(0.6,0.6,0.89,0.80);
1987     leg3->SetFillColor(kWhite);
1988     leg3->SetLineColor(kWhite);
1989     leg3->AddEntry(observed,"OSSF,JZB>peak","p");
1990     leg3->AddEntry(ttbarfunctions[1],"OSOF fit ('ttbar')","l");
1991     leg3->AddEntry(zjetsfunc,"OSSF,JZB<0 fit ('zjets')","l");
1992     leg3->Draw("same");
1993     TText *titleright = write_title("Summary of shapes and observed shape");
1994     titleright->Draw();
1995    
1996     c6->cd()->SetLogy(1);
1997     observed->Draw();
1998     zjetsc->Draw("same");zjetscn->Draw("same");zjetscp->Draw("same");
1999     logparc->Draw("same");
2000     logparcn->Draw("same");
2001     logparcp->Draw("same");
2002     leg3->Draw("same");
2003     titleright->Draw();
2004    
2005     if(scenario==0) {
2006     CompleteSave(c5,"Shapes2/Making_of___3jetsabove30"+writescenario);
2007     CompleteSave(c5->cd(1),"Shapes2/Making_of___3jetsabove30"+writescenario+"__cd1");
2008     CompleteSave(c5->cd(2),"Shapes2/Making_of___3jetsabove30"+writescenario+"__cd2");
2009     CompleteSave(c5->cd(3),"Shapes2/Making_of___3jetsabove30"+writescenario+"__cd3");
2010     CompleteSave(c6,"Shapes2/Background_Shapes___3jetsabove30"+writescenario);
2011     } else {
2012     CompleteSave(c5,"Shapes2/Making_of___2jetsabove50"+writescenario);
2013     CompleteSave(c5->cd(1),"Shapes2/Making_of___2jetsabove50"+writescenario+"__cd1");
2014     CompleteSave(c5->cd(2),"Shapes2/Making_of___2jetsabove50"+writescenario+"__cd2");
2015     CompleteSave(c5->cd(3),"Shapes2/Making_of___2jetsabove50"+writescenario+"__cd3");
2016     CompleteSave(c6,"Shapes2/Background_Shapes___2jetsabove50"+writescenario);
2017     }
2018     dout << "Statistics about our fits: " << endl;
2019     dout << "Z+Jets shape: Chi2/ndf = " << zjetsfunc->GetChisquare() << "/" << ossfleft->GetNbinsX() << endl;
2020     dout << "ttbar shape: Chi2/ndf = " << ttbarfunctions[1]->GetChisquare() << "/" << osof->GetNbinsX() << endl;
2021    
2022     c6->cd();
2023     TLegend *additionallegend = new TLegend(0.6,0.6,0.89,0.89);
2024     additionallegend->SetFillColor(kWhite);
2025     additionallegend->SetLineColor(kWhite);
2026     additionallegend->AddEntry(observed,"Data","p");
2027     additionallegend->AddEntry(ZplusJetsplusTTbar,"Fitted Z+jets & TTbar","l");
2028     additionallegend->AddEntry(zjetsc,"Z+jets","l");
2029     additionallegend->AddEntry(logparc,"TTbar","l");
2030     observed->Draw();
2031     ZplusJetsplusTTbar->SetLineColor(kGreen);
2032     ZplusJetsplusTTbarP->SetLineColor(kGreen);
2033     ZplusJetsplusTTbarN->SetLineColor(kGreen);
2034     ZplusJetsplusTTbarP->SetLineStyle(2);
2035     ZplusJetsplusTTbarN->SetLineStyle(2);
2036     TF1 *ZplusJetsplusTTbar2 = new TF1("ZplusJetsplusTTbar2",DoubleInvCrystalBall,0,1000,10);
2037     ZplusJetsplusTTbar2->SetParameters(ZplusJetsplusTTbar->GetParameters());
2038     ZplusJetsplusTTbar2->SetLineColor(kGreen);
2039     ZplusJetsplusTTbarP->SetFillColor(TColor::GetColor("#81F781"));
2040     ZplusJetsplusTTbarN->SetFillColor(kWhite);
2041     ZplusJetsplusTTbarP->Draw("fcsame");
2042     ZplusJetsplusTTbarN->Draw("fcsame");
2043     TH1F *hZplusJetsplusTTbar = (TH1F*)ZplusJetsplusTTbar2->GetHistogram();
2044     TH1F *hZplusJetsplusTTbarN = (TH1F*)ZplusJetsplusTTbarN->GetHistogram();
2045     TH1F *hZplusJetsplusTTbarP = (TH1F*)ZplusJetsplusTTbarP->GetHistogram();
2046     hZplusJetsplusTTbar->SetMarkerSize(0);
2047     hZplusJetsplusTTbarP->SetMarkerSize(0);
2048     hZplusJetsplusTTbarN->SetMarkerSize(0);
2049     for (int i=1;i<=hZplusJetsplusTTbar->GetNbinsX();i++) {
2050     float newerror=hZplusJetsplusTTbarP->GetBinContent(i)-hZplusJetsplusTTbar->GetBinContent(i);
2051     hZplusJetsplusTTbar->SetBinError(i,newerror);
2052     if(hZplusJetsplusTTbar->GetBinContent(i)<0.05) hZplusJetsplusTTbar->SetBinContent(i,0); //avoiding a displaying probolem
2053     }
2054     hZplusJetsplusTTbarP->SetFillColor(kGreen);
2055     hZplusJetsplusTTbarN->SetFillColor(kWhite);
2056     hZplusJetsplusTTbarN->Draw("same");
2057    
2058     ZplusJetsplusTTbar2->SetMarkerSize(0);
2059     ZplusJetsplusTTbar2->Draw("same");
2060    
2061     zjetsc->Draw("same");zjetscn->Draw("same");zjetscp->Draw("same");
2062     logparc->Draw("same");
2063     logparcn->Draw("same");
2064     logparcp->Draw("same");
2065     additionallegend->Draw("same");
2066     if(scenario==0) {
2067     CompleteSave(c6,"Shapes2/Background_Shapes___3jetsabove30__allfits__"+writescenario);
2068     } else {
2069     CompleteSave(c6,"Shapes2/Background_Shapes___2jetsabove50__allfits__"+writescenario);
2070     }
2071     //--------------------------------------------------------------------------------------------------------------------------------
2072     }
2073    
2074     void draw_ttbar_and_zjets_shape(string mcjzb, string datajzb) {
2075     int all_leptons=-1;
2076 buchmann 1.16 int threejetswith30gev=0;
2077     /*
2078     int twojetswith50gev=1;
2079 buchmann 1.1 int electrons_only=0;
2080     int mu_only=1;
2081 buchmann 1.16
2082 buchmann 1.1 draw_ttbar_and_zjets_shape_for_one_configuration(mcjzb,datajzb,all_leptons,twojetswith50gev);
2083     draw_ttbar_and_zjets_shape_for_one_configuration(mcjzb,datajzb,all_leptons,threejetswith30gev);
2084    
2085     draw_ttbar_and_zjets_shape_for_one_configuration(mcjzb,datajzb,electrons_only,twojetswith50gev);
2086     draw_ttbar_and_zjets_shape_for_one_configuration(mcjzb,datajzb,electrons_only,threejetswith30gev);
2087    
2088     draw_ttbar_and_zjets_shape_for_one_configuration(mcjzb,datajzb,mu_only,twojetswith50gev);
2089     draw_ttbar_and_zjets_shape_for_one_configuration(mcjzb,datajzb,mu_only,threejetswith30gev);
2090     */
2091    
2092     draw_ttbar_and_zjets_shape_for_one_configuration(mcjzb,datajzb,all_leptons,threejetswith30gev,true);
2093     }
2094    
2095     void find_correction_factors(string &jzbvardata,string &jzbvarmc) {
2096     //first: colorful plots
2097     TCanvas *cancorr = new TCanvas("cancorr","Canvas for Response Correction");
2098     cancorr->SetLogz();
2099     cancorr->SetRightMargin(0.13);
2100 buchmann 1.11 flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
2101 buchmann 1.1 TCut zptforresponsepresentation("pt<600"&&cutmass&&cutOSSF&&"((sumJetPt[1]/pt)<5.0)");
2102 buchmann 1.18 if(PlottingSetup::DoBTag) zptforresponsepresentation=zptforresponsepresentation&&PlottingSetup::bTagRequirement;
2103 buchmann 1.1 TH2F *niceresponseplotd = new TH2F("niceresponseplotd","",100,0,600,100,0,5);
2104     (allsamples.collection)[allsamples.FindSample("Data")[0]].events->Draw("sumJetPt[1]/pt:pt>>niceresponseplotd",zptforresponsepresentation);
2105     niceresponseplotd->SetStats(0);
2106     niceresponseplotd->GetXaxis()->SetTitle("Z p_{T} [GeV]");
2107     niceresponseplotd->GetYaxis()->SetTitle("Response");
2108     niceresponseplotd->GetXaxis()->CenterTitle();
2109     niceresponseplotd->GetYaxis()->CenterTitle();
2110     niceresponseplotd->Draw("COLZ");
2111     TProfile * profd = (TProfile*)niceresponseplotd->ProfileX();
2112     profd->SetMarkerSize(DataMarkerSize);
2113 buchmann 1.18 profd->Fit("pol0","","same,e1",100,400);
2114 buchmann 1.1 DrawPrelim();
2115     TText* title = write_text(0.5,0.7,"Data");
2116     title->SetTextAlign(12);
2117     title->Draw();
2118     TF1 *datapol=(TF1*)profd->GetFunction("pol0");
2119     float datacorrection=datapol->GetParameter(0);
2120     stringstream dataresstring;
2121     dataresstring<<"Response: "<<std::setprecision(2)<<100*datacorrection<<" %";
2122     TText* restitle = write_text(0.5,0.65,dataresstring.str());
2123     restitle->SetTextAlign(12);
2124     restitle->SetTextSize(0.03);
2125     restitle->Draw();
2126     CompleteSave(cancorr,"ResponseCorrection/Response_Correction_Illustration_Data");
2127    
2128     TH2F *niceresponseplotm = new TH2F("niceresponseplotm","",100,0,600,100,0,5);
2129 buchmann 1.22 (allsamples.collection)[allsamples.FindSample("DY")[0]].events->Draw("sumJetPt[1]/pt:pt>>niceresponseplotm",zptforresponsepresentation*cutWeight);
2130 buchmann 1.1 niceresponseplotm->SetStats(0);
2131     niceresponseplotm->GetXaxis()->SetTitle("Z p_{T} [GeV]");
2132     niceresponseplotm->GetYaxis()->SetTitle("Response");
2133     niceresponseplotm->GetXaxis()->CenterTitle();
2134     niceresponseplotm->GetYaxis()->CenterTitle();
2135     niceresponseplotm->Draw("COLZ");
2136     (allsamples.collection)[allsamples.FindSample("DY")[0]].events->Draw("sumJetPt[1]/pt:pt",zptforresponsepresentation,"PROF,same");
2137     TProfile * profm = (TProfile*)niceresponseplotm->ProfileX();
2138     profm->SetMarkerSize(DataMarkerSize);
2139 buchmann 1.18 profm->Fit("pol0","","same,e1",100,400);
2140 buchmann 1.1 DrawMCPrelim();
2141     title = write_text(0.5,0.7,"MC simulation");
2142     title->SetTextAlign(12);
2143     title->Draw();
2144     TF1 *mcpol=(TF1*)profm->GetFunction("pol0");
2145     float mccorrection=mcpol->GetParameter(0);
2146     stringstream mcresstring;
2147     mcresstring<<"Response: "<<std::setprecision(2)<<100*mccorrection<<" %";
2148     TText* mcrestitle = write_text(0.5,0.65,mcresstring.str());
2149     mcrestitle->SetTextAlign(12);
2150     mcrestitle->SetTextSize(0.03);
2151     mcrestitle->Draw();
2152     CompleteSave(cancorr,"ResponseCorrection/Response_Correction_Illustration_MC");
2153    
2154    
2155     //Step 2: Getting the result
2156     // TCut zptcutforresponse("pt>30&&pt<300&&TMath::Abs(91.2-mll)<20&&id1==id2&&(ch1*ch2<0)");
2157     stringstream jzbvardatas;
2158     if(datacorrection>1) jzbvardatas<<"(jzb[1]-"<<datacorrection-1<<"*pt)";
2159     if(datacorrection<1) jzbvardatas<<"(jzb[1]+"<<1-datacorrection<<"*pt)";
2160     jzbvardata=jzbvardatas.str();
2161     stringstream jzbvarmcs;
2162     if(mccorrection>1) jzbvarmcs<<"(jzb[1]-"<<mccorrection-1<<"*pt)";
2163     if(mccorrection<1) jzbvarmcs<<"(jzb[1]+"<<1-mccorrection<<"*pt)";
2164     jzbvarmc=jzbvarmcs.str();
2165     dout << "JZB Z pt correction summary : " << endl;
2166     dout << " Data: The response is " << datacorrection << " --> jzb variable is now : " << jzbvardata << endl;
2167     dout << " MC : The response is " << mccorrection << " --> jzb variable is now : " << jzbvarmc << endl;
2168     }
2169    
2170     void pick_up_events(string cut) {
2171     dout << "Picking up events with cut " << cut << endl;
2172     allsamples.PickUpEvents(cut);
2173     }
2174    
2175 buchmann 1.18 void save_template(string mcjzb, string datajzb,vector<float> jzb_cuts,float MCPeakError,float DataPeakError, vector<float> jzb_shape_limit_bins) {
2176 buchmann 1.1 dout << "Saving configuration template!" << endl;
2177     ofstream configfile;
2178     configfile.open("../DistributedModelCalculations/last_configuration.C");
2179     configfile<<"#include <iostream>\n";
2180     configfile<<"#include <vector>\n";
2181     configfile<<"#ifndef SampleClassLoaded\n";
2182     configfile<<"#include \"SampleClass.C\"\n";
2183     configfile<<"#endif\n";
2184     configfile<<"#define SetupLoaded\n";
2185     configfile<<"#ifndef ResultLibraryClassLoaded\n";
2186     configfile<<"#include \"ResultLibraryClass.C\"\n";
2187     configfile<<"#endif\n";
2188    
2189     configfile<<"\nusing namespace std;\n\n";
2190    
2191     configfile<<"namespace PlottingSetup { \n";
2192     configfile<<"string datajzb=\"datajzb_ERROR\";\n";
2193     configfile<<"string mcjzb=\"mcjzb_ERROR\";\n";
2194     configfile<<"vector<float>jzb_cuts;\n";
2195 buchmann 1.18 configfile<<"vector<float>jzb_shape_limit_bins;\n";
2196 buchmann 1.1 configfile<<"float MCPeakError=-999;\n";
2197 buchmann 1.11 configfile<<"float DataPeakError=-999;\n";
2198 buchmann 1.1 configfile<<"}\n\n";
2199    
2200     configfile<<"void read_config() {\n";
2201     configfile<<"datajzb=\""<<datajzb<<"\";\n";
2202     configfile<<"mcjzb=\""<<mcjzb<<"\";\n\n";
2203 buchmann 1.11 configfile<<"\n\nMCPeakError="<<MCPeakError<<";\n";
2204     configfile<<"DataPeakError="<<DataPeakError<<";\n\n";
2205 buchmann 1.16 for(int i=0;i<(int)jzb_cuts.size();i++) configfile<<"jzb_cuts.push_back("<<jzb_cuts[i]<<"); // JZB cut at " << jzb_cuts[i] << "\n";
2206 buchmann 1.1 configfile<<"\n\n";
2207 buchmann 1.16 for(int i=0;i<(int)Nobs.size();i++) configfile<<"Nobs.push_back("<<Nobs[i]<<"); // JZB cut at " << jzb_cuts[i] << "\n";
2208     for(int i=0;i<(int)Npred.size();i++) configfile<<"Npred.push_back("<<Npred[i]<<"); // JZB cut at " << jzb_cuts[i] << "\n";
2209     for(int i=0;i<(int)Nprederr.size();i++) configfile<<"Nprederr.push_back("<<Nprederr[i]<<"); // JZB cut at " << jzb_cuts[i] << "\n";
2210 buchmann 1.1 configfile<<"\n\n";
2211 buchmann 1.16 for(int i=0;i<(int)flippedNobs.size();i++) configfile<<"flippedNobs.push_back("<<flippedNobs[i]<<"); // JZB cut at " << jzb_cuts[i] << "\n";
2212     for(int i=0;i<(int)flippedNpred.size();i++) configfile<<"flippedNpred.push_back("<<flippedNpred[i]<<"); // JZB cut at " << jzb_cuts[i] << "\n";
2213 buchmann 1.21 for(int i=0;i<(int)flippedNprederr.size();i++) configfile<<"flippedNprederr.push_back("<<flippedNprederr[i]<<"); // JZB cut at " << jzb_cuts[i] << "\n";
2214 buchmann 1.18 for(int i=0;i<(int)jzb_shape_limit_bins.size();i++) configfile<<"jzb_shape_limit_bins.push_back("<<jzb_shape_limit_bins[i]<<"); // JZB shape bin boundary at " << jzb_shape_limit_bins[i] << "\n";
2215     configfile<<"\n\n";
2216 buchmann 1.1 configfile<<"\n\n";
2217     configfile<<"luminosity="<<luminosity<<";\n";
2218 buchmann 1.5 configfile<<"RestrictToMassPeak="<<RestrictToMassPeak<<";//defines the type of analysis we're running\n";
2219 buchmann 1.1
2220     configfile<<"\n\ncout << \"Configuration successfully loaded!\" << endl; \n \n } \n \n";
2221    
2222     configfile.close();
2223    
2224     }
2225    
2226     float get_nonzero_minimum(TH1F *histo) {
2227     float min=histo->GetMaximum();
2228     for(int ibin=1;ibin<=histo->GetNbinsX();ibin++) {
2229     float curcont=histo->GetBinContent(ibin);
2230     if(curcont<min&&curcont>0) min=curcont;
2231     }
2232     return min;
2233     }
2234    
2235     void draw_all_ttbar_histos(TCanvas *can, vector<TH1F*> histos, string drawoption="", float manualmin=-9) {
2236     can->cd();
2237     float min=1;
2238     float max=histos[0]->GetMaximum();
2239     if(manualmin>=0) min=manualmin;
2240     else {
2241 buchmann 1.16 for(int i=1;i<(int)histos.size();i++) {
2242 buchmann 1.1 float curmin=get_nonzero_minimum(histos[i]);
2243     float curmax=histos[i]->GetMaximum();
2244     if(curmin<min) min=curmin;
2245     if(curmax>max) max=curmax;
2246     }
2247     }
2248     histos[0]->GetYaxis()->SetRangeUser(min,4*max);
2249     histos[0]->Draw(drawoption.c_str());
2250     stringstream drawopt;
2251     drawopt << drawoption << ",same";
2252 buchmann 1.16 for(int i=1;i<(int)histos.size();i++) {
2253 buchmann 1.1 histos[i]->Draw(drawopt.str().c_str());
2254     }
2255     }
2256    
2257     void ttbar_sidebands_comparison(string mcjzb, vector<float> binning, string prestring) {
2258     //in the case of the on peak analysis, we compare the 3 control regions to the real value
2259     //in the case of the OFF peak analysis, we compare our control region to the real value
2260     TCut weightbackup=cutWeight;
2261 buchmann 1.21 // cutWeight="1.0";
2262 buchmann 1.1 float simulatedlumi = luminosity; //in pb please - adjust to your likings
2263    
2264    
2265     TH1F *TZem = systsamples.Draw("TZem",mcjzb,binning,"JZB [GeV]","events",cutmass&&cutOSOF&&cutnJets,mc,simulatedlumi,systsamples.FindSample("TTJets"));
2266     TH1F *nTZem = systsamples.Draw("nTZem","-"+mcjzb,binning,"JZB [GeV]","events",cutmass&&cutOSOF&&cutnJets,mc,simulatedlumi,systsamples.FindSample("TTJets"));
2267     TH1F *TSem;
2268     TH1F *nTSem;
2269     TH1F *TZeemm = systsamples.Draw("TZeemm",mcjzb,binning,"JZB [GeV]","events",cutmass&&cutOSSF&&cutnJets,mc,simulatedlumi,systsamples.FindSample("TTJets"));
2270     TH1F *nTZeemm = systsamples.Draw("nTZeemm","-"+mcjzb,binning,"JZB [GeV]","events",cutmass&&cutOSSF&&cutnJets,mc,simulatedlumi,systsamples.FindSample("TTJets"));
2271     TH1F *TSeemm;
2272     TH1F *nTSeemm;
2273    
2274     if(PlottingSetup::RestrictToMassPeak) {
2275     TSem = systsamples.Draw("TSem",mcjzb,binning,"JZB [GeV]","events",sidebandcut&&cutOSOF&&cutnJets,mc,simulatedlumi,systsamples.FindSample("TTJets"));
2276     nTSem = systsamples.Draw("nTSem","-"+mcjzb,binning,"JZB [GeV]","events",sidebandcut&&cutOSOF&&cutnJets,mc,simulatedlumi,systsamples.FindSample("TTJets"));
2277     TSeemm = systsamples.Draw("TSeemm",mcjzb,binning,"JZB [GeV]","events",sidebandcut&&cutOSSF&&cutnJets,mc,simulatedlumi,systsamples.FindSample("TTJets"));
2278     nTSeemm = systsamples.Draw("nTSeemm","-"+mcjzb,binning,"JZB [GeV]","events",sidebandcut&&cutOSSF&&cutnJets,mc,simulatedlumi,systsamples.FindSample("TTJets"));
2279     }
2280    
2281     TCanvas *tcan = new TCanvas("tcan","tcan");
2282    
2283     tcan->SetLogy(1);
2284    
2285     TZeemm->SetLineColor(kBlack);
2286     TZem->SetLineColor(kRed);
2287     TZeemm->SetMarkerColor(kBlack);
2288     TZem->SetMarkerColor(kRed);
2289    
2290    
2291     if(PlottingSetup::RestrictToMassPeak) {
2292     TSem->SetLineColor(TColor::GetColor("#00A616"));
2293     TSeemm->SetLineColor(kMagenta+2);
2294     TSem->SetMarkerColor(TColor::GetColor("#00A616"));
2295     TSeemm->SetMarkerColor(kMagenta+2);
2296     TSem->SetLineStyle(2);
2297     TSeemm->SetLineStyle(3);
2298     }
2299    
2300     vector<TH1F*> histos;
2301     histos.push_back(TZem);
2302     histos.push_back(TZeemm);
2303     if(PlottingSetup::RestrictToMassPeak) {
2304     histos.push_back(TSem);
2305     histos.push_back(TSeemm);
2306     }
2307     draw_all_ttbar_histos(tcan,histos,"histo",0.04);
2308    
2309     TLegend *leg = make_legend("MC t#bar{t}",0.6,0.65,false);
2310     leg->AddEntry(TZeemm,"SFZP","l");
2311     leg->AddEntry(TZem,"OFZP","l");
2312     if(PlottingSetup::RestrictToMassPeak) {
2313     leg->AddEntry(TSeemm,"SFSB","l");
2314     leg->AddEntry(TSem,"OFSB","l");
2315     }
2316     leg->Draw("same");
2317     DrawMCPrelim(simulatedlumi);
2318     CompleteSave(tcan,"Systematics/"+prestring+"/ttbar_shape_comparison");
2319     TH1F *TZemcopy = (TH1F*)TZem->Clone("TZemcopy");
2320     TH1F *TZeemmcopy = (TH1F*)TZeemm->Clone("TZeemmcopy");
2321 buchmann 1.11 TH1F *TSeemmcopy;
2322     TH1F *TSemcopy;
2323     if(PlottingSetup::RestrictToMassPeak) {
2324     TSeemmcopy = (TH1F*)TSeemm->Clone("TSeemmcopy");
2325     TSemcopy = (TH1F*)TSem->Clone("TSemcopy");
2326     }
2327 buchmann 1.1
2328     TZem->Divide(TZeemm);
2329     TZem->SetMarkerStyle(21);
2330     if(PlottingSetup::RestrictToMassPeak) {
2331     TSem->Divide(TZeemm);
2332     TSeemm->Divide(TZeemm);
2333     TSem->SetMarkerStyle(24);
2334     TSeemm->SetMarkerStyle(32);
2335     }
2336    
2337     tcan->SetLogy(0);
2338     TZem->GetYaxis()->SetRangeUser(0,2.5);
2339     TZem->GetYaxis()->SetTitle("ratio");
2340     TZem->Draw();
2341     if(PlottingSetup::RestrictToMassPeak) {
2342     TSem->Draw("same");
2343     TSeemm->Draw("same");
2344     }
2345    
2346     float linepos=0.25;
2347     if(PlottingSetup::RestrictToMassPeak) linepos=0.25;
2348 buchmann 1.11 if(!PlottingSetup::RestrictToMassPeak) linepos=0.1; //sys uncertainty for iJZB
2349 buchmann 1.1 TLine *top = new TLine(binning[0],1.0+linepos,binning[binning.size()-1],1.0+linepos);
2350     TLine *center = new TLine(binning[0],1.0,binning[binning.size()-1],1.0);
2351     TLine *bottom = new TLine(binning[0],1.0-linepos,binning[binning.size()-1],1.0-linepos);
2352    
2353 buchmann 1.11 /*write_warning(__FUNCTION__,"These two lines are to be removed!");
2354     TLine *topalt = new TLine(binning[0],1.0+0.1,binning[binning.size()-1],1.0+0.1);
2355     TLine *bottomalt = new TLine(binning[0],1.0-0.1,binning[binning.size()-1],1.0-0.1);
2356     topalt->SetLineColor(kRed);topalt->SetLineStyle(3);
2357     bottomalt->SetLineColor(kRed);bottomalt->SetLineStyle(3);
2358     topalt->Draw("same");bottomalt->Draw("same");*/
2359    
2360    
2361 buchmann 1.1 top->SetLineColor(kBlue);top->SetLineStyle(2);
2362     bottom->SetLineColor(kBlue);bottom->SetLineStyle(2);
2363     center->SetLineColor(kBlue);
2364    
2365     top->Draw("same");
2366     center->Draw("same");
2367     bottom->Draw("same");
2368    
2369     TLegend *leg2 = make_legend("MC t#bar{t}",0.55,0.75,false);
2370     leg2->AddEntry(TZem,"OFZP / SFZP","ple");
2371     if(PlottingSetup::RestrictToMassPeak) {
2372     leg2->AddEntry(TSeemm,"SFSB / SFZP","ple");
2373     leg2->AddEntry(TSem,"OFSB / SFZP","ple");
2374     }
2375     leg2->AddEntry(bottom,"syst. envelope","l");
2376     leg2->SetX1(0.25);leg2->SetX2(0.6);
2377     leg2->SetY1(0.65);
2378    
2379     leg2->Draw("same");
2380    
2381     DrawMCPrelim(simulatedlumi);
2382     CompleteSave(tcan,"Systematics/"+prestring+"/ttbar_shape_comparison_ratio");
2383    
2384    
2385     ///-------------- second part: only look at the quantity we actually care about!
2386     TH1F *leftsfzp = (TH1F*)nTZeemm->Clone("leftsfzp");
2387     TH1F *rightsfzp = (TH1F*)TZeemmcopy->Clone("rightsfzp");
2388     rightsfzp->Add(leftsfzp,-1);
2389     TH1F *leftofzp = (TH1F*)nTZem->Clone("leftofzp");
2390     TH1F *rightofzp = (TH1F*)TZemcopy->Clone("rightofzp");
2391     rightofzp->Add(leftofzp,-1);
2392     TH1F *leftofsb;
2393     TH1F *rightofsb;
2394     TH1F *leftsfsb;
2395     TH1F *rightsfsb;
2396     if(PlottingSetup::RestrictToMassPeak) {
2397     leftofsb = (TH1F*)nTSem->Clone("leftofsb");
2398     rightofsb = (TH1F*)TSemcopy->Clone("rightofsb");
2399     rightofsb->Add(leftofsb,-1);
2400     leftsfsb = (TH1F*)nTSeemm->Clone("leftsfsb");
2401     rightsfsb = (TH1F*)TSeemmcopy->Clone("rightsfsb");
2402     rightsfsb->Add(leftsfsb,-1);
2403     }
2404    
2405     tcan->SetLogy(1);
2406     rightsfzp->GetXaxis()->SetRangeUser(0,binning[binning.size()-1]);
2407     rightsfzp->GetYaxis()->SetTitle("#deltaJZB / 25 GeV");
2408     rightsfzp->Draw("histo");
2409     rightofzp->Draw("histo,same");
2410     if(PlottingSetup::RestrictToMassPeak) {
2411     rightofsb->Draw("histo,same");
2412     rightsfsb->Draw("histo,same");
2413     }
2414     DrawMCPrelim(simulatedlumi);
2415    
2416     TLegend *legA = make_legend("MC t#bar{t}",0.6,0.65,false);
2417     legA->AddEntry(rightsfzp,"SFZP","l");
2418     legA->AddEntry(rightofzp,"OFZP","l");
2419     if(PlottingSetup::RestrictToMassPeak) {
2420     legA->AddEntry(rightofsb,"SFSB","l");
2421     legA->AddEntry(rightsfsb,"OFSB","l");
2422     }
2423     legA->Draw();
2424    
2425     CompleteSave(tcan,"Systematics/"+prestring+"/ttbar_deltajzb_comparison");
2426    
2427     tcan->SetLogy(0);
2428     rightofzp->Divide(rightsfzp);
2429     rightofzp->GetXaxis()->SetRangeUser(0.0,binning[binning.size()-1]);
2430     rightofzp->GetYaxis()->SetRangeUser(0.0,2.0);
2431     rightofzp->GetYaxis()->SetTitle("#deltaJZB ratio");
2432     rightofzp->Draw();
2433     if(PlottingSetup::RestrictToMassPeak) {
2434     rightofsb->Divide(rightsfzp);
2435     rightsfsb->Divide(rightsfzp);
2436     rightofsb->Draw("same");
2437     rightsfsb->Draw("same");
2438     }
2439    
2440     TLine *top2 = new TLine(0.0,1.0+linepos,binning[binning.size()-1],1.0+linepos);
2441     TLine *center2 = new TLine(0.0,1.0,binning[binning.size()-1],1.0);
2442     TLine *bottom2 = new TLine(0.0,1.0-linepos,binning[binning.size()-1],1.0-linepos);
2443    
2444     top2->SetLineColor(kBlue);top2->SetLineStyle(2);
2445     bottom2->SetLineColor(kBlue);bottom2->SetLineStyle(2);
2446     center2->SetLineColor(kBlue);
2447    
2448     top2->Draw("same");
2449     center2->Draw("same");
2450     bottom2->Draw("same");
2451    
2452     rightofzp->SetMarkerStyle(21);
2453     if(PlottingSetup::RestrictToMassPeak) {
2454     rightofsb->SetMarkerStyle(24);
2455     rightsfsb->SetMarkerStyle(32);
2456     }
2457    
2458     TLegend *leg3 = make_legend("MC t#bar{t}",0.55,0.75,false);
2459     leg3->AddEntry(rightofzp,"OFZP / SFZP","ple");
2460     if(PlottingSetup::RestrictToMassPeak) {
2461     leg3->AddEntry(rightsfsb,"SFSB / SFZP","ple");
2462     leg3->AddEntry(rightofsb,"OFSB / SFZP","ple");
2463     }
2464     leg3->AddEntry(bottom,"syst. envelope","l");
2465     leg3->SetX1(0.25);leg3->SetX2(0.6);
2466     leg3->SetY1(0.65);
2467    
2468     leg3->Draw("same");
2469    
2470     DrawMCPrelim(simulatedlumi);
2471     CompleteSave(tcan,"Systematics/"+prestring+"/ttbar_deltajzb_comparison_ratio");
2472    
2473     delete TZem;
2474     delete nTZem;
2475     delete TZeemm;
2476     delete nTZeemm;
2477     if(PlottingSetup::RestrictToMassPeak) {
2478     delete TSem;
2479     delete nTSem;
2480     delete TSeemm;
2481     delete nTSeemm;
2482     }
2483    
2484     delete tcan;
2485     cutWeight=weightbackup;
2486     }
2487    
2488     void ttbar_sidebands_comparison(string mcjzb, vector<float> jzb_binning) {
2489     vector<float> nicer_binning;
2490     nicer_binning.push_back(-125);
2491     nicer_binning.push_back(-100);
2492     nicer_binning.push_back(-75);
2493     nicer_binning.push_back(-50);
2494     nicer_binning.push_back(-25);
2495     nicer_binning.push_back(0);
2496     nicer_binning.push_back(25);
2497     nicer_binning.push_back(50);
2498     nicer_binning.push_back(75);
2499     nicer_binning.push_back(100);
2500     nicer_binning.push_back(125);
2501     nicer_binning.push_back(150);
2502     nicer_binning.push_back(175);
2503     nicer_binning.push_back(200);
2504     nicer_binning.push_back(230);
2505     nicer_binning.push_back(280);
2506     nicer_binning.push_back(400);
2507     ttbar_sidebands_comparison(mcjzb,nicer_binning, "ttbar/");
2508     }
2509    
2510    
2511 buchmann 1.20 void zjets_prediction_comparison(string mcjzbWithPU) {
2512     TCanvas *zcan = new TCanvas("zcan","zcan");
2513     zcan->SetLogy(1);
2514     TCut weightbackup=cutWeight;
2515    
2516     /*
2517 buchmann 1.1 // Do it without PU re-weighting
2518     float MCPeakNoPU=0,MCPeakErrorNoPU=0,DataPeakNoPU=0,DataPeakErrorNoPU=0,MCSigma=0,DataSigma=0;
2519     stringstream resultsNoPU;
2520    
2521     stringstream mcjzbnoPU;
2522 buchmann 1.20 find_peaks(MCPeakNoPU,MCPeakErrorNoPU, DataPeakNoPU, DataPeakErrorNoPU,MCSigma,DataSigma,resultsNoPU,true);
2523 buchmann 1.1 if(MCPeakNoPU>0) mcjzbnoPU<<"("<<jzbvariablemc<<"-"<<TMath::Abs(MCPeakNoPU)<<")";
2524     else mcjzbnoPU<<"("<<jzbvariablemc<<"+"<<TMath::Abs(MCPeakNoPU)<<")";
2525    
2526     string mcjzb = mcjzbnoPU.str();
2527     dout << "The peak corrected JZB expression for MC without pileup is : " << mcjzb << endl;
2528    
2529     cutWeight="1.0";
2530 buchmann 1.20 */
2531     string mcjzb = mcjzbWithPU; // this is with PURW, if you want without it you have to uncomment the part above (and comment out this line)
2532 buchmann 1.1 float sbg_min=0.;
2533     float sbg_max=100.;
2534     int sbg_nbins=5;
2535     float simulatedlumi = luminosity;//in pb please - adjust to your likings
2536    
2537     TCut kPos((mcjzb+">0").c_str());
2538     TCut kNeg((mcjzb+"<0").c_str());
2539     string var( "abs("+mcjzb+")" );
2540    
2541 buchmann 1.20 TCut kcut(cutmass&&cutOSSF&&cutnJets);
2542     TH1F *hJZBpos = systsamples.Draw("hJZBpos",var,sbg_nbins,sbg_min,sbg_max, "JZB [GeV]", "events",kcut&&kPos,mc,simulatedlumi,systsamples.FindSample("/DY"));
2543     TH1F *hJZBneg = systsamples.Draw("hJZBneg",var,sbg_nbins,sbg_min,sbg_max, "JZB [GeV]", "events",kcut&&kNeg,mc,simulatedlumi,systsamples.FindSample("/DY"));
2544 buchmann 1.1 hJZBpos->SetLineColor(kBlack);
2545     hJZBneg->SetLineColor(kRed);
2546    
2547    
2548     hJZBpos->Draw("e1");
2549     hJZBneg->Draw("same,hist");
2550     hJZBpos->Draw("same,e1"); // So it's on top...
2551    
2552     TLegend *leg = make_legend("MC Z+jets",0.55,0.75,false);
2553     leg->AddEntry(hJZBpos,"Observed","pe");
2554     leg->AddEntry(hJZBneg,"Predicted","l");
2555     leg->Draw("same");
2556     DrawMCPrelim(simulatedlumi);
2557     CompleteSave(zcan,"Systematics/zjets_prediction");
2558    
2559     TH1F* hratio = (TH1F*)hJZBpos->Clone("hratio");
2560     hratio->Divide(hJZBneg);
2561    
2562     zcan->SetLogy(0);
2563     hratio->GetYaxis()->SetRangeUser(0,2.5);
2564     hratio->GetYaxis()->SetTitle("Observed/Predicted");
2565     hratio->Draw("e1");
2566    
2567     TLine *top = new TLine(sbg_min,1.25,sbg_max,1.25);
2568     TLine *center = new TLine(sbg_min,1.0,sbg_max,1.0);
2569     TLine *bottom = new TLine(sbg_min,0.75,sbg_max,0.75);
2570    
2571    
2572     top->SetLineColor(kBlue);top->SetLineStyle(2);
2573     bottom->SetLineColor(kBlue);bottom->SetLineStyle(2);
2574     center->SetLineColor(kBlue);
2575    
2576     top->Draw("same");
2577     center->Draw("same");
2578     bottom->Draw("same");
2579    
2580     TLegend *leg2 = make_legend("MC Z+jets",0.25,0.75,false);
2581     leg2->AddEntry(hratio,"obs / pred","pe");
2582     leg2->AddEntry(bottom,"syst. envelope","l");
2583     leg2->Draw("same");
2584     DrawMCPrelim(simulatedlumi);
2585     CompleteSave(zcan,"Systematics/zjets_prediction_ratio");
2586    
2587     delete zcan;
2588     cutWeight=weightbackup;
2589    
2590     }
2591    
2592    
2593    
2594     void sideband_assessment(string datajzb, float min=30.0, float max=50.0) {
2595     tout << endl << endl;
2596     stringstream bordercut;
2597     bordercut << "(TMath::Abs(" << datajzb << ")<" << max << ")&&(TMath::Abs(" << datajzb << ")>" << min << ")";
2598     tout << bordercut.str().c_str() << endl;
2599     TH1F *ofsb = allsamples.Draw("ofsb",datajzb,100, 0,100, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets&&bordercut.str().c_str(),data, luminosity);
2600     TH1F *ofzp = allsamples.Draw("ofzp",datajzb,100, 0,100, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets&&bordercut.str().c_str(),data, luminosity);
2601     float OFSB = ofsb->Integral();
2602     float OFZP = ofzp->Integral();
2603    
2604     tout << "\\begin{table}[hbtp]" << endl;
2605     tout << "\\renewcommand{\\arraystretch}{1.3}" << endl;
2606     tout << "\\begin{center}" << endl;
2607     tout << "\\caption{Total number of events observed in the auxiliary region $|JZB|>"<<min<<"\\GeV \\cup |JZB|<"<<max<<"\\GeV$ for the {\\OFZP} and {\\OFSB}.}\\label{tab:auxcounts}" << endl;
2608     tout << "\\begin{tabular}{l|cc}" << endl;
2609     tout << "\\hline" << endl;
2610     tout << "& {\\OFZP} & {\\OFSB} \\\\\\hline" << endl;
2611     tout << "\\#(events) & "<<OFSB<<" & "<<OFZP<<"\\\\ \\hline" << endl;
2612     tout << "\\end{tabular}" << endl;
2613     tout << "\\end{center}" << endl;
2614     tout << "\\end{table}" << endl;
2615    
2616    
2617     }
2618    
2619     void make_table(samplecollection &coll, string jzbexpr, bool is_data, vector<float> jzb_cuts, string subselection="none") {
2620    
2621     vector<float> jzbcutprediction;
2622     vector<float> metcutprediction;
2623    
2624     vector<float> jzbcutobservation;
2625     vector<float> metcutobservation;
2626    
2627     TCanvas *cannie = new TCanvas("cannie","cannie");
2628    
2629 buchmann 1.16 for(int icut=0;icut<(int)jzb_cuts.size();icut++) {
2630 buchmann 1.1 float currcut=jzb_cuts[icut];
2631     int nbins=1;float low=currcut;
2632     vector<int> mysample;
2633     if(subselection!="none") mysample=coll.FindSample(subselection);
2634     TH1F *RcorrJZBeemm = coll.Draw("RcorrJZBeemm",jzbexpr,1,currcut,14000, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,is_data, luminosity, mysample);
2635     TH1F *LcorrJZBeemm = coll.Draw("LcorrJZBeemm",("-("+jzbexpr+")").c_str(),nbins,low,14000, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,is_data, luminosity, mysample);
2636     TH1F *RcorrJZBem = coll.Draw("RcorrJZBem",jzbexpr.c_str(),nbins,low,14000, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,is_data, luminosity, mysample);
2637     TH1F *LcorrJZBem = coll.Draw("LcorrJZBem",("-("+jzbexpr+")").c_str(),nbins,low,14000, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,is_data, luminosity, mysample);
2638    
2639     TH1F *RcorrJZBSBem;
2640     TH1F *LcorrJZBSBem;
2641     TH1F *RcorrJZBSBeemm;
2642     TH1F *LcorrJZBSBeemm;
2643    
2644     if(PlottingSetup::RestrictToMassPeak) {
2645     RcorrJZBSBem = coll.Draw("RcorrJZBSBem",jzbexpr.c_str(),nbins,low,14000, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets,is_data, luminosity, mysample);
2646     LcorrJZBSBem = coll.Draw("LcorrJZBSBem",("-("+jzbexpr+")").c_str(),nbins,low,14000, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets,is_data, luminosity, mysample);
2647     RcorrJZBSBeemm = coll.Draw("RcorrJZBSBeemm",jzbexpr.c_str(),nbins,low,14000, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets,is_data, luminosity, mysample);
2648     LcorrJZBSBeemm = coll.Draw("LcorrJZBSBeemm",("-("+jzbexpr+")").c_str(),nbins,low,14000, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets,is_data, luminosity, mysample);
2649     }
2650    
2651     TH1F *Bpred = (TH1F*)LcorrJZBeemm->Clone("Bpred");
2652     if(PlottingSetup::RestrictToMassPeak) {
2653     Bpred->Add(RcorrJZBem,1.0/3.);
2654     Bpred->Add(LcorrJZBem,-1.0/3.);
2655     Bpred->Add(RcorrJZBSBem,1.0/3.);
2656     Bpred->Add(LcorrJZBSBem,-1.0/3.);
2657     Bpred->Add(RcorrJZBSBeemm,1.0/3.);
2658     Bpred->Add(LcorrJZBSBeemm,-1.0/3.);
2659     } else {
2660     Bpred->Add(RcorrJZBem,1.0);
2661     Bpred->Add(LcorrJZBem,-1.0);
2662     }
2663    
2664     jzbcutobservation.push_back(RcorrJZBeemm->Integral());
2665     jzbcutprediction.push_back(Bpred->Integral());
2666    
2667     delete RcorrJZBeemm;
2668     delete LcorrJZBeemm;
2669     delete RcorrJZBem;
2670     delete LcorrJZBem;
2671     delete RcorrJZBSBem;
2672     delete LcorrJZBSBem;
2673     delete RcorrJZBSBeemm;
2674     delete LcorrJZBSBeemm;
2675    
2676     TH1F *MetObs = coll.Draw("MetObs",("met[4]"),nbins,low,14000, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,is_data, luminosity, mysample);
2677     TH1F *MetPred = coll.Draw("MetPred",("met[4]"),nbins,low,14000, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,is_data, luminosity, mysample);
2678    
2679     metcutobservation.push_back(MetObs->Integral());
2680     metcutprediction.push_back(MetPred->Integral());
2681     delete MetObs;
2682     delete MetPred;
2683     }//end of cut loop
2684    
2685     //prediction part
2686     if(is_data) cout << "Data prediction & ";
2687     if(subselection!="none") cout << subselection << " prediction &";
2688 buchmann 1.16 for(int ij=0;ij<(int)jzb_cuts.size();ij++) cout << jzbcutprediction[ij] << " vs " << metcutprediction[ij] << " & ";
2689 buchmann 1.1
2690     cout << endl;
2691     //observation part
2692     if(is_data) cout << "Data observation & ";
2693     if(subselection!="none") cout << subselection << " observation &";
2694 buchmann 1.16 for(int ij=0;ij<(int)jzb_cuts.size();ij++) cout << jzbcutobservation[ij] << " vs " << metcutobservation[ij] << " & ";
2695 buchmann 1.1 cout << endl;
2696     cout << "_________________________________________________________________" << endl;
2697     delete cannie;
2698     }
2699    
2700     void met_jzb_cut(string datajzb, string mcjzb, vector<float> jzb_cut) {
2701     /*we want a table like this:
2702     __________________ 50 | 100 | ...
2703     | Data prediction | a vs b | a vs b | ...
2704     | Data observed | a vs b | a vs b | ...
2705     --------------------------------------
2706     --------------------------------------
2707     | LM4 prediction | a vs b | a vs b | ...
2708     | LM4 observed | a vs b | a vs b | ...
2709     | LM8 prediction | a vs b | a vs b | ...
2710     | LM8 observed | a vs b | a vs b | ...
2711    
2712     where a is the result for a jzb cut at X, and b is the result for a met cut at X
2713     */
2714     make_table(allsamples,datajzb, true,jzb_cut,"none");
2715     make_table(signalsamples,mcjzb, false,jzb_cut,"LM4");
2716     make_table(signalsamples,mcjzb, false,jzb_cut,"LM8");
2717     }
2718    
2719    
2720     //________________________________________________________________________________________
2721     // JZB Efficiency plot (efficiency of passing reco. JZB cut as function of generator JZB cut)
2722     void JZBSelEff(string mcjzb, TTree* events, string informalname, vector<float> jzb_bins) {
2723    
2724     float min = 0, max = 400;
2725     float xbins[] = {min,10,20,30,40,50,60,70,80,90,100,110,120,140,150,160,170,180,190,200,210,220,250,300,max,max+100};
2726     int nbins = sizeof(xbins)/sizeof(float)-1;
2727     int markers[] = { 20, 26, 21, 24, 22, 25, 28 };
2728    
2729 buchmann 1.14
2730     TH1F* heff = new TH1F("heff", "JZB eff ; generator JZB [GeV]; efficiency",nbins,xbins);
2731     TH1F* hgen = new TH1F("hgen", "JZB gen ; generator JZB [GeV]; efficiency",nbins,xbins);
2732     TH1F* hreco = new TH1F("hreco","JZB reco ; generator JZB [GeV]; efficiency",nbins,xbins);
2733 buchmann 1.1
2734     TCut kgen(genMassCut&&"genZPt>0&&genNjets>2&&abs(genMID)==23"&&cutOSSF);
2735     TCut kreco(cutmass);
2736    
2737     TF1* func = new TF1("func","0.5*[2]*(TMath::Erf((x-[1])/[0])+1)",min,max);
2738     func->SetParNames("epsilon","x_{1/2}","sigma");
2739     func->SetParameter(0,50.);
2740     func->SetParameter(1,0.);
2741     func->SetParameter(2,1.);
2742     gStyle->SetOptStat(0);
2743     gStyle->SetOptFit(0);
2744    
2745     TCanvas *can = new TCanvas("can","Canvas for JZB Efficiency",600,600);
2746     can->SetGridx(1);
2747     can->SetGridy(1);
2748 buchmann 1.14 can->SetLeftMargin(0.16);
2749     can->SetRightMargin(0.05);
2750 buchmann 1.1 TLegend *leg = make_legend("",0.6,0.2,false,0.89,0.5);
2751 buchmann 1.14 leg->SetBorderSize(1);
2752     leg->SetLineColor(kBlack);
2753     leg->SetTextFont(62);
2754 buchmann 1.1
2755 buchmann 1.16 for ( int icut=0; icut<(int)jzb_bins.size(); ++icut ) {
2756 buchmann 1.1
2757     ostringstream selection;
2758     selection << mcjzb << ">" << jzb_bins[icut];
2759     TCut ksel(selection.str().c_str());
2760     events->Draw("genJZB>>hgen", kgen&&kreco, "goff");
2761     events->Draw("genJZB>>hreco",kgen&&kreco&&ksel,"goff");
2762    
2763     // Loop over steps to get efficiency curve
2764     for ( Int_t iBin = 0; iBin<nbins-1; ++iBin ) {
2765     Float_t eff = hreco->GetBinContent(iBin+1)/hgen->GetBinContent(iBin+1);
2766     heff->SetBinContent(iBin+1,eff);
2767     heff->SetBinError(iBin+1,TMath::Sqrt(eff*(1-eff)/hgen->GetBinContent(iBin+1)));
2768     }
2769    
2770     heff->GetXaxis()->SetRangeUser(min, max);
2771 buchmann 1.14 // heff->GetXaxis()->SetLabelSize(0.05); // paper style. overruled.
2772     // heff->GetYaxis()->SetLabelSize(0.05); // paper style. overruled.
2773     // heff->GetXaxis()->SetTitleSize(0.06); // paper style. overruled.
2774     // heff->GetYaxis()->SetTitleSize(0.06); // paper style. overruled.
2775 buchmann 1.1 heff->SetMarkerStyle(markers[icut]);
2776     heff->Fit("func","Q+","same");
2777    
2778     // Print values
2779     dout << "+++ For " << selection.str() << std::endl;
2780     for ( int i=0; i<func->GetNpar(); ++i )
2781     dout << " " << func->GetParName(i) << " " << func->GetParameter(i) << " \\pm " << func->GetParError(i) << std::endl;
2782     char hname[256]; sprintf(hname,"heff%d",icut);
2783    
2784     // Store plot
2785     TH1F* h = (TH1F*)heff->Clone(hname);
2786 buchmann 1.14 h->SetNdivisions(505,"X");
2787 buchmann 1.1 if ( icut) h->Draw("same");
2788     else h->Draw();
2789     char htitle[256]; sprintf(htitle,"JZB > %3.0f GeV", jzb_bins[icut]);
2790     leg->AddEntry(h,htitle,"p");
2791    
2792     }
2793    
2794     leg->Draw();
2795     DrawMCPrelim(0.0);
2796     CompleteSave(can, "Systematics/jzb_efficiency_curve"+informalname );
2797    
2798     delete hgen;
2799     delete hreco;
2800     delete heff;
2801     }
2802    
2803     //________________________________________________________________________________________
2804     // Calls the above function for each signal sample
2805     void plot_jzb_sel_eff(string mcjzb, samplecollection &signalsamples, vector<float> bins )
2806     {
2807 buchmann 1.16 for (int isignal=0; isignal<(int)signalsamples.collection.size();isignal++) {
2808 buchmann 1.1 dout << "JZB selection efficiency curve: " << std::endl;
2809     JZBSelEff(mcjzb,(signalsamples.collection)[isignal].events,(signalsamples.collection)[isignal].samplename,bins); // Only for some selected samples
2810     }
2811     }
2812 buchmann 1.3
2813     void qcd_plots(string datajzb, string mcjzb, vector<float> bins) {
2814     // What this function aims to do:
2815     // Illustrate cut flow for QCD (requiring only one lepton, requiring etc.)
2816     // Illustrate how little QCD is left over! i.e. make some pointless JZB plot with only QCD to visualize the fact that there's not much really.
2817     TCanvas *can = new TCanvas("can","can");
2818     TPad *kinpad = new TPad("kinpad","kinpad",0,0,1,1);
2819     kinpad->cd();
2820    
2821     string jzb=mcjzb;
2822    
2823     float hi=400;
2824     bool use_signal=false;
2825     bool use_data=false;
2826    
2827     bool is_data=false;
2828     int nbins=50;//100;
2829     float low=0;
2830     float high=500;
2831     int versok=false;
2832     if(gROOT->GetVersionInt()>=53000) versok=true;
2833    
2834     TH1F *blankback = new TH1F("blankback","blankback",int(high/10),0,high);
2835     TH1F *RcorrJZBeemm = qcdsamples.Draw("RcorrJZBeemm",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,is_data, luminosity,use_signal);
2836     TH1F *LcorrJZBeemm = qcdsamples.Draw("LcorrJZBeemm",("-"+jzb).c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,is_data, luminosity,use_signal);
2837     TH1F *RcorrJZBem = qcdsamples.Draw("RcorrJZBem",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,is_data, luminosity,use_signal);
2838     TH1F *LcorrJZBem = qcdsamples.Draw("LcorrJZBem",("-"+jzb).c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,is_data, luminosity,use_signal);
2839     blankback->GetXaxis()->SetTitle(RcorrJZBeemm->GetXaxis()->GetTitle());
2840     blankback->GetYaxis()->SetTitle(RcorrJZBeemm->GetYaxis()->GetTitle());
2841     blankback->GetXaxis()->CenterTitle();
2842     blankback->GetYaxis()->CenterTitle();
2843    
2844     flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak ---- prediction changed.
2845     TH1F *RcorrJZBSBem;
2846     TH1F *LcorrJZBSBem;
2847     TH1F *RcorrJZBSBeemm;
2848     TH1F *LcorrJZBSBeemm;
2849    
2850 buchmann 1.16 // TH1F *RcorrJZBeemmNoS;
2851 buchmann 1.3
2852     //these are for the ratio
2853     TH1F *JRcorrJZBeemm = qcdsamples.Draw("JRcorrJZBeemm",jzb.c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,is_data, luminosity,use_signal);
2854     TH1F *JLcorrJZBeemm = qcdsamples.Draw("JLcorrJZBeemm",("-"+jzb).c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,is_data, luminosity,use_signal);
2855     TH1F *JRcorrJZBem = qcdsamples.Draw("JRcorrJZBem",jzb.c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,is_data, luminosity,use_signal);
2856     TH1F *JLcorrJZBem = qcdsamples.Draw("JLcorrJZBem",("-"+jzb).c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,is_data, luminosity,use_signal);
2857    
2858     TH1F *JRcorrJZBSBem;
2859     TH1F *JLcorrJZBSBem;
2860     TH1F *JRcorrJZBSBeemm;
2861     TH1F *JLcorrJZBSBeemm;
2862    
2863     if(PlottingSetup::RestrictToMassPeak) {
2864     RcorrJZBSBem = qcdsamples.Draw("RcorrJZBSBem",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets,is_data, luminosity,use_signal);
2865     LcorrJZBSBem = qcdsamples.Draw("LcorrJZBSBem",("-"+jzb).c_str(),nbins,low,hi, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets,is_data, luminosity,use_signal);
2866     RcorrJZBSBeemm = qcdsamples.Draw("RcorrJZBSBeemm",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets,is_data, luminosity,use_signal);
2867     LcorrJZBSBeemm = qcdsamples.Draw("LcorrJZBSBeemm",("-"+jzb).c_str(),nbins,low,hi, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets,is_data, luminosity,use_signal);
2868    
2869     //these are for the ratio
2870     JRcorrJZBSBem = qcdsamples.Draw("JRcorrJZBSBem",jzb.c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets,is_data, luminosity,use_signal);
2871     JLcorrJZBSBem = qcdsamples.Draw("JLcorrJZBSBem",("-"+jzb).c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets,is_data, luminosity,use_signal);
2872     JRcorrJZBSBeemm = qcdsamples.Draw("JRcorrJZBSBeemm",jzb.c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets,is_data, luminosity,use_signal);
2873     JLcorrJZBSBeemm = qcdsamples.Draw("JLcorrJZBSBeemm",("-"+jzb).c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets,is_data, luminosity,use_signal);
2874    
2875     }
2876    
2877     flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak ---- prediction changed.
2878    
2879     TH1F *Zjetspred = (TH1F*)LcorrJZBeemm->Clone("Zjetspred");
2880     TH1F *TTbarpred = (TH1F*)RcorrJZBem->Clone("TTbarpred");
2881    
2882     TH1F *Bpred = (TH1F*)LcorrJZBeemm->Clone("Bpred");
2883     TH1F *JBpred = (TH1F*)JLcorrJZBeemm->Clone("Bpred");
2884     if(PlottingSetup::RestrictToMassPeak) {
2885     Bpred->Add(RcorrJZBem,1.0/3.);
2886     Bpred->Add(LcorrJZBem,-1.0/3.);
2887     Bpred->Add(RcorrJZBSBem,1.0/3.);
2888     Bpred->Add(LcorrJZBSBem,-1.0/3.);
2889     Bpred->Add(RcorrJZBSBeemm,1.0/3.);
2890     Bpred->Add(LcorrJZBSBeemm,-1.0/3.);
2891    
2892     TTbarpred->Scale(1.0/3);
2893     Zjetspred->Add(LcorrJZBem,-1.0/3.);
2894     Zjetspred->Add(LcorrJZBSBem,-1.0/3.);
2895     TTbarpred->Add(RcorrJZBSBem,1.0/3.);
2896     Zjetspred->Add(LcorrJZBSBeemm,-1.0/3.);
2897     TTbarpred->Add(RcorrJZBSBeemm,1.0/3.);
2898    
2899     //these are for the ratio
2900     JBpred->Add(JRcorrJZBem,1.0/3.);
2901     JBpred->Add(JLcorrJZBem,-1.0/3.);
2902     JBpred->Add(JRcorrJZBSBem,1.0/3.);
2903     JBpred->Add(JLcorrJZBSBem,-1.0/3.);
2904     JBpred->Add(JRcorrJZBSBeemm,1.0/3.);
2905     JBpred->Add(JLcorrJZBSBeemm,-1.0/3.);
2906     } else {
2907     Bpred->Add(RcorrJZBem,1.0);
2908     Bpred->Add(LcorrJZBem,-1.0);
2909    
2910     Zjetspred->Add(LcorrJZBem,-1.0);
2911    
2912     //these are for the ratio
2913     JBpred->Add(JRcorrJZBem,1.0);
2914     JBpred->Add(JLcorrJZBem,-1.0);
2915     }
2916    
2917    
2918     flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak ---- prediction changed
2919 buchmann 1.1
2920 buchmann 1.3 TLegend *legBpred = make_legend("",0.6,0.55,false);
2921     RcorrJZBeemm->Draw("e1x0,same");
2922     Bpred->Draw("hist,same");
2923     RcorrJZBeemm->Draw("e1x0,same");//HAVE IT ON TOP!
2924     legBpred->AddEntry(RcorrJZBeemm,"MC observed","p");
2925     legBpred->AddEntry(Bpred,"MC predicted","l");
2926     if(versok) legBpred->AddEntry((TObject*)0,"",""); // Just for alignment // causes seg fault on root v5.18
2927     if(versok) legBpred->AddEntry((TObject*)0,"",""); // causes seg fault on root v5.18
2928     legBpred->Draw();
2929     DrawMCPrelim();
2930    
2931     //3rd last argument: do special bpred ratio, 2nd last argument: extended range!, last: y-axis title
2932     string ytitle("ratio");
2933     if ( use_data==1 ) ytitle = "data/pred";
2934     save_with_ratio(JRcorrJZBeemm,JBpred,kinpad,"QCD/Bpred",true,use_data!=1,ytitle);
2935    
2936     TH1F *allevents = qcdsamples.Draw("allevents","pfJetGoodNum",1,0,100, "internal code", "events", "" ,mc, luminosity);
2937     TH1F *ossf = qcdsamples.Draw("ossf","pfJetGoodNum",1,0,100, "internal code", "events", cutOSSF ,mc, luminosity);
2938     TH1F *osof = qcdsamples.Draw("osof","pfJetGoodNum",1,0,100, "internal code", "events", cutOSOF ,mc, luminosity);
2939     TH1F *njossf = qcdsamples.Draw("njossf","pfJetGoodNum",1,0,100, "internal code", "events", cutnJets&&cutOSSF ,mc, luminosity);
2940     TH1F *njosof = qcdsamples.Draw("njosof","pfJetGoodNum",1,0,100, "internal code", "events", cutnJets&&cutOSOF ,mc, luminosity);
2941    
2942     dout << "______________________________________________" << endl;
2943     dout << "QCD contribution: " << endl;
2944     dout << "Total number of events: " << allevents->Integral() << endl;
2945     dout << "OSSF events: " << ossf->Integral() << endl;
2946     dout << "OSOF events: " << osof->Integral() << endl;
2947     dout << "OSSF events with >=3 jets:" << njossf->Integral() << endl;
2948     dout << "OSOF events with >=3 jets:" << njosof->Integral() << endl;
2949     dout << "(note that no mass requirement has been imposed)" << endl;
2950    
2951     dout << "______________________________________________" << endl;
2952     dout << "How QCD shows up in the different regions: " << endl;
2953     dout << "OSSF: " << endl;
2954     if(PlottingSetup::RestrictToMassPeak) {
2955     dout << " Z window: \t" << RcorrJZBeemm->Integral() << " (JZB>0) , " << LcorrJZBeemm->Integral() << " (JZB<0) --> total: " << RcorrJZBeemm->Integral() + LcorrJZBeemm->Integral() << endl;
2956     dout << " sideband: \t" << RcorrJZBSBeemm->Integral() << " (JZB>0) , " << LcorrJZBSBeemm->Integral() << " (JZB<0) --> total: " << RcorrJZBSBeemm->Integral() + LcorrJZBSBeemm->Integral() << endl;
2957     } else {
2958     dout << " " << RcorrJZBeemm->Integral() << " (JZB>0) , " << LcorrJZBeemm->Integral() << " (JZB<0) --> total: " << RcorrJZBeemm->Integral() + LcorrJZBeemm->Integral() << endl;
2959     }
2960     dout << "OSOF: " << endl;
2961     if(PlottingSetup::RestrictToMassPeak) {
2962     dout << " Z window: \t" << RcorrJZBem->Integral() << " (JZB>0) , " << LcorrJZBem->Integral() << " (JZB<0) --> total: " << RcorrJZBem->Integral() + LcorrJZBem->Integral() << endl;
2963     dout << " sideband: \t" << RcorrJZBSBem->Integral() << " (JZB>0) , " << LcorrJZBSBem->Integral() << " (JZB<0) --> total: " << RcorrJZBSBem->Integral() + LcorrJZBSBem->Integral() << endl;
2964     } else {
2965     dout << " Z window: \t" << RcorrJZBem->Integral() << " (JZB>0) , " << LcorrJZBem->Integral() << " (JZB<0) --> total: " << RcorrJZBem->Integral() + LcorrJZBem->Integral() << endl;
2966     }
2967     dout << "Therefore: " << endl;
2968     if(PlottingSetup::RestrictToMassPeak) {
2969     dout << " Prediction increases by : " << LcorrJZBeemm->Integral() << " + (1.0/3)*(" << RcorrJZBSBeemm->Integral() <<"-"<< LcorrJZBSBeemm->Integral() << ") (SFSB) ";
2970     dout << " + (1.0/3)*(" << RcorrJZBem->Integral() <<"-"<< LcorrJZBem->Integral() << ") (OFZP) ";
2971     dout << " + (1.0/3)*(" << RcorrJZBSBem->Integral() <<"-"<< LcorrJZBSBem->Integral() << ") (OFSB) ";
2972     dout << " = " << LcorrJZBeemm->Integral() + (1.0/3)*(RcorrJZBSBeemm->Integral() - LcorrJZBSBeemm->Integral() + RcorrJZBem->Integral() - LcorrJZBem->Integral() + RcorrJZBSBem->Integral() - LcorrJZBSBem->Integral()) << endl;
2973     } else {
2974     dout << " Prediction increases by : " << LcorrJZBeemm->Integral();
2975     dout << " + (" << RcorrJZBem->Integral() <<"-"<< LcorrJZBem->Integral() << ") (OFZP) ";
2976     dout << " = " << LcorrJZBeemm->Integral() + RcorrJZBem->Integral() - LcorrJZBem->Integral() << endl;
2977     }
2978     dout << " Observation increases by : " << RcorrJZBeemm->Integral() << endl;
2979    
2980     dout << endl;
2981 buchmann 1.16 for(int i=0;i<(int)bins.size();i++) {
2982 buchmann 1.3 dout << " JZB > " << bins[i] << " : " << endl;
2983     dout << " Observation increases by : " << RcorrJZBeemm->Integral(RcorrJZBeemm->FindBin(bins[i]),RcorrJZBeemm->GetNbinsX()) << endl;
2984     if(PlottingSetup::RestrictToMassPeak) {
2985     dout << " Prediction increases by : " << LcorrJZBeemm->Integral(RcorrJZBeemm->FindBin(bins[i]),RcorrJZBeemm->GetNbinsX()) + (1.0/3)*(RcorrJZBSBeemm->Integral(RcorrJZBeemm->FindBin(bins[i]),RcorrJZBeemm->GetNbinsX()) - LcorrJZBSBeemm->Integral(RcorrJZBeemm->FindBin(bins[i]),RcorrJZBeemm->GetNbinsX()) + RcorrJZBem->Integral(RcorrJZBeemm->FindBin(bins[i]),RcorrJZBeemm->GetNbinsX()) - LcorrJZBem->Integral(RcorrJZBeemm->FindBin(bins[i]),RcorrJZBeemm->GetNbinsX()) + RcorrJZBSBem->Integral(RcorrJZBeemm->FindBin(bins[i]),RcorrJZBeemm->GetNbinsX()) - LcorrJZBSBem->Integral(RcorrJZBeemm->FindBin(bins[i]),RcorrJZBeemm->GetNbinsX())) << endl;
2986     } else {
2987     dout << " Prediction increases by : " << LcorrJZBeemm->Integral(RcorrJZBeemm->FindBin(bins[i]),RcorrJZBeemm->GetNbinsX()) + RcorrJZBem->Integral(RcorrJZBeemm->FindBin(bins[i]),RcorrJZBeemm->GetNbinsX()) - LcorrJZBem->Integral(RcorrJZBeemm->FindBin(bins[i]),RcorrJZBeemm->GetNbinsX()) << endl;
2988     }
2989     }
2990    
2991     delete can;
2992     delete allevents;
2993     if(ossf) delete ossf;
2994     if(RcorrJZBem) delete RcorrJZBem;
2995     if(LcorrJZBem) delete LcorrJZBem;
2996     if(RcorrJZBeemm) delete RcorrJZBeemm;
2997     if(LcorrJZBeemm) delete LcorrJZBeemm;
2998     if(PlottingSetup::RestrictToMassPeak&&RcorrJZBSBem) delete RcorrJZBSBem;
2999     if(PlottingSetup::RestrictToMassPeak&&LcorrJZBSBem) delete LcorrJZBSBem;
3000     if(PlottingSetup::RestrictToMassPeak&&RcorrJZBSBeemm) delete RcorrJZBSBeemm;
3001     if(PlottingSetup::RestrictToMassPeak&&LcorrJZBSBeemm) delete LcorrJZBSBeemm;
3002     }
3003 buchmann 1.1
3004 buchmann 1.4 void check_ptsanity() {
3005     TCanvas *ptsancan = new TCanvas("ptsancan","ptsancan",600,1800);
3006     TH1F *individualpt1histos[allsamples.collection.size()];
3007     TH1F *individualpt2histos[allsamples.collection.size()];
3008     TH1F *fpt1 = new TH1F("fpt1","fpt1",50,0,50);
3009     fpt1->GetYaxis()->SetRangeUser(0,1);
3010     fpt1->GetXaxis()->SetTitle("p_{T,1}");
3011     fpt1->GetXaxis()->CenterTitle();
3012    
3013     TH1F *fpt2 = new TH1F("fpt2","fpt2",50,0,50);
3014     fpt2->GetXaxis()->SetTitle("p_{T,2}");
3015     fpt2->GetXaxis()->CenterTitle();
3016    
3017     ptsancan->Divide(1,3);
3018     ptsancan->cd(1);
3019     float maxpt1entry=0;
3020     float maxpt2entry=0;
3021    
3022     TLegend *leg = make_legend();
3023     leg->SetX1(0.0);
3024     leg->SetY1(0.0);
3025     leg->SetX2(1.0);
3026     leg->SetY2(1.0);
3027    
3028    
3029 buchmann 1.16 for(int isample=0;isample<(int)allsamples.collection.size();isample++) {
3030 buchmann 1.4 string nowname=(allsamples.collection)[isample].filename;
3031     cout << "Drawing: " << nowname << " (sample " << isample+1 << " / " << allsamples.collection.size() << ")" << endl;
3032     individualpt1histos[isample] = allsamples.Draw(GetNumericHistoName(),"pt1",50,0,50, "p_{T,1}", "events",cutOSSF&&cutnJets,mc,luminosity,allsamples.FindSample(nowname));
3033     individualpt2histos[isample] = allsamples.Draw(GetNumericHistoName(),"pt2",50,0,50, "p_{T,2}", "events",cutOSSF&&cutnJets,mc,luminosity,allsamples.FindSample(nowname));
3034     individualpt1histos[isample]->SetLineColor(isample+1);
3035     individualpt2histos[isample]->SetLineColor(isample+1);
3036     float currmaxpt1entry=individualpt1histos[isample]->GetMaximum()/individualpt1histos[isample]->Integral();
3037     float currmaxpt2entry=individualpt2histos[isample]->GetMaximum()/individualpt2histos[isample]->Integral();
3038     cout << " pt 1 histo contains; " << individualpt1histos[isample]->Integral() << endl;
3039     cout << " pt 2 histo contains; " << individualpt2histos[isample]->Integral() << endl;
3040     if(currmaxpt1entry>maxpt1entry)maxpt1entry=currmaxpt1entry;
3041     if(currmaxpt2entry>maxpt2entry)maxpt2entry=currmaxpt2entry;
3042     leg->AddEntry(individualpt2histos[isample],((allsamples.collection)[isample].filename).c_str(),"f");
3043     }
3044    
3045     fpt1->GetYaxis()->SetRangeUser(0,maxpt1entry);
3046     fpt2->GetYaxis()->SetRangeUser(0,maxpt2entry);
3047    
3048     ptsancan->cd(1);
3049     fpt1->Draw();
3050     ptsancan->cd(2);
3051     fpt2->Draw();
3052    
3053 buchmann 1.16 for(int isample=0;isample<(int)allsamples.collection.size();isample++) {
3054 buchmann 1.4 ptsancan->cd(1);
3055     individualpt1histos[isample]->DrawNormalized("same,histo");
3056     ptsancan->cd(2);
3057     individualpt2histos[isample]->DrawNormalized("same,histo");
3058     }
3059     ptsancan->cd(3);
3060     leg->Draw();
3061     CompleteSave(ptsancan,"PtSanityCheck");
3062    
3063     delete ptsancan;
3064     }
3065    
3066 buchmann 1.6 void do_mlls_plot(string mcjzb) {
3067     cout << "At this point we'd plot the mll distribution" << endl;
3068     TCanvas *sigcan = new TCanvas("sigcan","sigcan");
3069 buchmann 1.16 for(int isig=0;isig<(int)(signalsamples.collection).size();isig++) {
3070 buchmann 1.6 if(!(signalsamples.collection)[isig].events) continue;
3071     string nowname=(signalsamples.collection)[isig].filename;
3072     TH1F *mll = signalsamples.Draw("mllhisto","mll",150,0,150, "m_{ll}", "events",cutOSSF&&cutnJets,mc,luminosity,signalsamples.FindSample(nowname));
3073     // TH1F *mll = signalsamples.Draw("mllhisto","mll",150,0,150, "m_{ll}", "events","",mc,luminosity,signalsamples.FindSample(nowname));
3074     mll->SetLineColor(TColor::GetColor("#04B404"));
3075     stringstream poscutS;
3076     poscutS << "((" << mcjzb <<")>50)";
3077     TCut poscut(poscutS.str().c_str());
3078     TH1F *mllP = signalsamples.Draw("mllhistoP","mll",150,0,150, "m_{ll}", "events",cutOSSF&&cutnJets&&poscut,mc,luminosity,signalsamples.FindSample(nowname));
3079     mllP->SetLineColor(TColor::GetColor("#0040FF"));
3080     mll->Draw("histo");
3081     mllP->Draw("histo,same");
3082     TLegend *leg = make_legend();
3083     leg->SetY1(0.8);
3084     leg->AddEntry(mll,(signalsamples.collection)[isig].samplename.c_str(),"L");
3085     leg->AddEntry(mllP,((signalsamples.collection)[isig].samplename+", JZB>50").c_str(),"L");
3086     leg->Draw();
3087     TLine *lin = new TLine(71.2,0,71.2,mll->GetMaximum());
3088     TLine *lin2 = new TLine(111.2,0,111.2,mll->GetMaximum());
3089     lin->Draw("same");
3090     lin2->Draw("same");
3091    
3092     CompleteSave(sigcan,"MllShape/"+(signalsamples.collection)[isig].samplename);
3093     delete mll;
3094     delete mllP;
3095     }
3096     }
3097    
3098 buchmann 1.9 void met_vs_jzb_plots() {
3099    
3100     TCanvas *canmetjzb = new TCanvas("canmet","MET vs JZB canvas");
3101     canmetjzb->SetRightMargin(0.16);
3102    
3103     vector<string> findme;
3104     findme.push_back("DY");
3105     findme.push_back("TTJets");
3106     findme.push_back("LM");
3107    
3108 buchmann 1.16 for(int ifind=0;ifind<(int)findme.size();ifind++) {
3109 buchmann 1.9 vector<int> selsamples = allsamples.FindSample(findme[ifind]);
3110     TH2F *metvsjzb = new TH2F("metvsjzb","metvsjzb",200,0,100,400,-100,100);
3111 buchmann 1.16 for(int isel=0;isel<(int)selsamples.size();isel++) {
3112 buchmann 1.9 cout << "Producing MET:JZB plot ... working on sample: " << allsamples.collection[selsamples[isel]].filename << endl;
3113     allsamples.collection[selsamples[isel]].events->Draw("jzb[1]:met[4]>>+metvsjzb",cutmass&&cutOSSF);
3114     }
3115     metvsjzb->Scale(allsamples.collection[selsamples[0]].weight);
3116     metvsjzb->SetStats(0);
3117     metvsjzb->GetXaxis()->SetTitle("MET (GeV)");
3118     metvsjzb->GetYaxis()->SetTitle("JZB (GeV)");
3119     metvsjzb->GetXaxis()->CenterTitle();
3120     metvsjzb->GetYaxis()->CenterTitle();
3121     metvsjzb->Draw("COLZ");
3122     TText* title = write_text(0.5,0.95,allsamples.collection[selsamples[0]].samplename);
3123     title->SetTextAlign(12);
3124     title->Draw();
3125     CompleteSave(canmetjzb,(string)"METvsJZBplots/"+findme[ifind]);
3126     }
3127     }
3128    
3129    
3130 buchmann 1.1 void test() {
3131    
3132     TCanvas *testcanv = new TCanvas("testcanv","testcanv");
3133     testcanv->cd();
3134     switch_overunderflow(true);
3135     TH1F *ptdistr = allsamples.Draw("ptdistr","pt1",100,30,200, "p_{T} [GeV]", "events", cutOSSF,data,luminosity);
3136     switch_overunderflow(false);
3137     ptdistr->Draw();
3138     testcanv->SaveAs("test.png");
3139     dout << "HELLO there!" << endl;
3140    
3141     }