ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/ExperimentalModule.C
Revision: 1.3
Committed: Wed Aug 15 13:52:12 2012 UTC (12 years, 8 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.2: +1 -1 lines
Log Message:
Added option to switch off sidebands for classic JZB

File Contents

# User Rev Content
1 buchmann 1.1 /****
2    
3     Off peak status (RestrictToMassPeak) :
4    
5     x Necessary adaptations identified
6     x Started working on necessary adaptations
7     x Necessary adaptations implemented
8     x Necessary adaptations tested
9    
10     DONE!
11    
12     ****/
13    
14     #include <iostream>
15     #include <vector>
16     #include <sys/stat.h>
17    
18     #include <TCut.h>
19     #include <TROOT.h>
20     #include <TCanvas.h>
21     #include <TMath.h>
22     #include <TColor.h>
23     #include <TPaveText.h>
24     #include <TRandom.h>
25     #include <TH1.h>
26     #include <TH2.h>
27     #include <TF1.h>
28     #include <TSQLResult.h>
29     #ifndef GeneralToolBoxLoaded
30     #include "GeneralToolBox.C"
31     #endif
32    
33     using namespace std;
34    
35     using namespace PlottingSetup;
36    
37     //--------------------------------------------- below: ttbar
38    
39     void ttbar_limit_shapes_for_systematic_effect(TFile *limfile, string identifier, string mcjzb, string datajzb, int JES,vector<float> binning, TCanvas *limcan) {
40     dout << "Creatig shape templates ... ";
41     if(identifier!="") dout << "for systematic called "<<identifier;
42     dout << endl;
43     int dataormc=mc;//this is only for tests - for real life you want dataormc=data !!!
44     if(dataormc!=data) write_warning("limit_shapes_for_systematic_effect","WATCH OUT! Not using data for limits!!!! this is ok for tests, but not ok for anything official!");
45    
46     TCut limitnJetcut;
47     if(JES==noJES) limitnJetcut=cutnJets;
48     else {
49     if(JES==JESdown) limitnJetcut=cutnJetsJESdown;
50     if(JES==JESup) limitnJetcut=cutnJetsJESup;
51     }
52     TH1F *ZOSSFP = allsamples.Draw("ZOSSFP",datajzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,dataormc,luminosity);
53     TH1F *ZOSOFP = allsamples.Draw("ZOSOFP",datajzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,dataormc,luminosity);
54     TH1F *ZOSSFN = allsamples.Draw("ZOSSFN","-"+datajzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,dataormc,luminosity);
55     TH1F *ZOSOFN = allsamples.Draw("ZOSOFN","-"+datajzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,dataormc,luminosity);
56    
57     TH1F *LZOSSFP = allsamples.Draw("LZOSSFP",mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("TTJets"));
58     TH1F *LZOSOFP = allsamples.Draw("LZOSOFP",mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("TTJets"));
59     TH1F *LZOSSFN = allsamples.Draw("LZOSSFN","-"+mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("TTJets"));
60     TH1F *LZOSOFN = allsamples.Draw("LZOSOFN","-"+mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("TTJets"));
61    
62     string obsname="data_obs";
63     string predname="background";
64     string signalname="signal";
65     if(identifier!="") {
66     obsname=("data_"+identifier);
67     predname=("background_"+identifier);
68     signalname="signal_"+identifier;
69     }
70    
71     TH1F *obs = (TH1F*)ZOSSFP->Clone();
72     obs->SetName(obsname.c_str());
73     obs->Write();
74     TH1F *pred = (TH1F*)ZOSSFN->Clone();
75     pred->Add(ZOSOFN,-1.0);
76     pred->SetName(predname.c_str());
77     pred->Write();
78    
79     TH1F *Lobs = new TH1F("Lobs","Lobs",binning.size()-1,&binning[0]);
80     TH1F *Lpred = new TH1F("Lpred","Lpred",binning.size()-1,&binning[0]);
81     Lobs->Add(LZOSSFP);
82     Lpred->Add(LZOSSFN);
83     TH1F *signal = (TH1F*)Lobs->Clone();
84     signal->Add(Lpred,-1);
85     signal->SetName(signalname.c_str());
86     signal->Write();
87    
88     delete Lobs;
89     delete Lpred;
90    
91     delete ZOSSFP;
92     delete ZOSOFP;
93     delete ZOSSFN;
94     delete ZOSOFN;
95    
96     delete LZOSSFP;
97     delete LZOSOFP;
98     delete LZOSSFN;
99     delete LZOSOFN;
100    
101     }
102    
103     void prepare_ttbar_datacard(TFile *f) {
104     TH1F *dataob = (TH1F*)f->Get("data_obs");
105     TH1F *signal = (TH1F*)f->Get("signal");
106     TH1F *background = (TH1F*)f->Get("background");
107    
108     ofstream datacard;
109     ensure_directory_exists(get_directory()+"/limits");
110     datacard.open ((get_directory()+"/limits/ttbardatacard.txt").c_str());
111     datacard << "Writing this to a file.\n";
112     datacard << "imax 1\n";
113     datacard << "jmax 1\n";
114     datacard << "kmax *\n";
115     datacard << "---------------\n";
116     datacard << "shapes * * ttbarlimitfile.root $PROCESS $PROCESS_$SYSTEMATIC\n";
117     datacard << "---------------\n";
118     datacard << "bin 1\n";
119     datacard << "observation "<<dataob->Integral()<<"\n";
120     datacard << "------------------------------\n";
121     datacard << "bin 1 1\n";
122     datacard << "process signal background\n";
123     datacard << "process 0 1\n";
124     datacard << "rate "<<signal->Integral()<<" "<<background->Integral()<<"\n";
125     datacard << "--------------------------------\n";
126     datacard << "lumi lnN 1.10 1.0\n";
127     datacard << "bgnorm lnN 1.00 1.4 uncertainty on our prediction (40%)\n";
128     datacard << "JES shape 1 1 uncertainty on background shape and normalization\n";
129     datacard << "peak shape 1 1 uncertainty on signal resolution. Assume the histogram is a 2 sigma shift, \n";
130     datacard << "# so divide the unit gaussian by 2 before doing the interpolation\n";
131     datacard.close();
132    
133    
134    
135     }
136    
137     void prepare_ttbar_limits(string mcjzb, string datajzb, float jzbpeakerrordata, float jzbpeakerrormc, vector<float> jzbbins) {
138     ensure_directory_exists(get_directory()+"/limits");
139     TFile *limfile = new TFile((get_directory()+"/limits/ttbarlimitfile.root").c_str(),"RECREATE");
140     TCanvas *limcan = new TCanvas("limcan","Canvas for calculating limits");
141     ttbar_limit_shapes_for_systematic_effect(limfile,"",mcjzb,datajzb,noJES,jzbbins,limcan);
142     ttbar_limit_shapes_for_systematic_effect(limfile,"peakUp",newjzbexpression(mcjzb,jzbpeakerrormc),newjzbexpression(datajzb,jzbpeakerrordata),noJES,jzbbins,limcan);
143     ttbar_limit_shapes_for_systematic_effect(limfile,"peakDown",newjzbexpression(mcjzb,-jzbpeakerrormc),newjzbexpression(datajzb,-jzbpeakerrordata),noJES,jzbbins,limcan);
144     ttbar_limit_shapes_for_systematic_effect(limfile,"JESUp",mcjzb,datajzb,JESup,jzbbins,limcan);
145     ttbar_limit_shapes_for_systematic_effect(limfile,"JESDown",mcjzb,datajzb,JESdown,jzbbins,limcan);
146    
147     prepare_ttbar_datacard(limfile);
148    
149 buchmann 1.2 write_info(__FUNCTION__,"limitfile.root and datacard.txt have been generated. You can now use them to calculate limits!");
150 buchmann 1.1 limfile->Close();
151    
152     }
153    
154     void Poisson_ratio_plot(int is_data,vector<float> binning, string jzb, TCanvas *can, float high=-9999, float jzbpearerrorMC=-999, float jzbpeakerrorData=-999) {
155     can->SetLogy(0);
156     cout << "About to prepare ratio plots with Poisson errors and obs/pred plots, in the new style!" << endl;
157     TH1F *pred = new TH1F("pred","",binning.size()-1,&binning[0]);
158     TH1F *obs = new TH1F("obs","",binning.size()-1,&binning[0]);
159     for(int ibin=0;ibin<binning.size()-1;ibin++) {
160     // cout << "________________________________________________________________________" << endl;
161     cout << "Ratio plot: Working on bin " << ibin+1 << " of " << binning.size()-1 << " for step " << is_data+1 << " of 3" << endl;
162     float binerr=0;
163     vector<float> contentanderror=get_result_between_two_fixed_jzb_values(true,binning[ibin],binning[ibin+1], jzb,jzb, is_data,jzbpearerrorMC, jzbpeakerrorData, can,false, false,false);
164     //[0] Bpred [1] Bpred uncert [2] Observed [3] Observed error
165     // cout << "Bpred : " << contentanderror[0] << " +/- " << contentanderror[1] << endl;
166     // cout << "Obs : " << contentanderror[2] << " +/- " << contentanderror[3] << endl;
167     pred->SetBinContent(ibin+1,contentanderror[0]);
168     pred->SetBinError(ibin+1,contentanderror[1]);
169     obs->SetBinContent(ibin+1,contentanderror[2]);
170     obs->SetBinError(ibin+1,contentanderror[3]);
171     }
172     TH1F *JZBratioforfitting=(TH1F*)obs->Clone("JZBratioforfitting");
173     JZBratioforfitting->Divide(pred);
174     TGraphAsymmErrors *JZBratio = histRatio(obs,pred,is_data,binning,true);//the last argument basically forces the function to give us the real error, instead of the stat one
175    
176     JZBratio->SetTitle("");
177    
178     TF1 *pol0 = new TF1("pol0","[0]",0,1000);
179     TF1 *pol0d = new TF1("pol0","[0]",0,1000);
180     JZBratioforfitting->Fit(pol0,"Q0R","",0,30);
181     pol0d->SetParameter(0,pol0->GetParameter(0));
182    
183     can->cd();
184     JZBratio->GetYaxis()->SetTitle("Observed/Predicted");
185     JZBratio->GetXaxis()->SetTitle("JZB [GeV]");
186     if ( high>0 ) JZBratio->GetXaxis()->SetRangeUser(0.0,high);
187     JZBratio->GetYaxis()->SetNdivisions(519);
188     JZBratio->GetYaxis()->SetRangeUser(0.0,4.0);
189     JZBratio->GetYaxis()->CenterTitle();
190     JZBratio->GetXaxis()->CenterTitle();
191     JZBratio->SetMarkerSize(DataMarkerSize);
192     JZBratio->Draw("AP");
193     /////----------------------------
194     TPaveText *writeresult = new TPaveText(0.15,0.78,0.49,0.91,"blNDC");
195     writeresult->SetFillStyle(4000);
196     writeresult->SetFillColor(kWhite);
197     writeresult->SetTextFont(42);
198     writeresult->SetTextSize(0.03);
199     writeresult->SetTextAlign(12);
200     ostringstream jzb_agreement_data_text;
201     jzb_agreement_data_text<< setprecision(2) << "mean =" << pol0->GetParameter(0) << " #pm " << setprecision(1) << pol0->GetParError(0);
202     if(is_data==1) fitresultconstdata=pol0->GetParameter(0);// data
203     if(is_data==0) fitresultconstmc=pol0->GetParameter(0); // monte carlo, no signal
204    
205     writeresult->AddText(jzb_agreement_data_text.str().c_str());
206     // writeresult->Draw("same");
207     // pol0d->Draw("same");
208     TF1 *topline = new TF1("","1.5",0,1000);
209     TF1 *bottomline = new TF1("","0.5",0,1000);
210     topline->SetLineColor(kBlue);
211     topline->SetLineStyle(2);
212     bottomline->SetLineColor(kBlue);
213     bottomline->SetLineStyle(2);
214     // topline->Draw("same");
215     // bottomline->Draw("same");
216     TF1 *oneline = new TF1("","1.0",0,1000);
217     oneline->SetLineColor(kBlue);
218     oneline->SetLineStyle(1);
219     oneline->Draw("same");
220     if(is_data==1) DrawPrelim();
221     else DrawMCPrelim();
222     TLegend *leg = new TLegend(0.5,0.55,0.94,0.89);
223     leg->SetTextFont(42);
224     leg->SetTextSize(0.04);
225    
226     TString MCtitle("MC ");
227     if (is_data==1) MCtitle = "";
228    
229     leg->SetFillStyle(4000);
230     leg->SetFillColor(kWhite);
231     leg->SetTextFont(42);
232     leg->AddEntry(JZBratio,MCtitle+"obs / "+MCtitle+"pred","p");
233     leg->AddEntry(oneline,"ratio = 1","l");
234     leg->AddEntry(pol0d,"fit in [0,30] GeV","l");
235     leg->AddEntry(bottomline,"#pm50% envelope","l");
236     // leg->Draw("same");
237     if(is_data==1) CompleteSave(can, "precise_jzb_ratio_data");
238     if(is_data==0) CompleteSave(can, "precise_jzb_ratio_mc");
239     if(is_data==2) CompleteSave(can, "precise_jzb_ratio_mc_BandS");//special case, MC with signal!
240    
241     can->SetLogy(1);
242     if(is_data==1) DrawPrelim();
243     else DrawMCPrelim();
244     TLegend *secondleg = make_legend("Observed vs Predicted");
245     secondleg->AddEntry(pred,"Predicted (stat&syst band)","f");
246     secondleg->AddEntry(obs,"Observed","p");
247     secondleg->SetY1(0.7);
248     secondleg->SetX1(0.4);
249    
250     pred->GetXaxis()->SetTitle("JZB [GeV]");
251     pred->GetXaxis()->CenterTitle();
252     pred->GetYaxis()->SetTitle("events");
253     pred->GetYaxis()->CenterTitle();
254     pred->SetFillColor(kGreen);
255     pred->SetMarkerColor(kWhite);
256     pred->SetMarkerSize(0.000001);
257     pred->SetLineColor(pred->GetFillColor());
258     pred->Draw("e3");
259     DrawPrelim();
260     obs->Draw("e1x0,same");
261     secondleg->Draw("same");
262    
263     if(is_data==1) CompleteSave(can, "Exp/precise_obs_pred_data");
264     if(is_data==0) CompleteSave(can, "Exp/precise_obs_pred_mc");
265     if(is_data==2) CompleteSave(can, "Exp/precise_obs_pred_mc_BandS");//special case, MC with signal!
266    
267     delete obs;
268     delete pred;
269     }
270    
271     void Poisson_ratio_plots(string mcjzb,string datajzb,vector<float> ratio_binning,float jzbpearerrorMC,float jzbpeakerrorData) {
272     TCanvas *globalc = new TCanvas("globalc","Ratio Plot Canvas");
273     globalc->SetLogy(0);
274    
275     Poisson_ratio_plot(data,ratio_binning,datajzb,globalc, jzbHigh,jzbpearerrorMC,jzbpeakerrorData);
276     Poisson_ratio_plot(mc,ratio_binning,mcjzb,globalc, jzbHigh,jzbpearerrorMC,jzbpeakerrorData);
277     Poisson_ratio_plot(mcwithsignal,ratio_binning,mcjzb,globalc, jzbHigh,jzbpearerrorMC,jzbpeakerrorData);
278     }
279    
280    
281    
282    
283    
284     TF1* do_new_logpar_fit_to_plot(TH1F *osof) {
285     TCanvas *logpar_fit_can = new TCanvas("logpar_fit_can","Fit canvas for LogPar");
286     TF1 *logparfunc = new TF1("logparfunc",LogParabola,60,200,3);
287     TF1 *logparfunc2 = new TF1("logparfunc2",LogParabola,0,(osof->GetXaxis()->GetBinLowEdge(osof->GetNbinsX())+osof->GetXaxis()->GetBinWidth(osof->GetNbinsX())),3);
288     TF1 *logparfuncN = new TF1("logparfuncN",LogParabolaN,0,(osof->GetXaxis()->GetBinLowEdge(osof->GetNbinsX())+osof->GetXaxis()->GetBinWidth(osof->GetNbinsX())),3);
289     TF1 *logparfuncP = new TF1("logparfuncP",LogParabolaP,0,(osof->GetXaxis()->GetBinLowEdge(osof->GetNbinsX())+osof->GetXaxis()->GetBinWidth(osof->GetNbinsX())),3);
290     osof->SetMinimum(0);
291     osof->Fit(logparfunc,"QR");
292     osof->Draw();
293     logparfunc->SetLineWidth(2);
294     logparfunc2->SetParameters(logparfunc->GetParameters());
295     logparfuncN->SetParameters(logparfunc->GetParameters());
296     logparfuncP->SetParameters(logparfunc->GetParameters());
297     stringstream fitinfo;
298     fitinfo << "#Chi^{2} / NDF : " << logparfunc->GetChisquare() << " / " << logparfunc->GetNDF();
299     TText *writefitinfo = write_text(0.8,0.8,fitinfo.str());
300     writefitinfo->SetTextSize(0.03);
301     DrawPrelim();
302     writefitinfo->Draw();
303     logparfunc->Draw("same");
304     logparfunc2->Draw("same");
305     logparfuncN->SetLineStyle(2);
306     logparfuncP->SetLineStyle(2);
307     logparfuncN->Draw("same");
308     logparfuncP->Draw("same");
309     CompleteSave(logpar_fit_can,"Bpred_ttbar_LogPar_Fit_To_OSOF");
310     delete logpar_fit_can;
311     return logparfunc2;
312     }
313    
314     vector<TF1*> do_new_extended_fit_to_plot(TH1F *prediction, TH1F *ossf, TH1F *osof) {
315     /* there are mainly two background contributions: Z+Jets (a) and ttbar (b). So:
316     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.
317     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.
318     Once we have these two components, we use the combined parameters to get the final function and we're done.
319     */
320     //Step 1: take the OSSF prediction - OSOF prediction, and fit a crystal ball function to it
321     TH1F *step1cb = (TH1F*)ossf->Clone("step1cb");
322     step1cb->Add(osof,-1);
323     vector<TF1*> functions = do_cb_fit_to_plot(step1cb,PlottingSetup::JZBPeakWidthData);
324     TF1 *zjetscrystalball = functions[1];
325    
326     //Step 2: use the OSOF distribution and look at the [10,100] GeV JZB range, and fit our log parabola
327     TH1F *ttbarprediction=(TH1F*)prediction->Clone("ttbarprediction");
328     ttbarprediction->Add(ossf,-1);//without the Z+Jets estimate, this is really just the ttbar estimate!
329     TF1 *ttbarlogpar = do_new_logpar_fit_to_plot(ttbarprediction);
330     TF1 *ttbarlogparP = new TF1("ttbarlogparP",LogParabolaP,0,(prediction->GetXaxis()->GetBinLowEdge(prediction->GetNbinsX())+prediction->GetXaxis()->GetBinWidth(prediction->GetNbinsX())),8);
331     TF1 *ttbarlogparN = new TF1("ttbarlogparN",LogParabolaN,0,(prediction->GetXaxis()->GetBinLowEdge(prediction->GetNbinsX())+prediction->GetXaxis()->GetBinWidth(prediction->GetNbinsX())),8);
332    
333     //and now fuse the two!
334     TF1 *kmlp = new TF1("kmlp", CrystalBallPlusLogParabola, 0,(prediction->GetXaxis()->GetBinLowEdge(prediction->GetNbinsX())+prediction->GetXaxis()->GetBinWidth(prediction->GetNbinsX())),8);
335     TF1 *kmlpP= new TF1("kmlpP",CrystalBallPlusLogParabolaP,0,(prediction->GetXaxis()->GetBinLowEdge(prediction->GetNbinsX())+prediction->GetXaxis()->GetBinWidth(prediction->GetNbinsX())),8);
336     TF1 *kmlpN= new TF1("kmlpN",CrystalBallPlusLogParabolaN,0,(prediction->GetXaxis()->GetBinLowEdge(prediction->GetNbinsX())+prediction->GetXaxis()->GetBinWidth(prediction->GetNbinsX())),8);
337     double kmlp_pars[10];
338     for(int i=0;i<5;i++) kmlp_pars[i]=zjetscrystalball->GetParameter(i);
339     for(int i=0;i<3;i++) kmlp_pars[5+i]=ttbarlogpar->GetParameter(i);
340     ttbarlogparP->SetParameters(ttbarlogpar->GetParameters());
341     ttbarlogparN->SetParameters(ttbarlogpar->GetParameters());
342     kmlp->SetParameters(kmlp_pars);
343     prediction->Fit(kmlp);
344    
345     kmlpP->SetParameters(kmlp->GetParameters());
346     kmlpN->SetParameters(kmlp->GetParameters());
347    
348     // now that we're done, let's save all of this so we can have a look at it afterwards.
349     TCanvas *can = new TCanvas("can","Prediction Fit Canvas");
350     can->SetLogy(1);
351     prediction->SetMarkerColor(kRed);
352     prediction->Draw();
353    
354     kmlp->SetLineColor(TColor::GetColor("#04B404"));
355     kmlpP->SetLineColor(TColor::GetColor("#04B404"));
356     kmlpN->SetLineColor(TColor::GetColor("#04B404"));
357     kmlp->Draw("same");
358     kmlpN->SetLineStyle(2);
359     kmlpP->SetLineStyle(2);
360     kmlpN->Draw("same");
361     kmlpP->Draw("same");
362    
363     ttbarlogpar->SetLineColor(TColor::GetColor("#CC2EFA"));//purple
364     ttbarlogpar->Draw("same");
365     ttbarlogparP->SetLineColor(TColor::GetColor("#CC2EFA"));//purple
366     ttbarlogparN->SetLineColor(TColor::GetColor("#CC2EFA"));//purple
367     ttbarlogparP->SetLineStyle(2);
368     ttbarlogparN->SetLineStyle(2);
369     ttbarlogparP->Draw("same");
370     ttbarlogparN->Draw("same");
371    
372     functions[0]->Draw("same");functions[1]->Draw("same");functions[2]->Draw("same");
373    
374     TLegend *analyticalBpredLEG = make_legend("",0.5,0.55);
375     analyticalBpredLEG->AddEntry(prediction,"predicted","p");
376     analyticalBpredLEG->AddEntry(functions[1],"Crystal Ball fit","l");
377     analyticalBpredLEG->AddEntry(functions[0],"1#sigma Crystal Ball fit","l");
378     analyticalBpredLEG->AddEntry(ttbarlogparN,"TTbar fit","l");
379     analyticalBpredLEG->AddEntry(ttbarlogpar,"1#sigma TTbar fit","l");
380     analyticalBpredLEG->AddEntry(kmlp,"Combined function","l");
381     analyticalBpredLEG->AddEntry(kmlpN,"1#sigma combined function","l");
382     analyticalBpredLEG->Draw("same");
383    
384     CompleteSave(can,"Bpred_MC_Analytical_Function_Composition");
385     delete can;
386    
387     //and finally: prep return functions
388     vector<TF1*> return_functions;
389     return_functions.push_back(kmlpN);
390     return_functions.push_back(kmlp);
391     return_functions.push_back(kmlpP);
392    
393     return_functions.push_back(ttbarlogparN);
394     return_functions.push_back(ttbarlogpar);
395     return_functions.push_back(ttbarlogparP);
396    
397     return_functions.push_back(functions[0]);
398     return_functions.push_back(functions[1]);
399     return_functions.push_back(functions[2]);
400    
401     return return_functions;
402     }
403     void do_new_prediction_plot(string jzb, TCanvas *globalcanvas, float sigma, float high, bool is_data, bool overlay_signal = false )
404     {
405     int nbins=100;
406     if(is_data) nbins=50;
407     float low=0;
408     float hi=500;
409    
410    
411     TH1F *RcorrJZBeemm = allsamples.Draw("RcorrJZBeemm",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,is_data, luminosity,allsamples.FindSample("TTJet"));
412     TH1F *LcorrJZBeemm = allsamples.Draw("LcorrJZBeemm",("-"+jzb).c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,is_data, luminosity,allsamples.FindSample("TTJet"));
413     TH1F *RcorrJZBem = allsamples.Draw("RcorrJZBem",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,is_data, luminosity,allsamples.FindSample("TTJet"));
414     TH1F *LcorrJZBem = allsamples.Draw("LcorrJZBem",("-"+jzb).c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,is_data, luminosity,allsamples.FindSample("TTJet"));
415    
416     TH1F *RcorrJZBSBem;
417     TH1F *LcorrJZBSBem;
418     TH1F *RcorrJZBSBeemm;
419     TH1F *LcorrJZBSBeemm;
420    
421     TH1F *lm4RcorrJZBeemm = allsamples.Draw("lm4RcorrJZBeemm",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,is_data, luminosity,allsamples.FindSample("LM4"));
422    
423 buchmann 1.3 if(RestrictToMassPeak||!PlottingSetup::UseSidebandsForcJZB) {
424 buchmann 1.1 RcorrJZBSBem = allsamples.Draw("RcorrJZBSBem",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets,is_data, luminosity,allsamples.FindSample("TTJet"));
425     LcorrJZBSBem = allsamples.Draw("LcorrJZBSBem",("-"+jzb).c_str(),nbins,low,hi, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets,is_data, luminosity,allsamples.FindSample("TTJet"));
426     RcorrJZBSBeemm = allsamples.Draw("RcorrJZBSBeemm",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets,is_data, luminosity,allsamples.FindSample("TTJet"));
427     LcorrJZBSBeemm = allsamples.Draw("LcorrJZBSBeemm",("-"+jzb).c_str(),nbins,low,hi, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets,is_data, luminosity,allsamples.FindSample("TTJet"));
428     }
429    
430     TH1F *conta = allsamples.Draw("conta",jzb.c_str(),2*nbins,-hi,hi, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,is_data, luminosity,allsamples.FindSample("TTJet"));
431    
432     TFile *ttbarfile=new TFile("ttbarfile.root","RECREATE");
433     conta->Write();
434     ttbarfile->Close();
435    
436     TH1F *Bpred = (TH1F*)LcorrJZBeemm->Clone("Bpred");
437     if(RestrictToMassPeak) {
438     Bpred->Add(RcorrJZBem,1.0/3.);
439     Bpred->Add(LcorrJZBem,-1.0/3.);
440     Bpred->Add(RcorrJZBSBem,1.0/3.);
441     Bpred->Add(LcorrJZBSBem,-1.0/3.);
442     Bpred->Add(RcorrJZBSBeemm,1.0/3.);
443     Bpred->Add(LcorrJZBSBeemm,-1.0/3.);
444     } else {
445     Bpred->Add(RcorrJZBem,1.0);
446     Bpred->Add(LcorrJZBem,-1.0);
447     }
448    
449     globalcanvas->cd();
450     globalcanvas->SetLogy(1);
451    
452     RcorrJZBeemm->SetMarkerStyle(20);
453     RcorrJZBeemm->GetXaxis()->SetRangeUser(0,high);
454     RcorrJZBeemm->Draw("e1x0");
455    
456     Bpred->SetLineColor(kRed);
457     Bpred->SetStats(0);
458     Bpred->Draw("hist,same");
459    
460     if ( overlay_signal ) {
461     lm4RcorrJZBeemm->SetLineColor(TColor::GetColor("#088A08"));
462     lm4RcorrJZBeemm->Draw("histo,same");
463     }
464     RcorrJZBeemm->SetMarkerSize(DataMarkerSize);
465    
466     TLegend *legBpred = make_legend("",0.6,0.55);
467     if(is_data)
468     {
469     vector<TF1*> analytical_function = do_new_extended_fit_to_plot(Bpred,LcorrJZBeemm,LcorrJZBem);
470     globalcanvas->cd();
471     analytical_function[0]->Draw("same"); analytical_function[1]->Draw("same");analytical_function[2]->Draw("same");
472     legBpred->AddEntry(RcorrJZBeemm,"observed","p");
473     legBpred->AddEntry(Bpred,"predicted","l");
474     legBpred->AddEntry(analytical_function[1],"predicted fit","l");
475     legBpred->AddEntry(analytical_function[2],"stat. uncert.","l");
476     if ( overlay_signal ) legBpred->AddEntry(lm4RcorrJZBeemm,"LM4","l");
477     legBpred->Draw();
478     CompleteSave(globalcanvas,"Bpred_Data");
479     }
480     else {
481     vector<TF1*> analytical_function = do_new_extended_fit_to_plot(Bpred,LcorrJZBeemm,LcorrJZBem);
482     globalcanvas->cd();
483     analytical_function[0]->Draw("same"); analytical_function[1]->Draw("same");analytical_function[2]->Draw("same");
484     legBpred->AddEntry(RcorrJZBeemm,"MC observed","p");
485     legBpred->AddEntry(Bpred,"MC predicted","l");
486     legBpred->AddEntry(analytical_function[1],"predicted fit","l");
487     legBpred->AddEntry(analytical_function[2],"stat. uncert.","l");
488     if ( overlay_signal ) legBpred->AddEntry(lm4RcorrJZBeemm,"LM4","l");
489     legBpred->Draw();
490     CompleteSave(globalcanvas,"Bpred_MC");
491     }
492    
493     TH1F *Bpredem = (TH1F*)LcorrJZBeemm->Clone("Bpredem");
494     Bpredem->Add(RcorrJZBem);
495     Bpredem->Add(LcorrJZBem,-1);
496     TH1F *BpredSBem = (TH1F*)LcorrJZBeemm->Clone("BpredSBem");
497     BpredSBem->Add(RcorrJZBSBem);
498     Bpred->Add(LcorrJZBSBem,-1);
499     TH1F *BpredSBeemm = (TH1F*)LcorrJZBeemm->Clone("BpredSBeemm");
500     BpredSBeemm->Add(RcorrJZBSBeemm);
501     BpredSBeemm->Add(LcorrJZBSBeemm,-1.0);
502     globalcanvas->cd();
503     globalcanvas->SetLogy(1);
504    
505     RcorrJZBeemm->SetMarkerStyle(20);
506     RcorrJZBeemm->GetXaxis()->SetRangeUser(0,high);
507     RcorrJZBeemm->Draw("e1x0");
508     RcorrJZBeemm->SetMarkerSize(DataMarkerSize);
509    
510     Bpredem->SetLineColor(kRed+1);
511     Bpredem->SetStats(0);
512     Bpredem->Draw("hist,same");
513    
514     BpredSBem->SetLineColor(kGreen+2);//TColor::GetColor("#0B6138"));
515     BpredSBem->SetLineStyle(2);
516     BpredSBem->Draw("hist,same");
517    
518     BpredSBeemm->SetLineColor(kBlue+1);
519     BpredSBeemm->SetLineStyle(3);
520     BpredSBeemm->Draw("hist,same");
521    
522     TLegend *legBpredc = make_legend("",0.6,0.55);
523     if(is_data)
524     {
525     legBpredc->AddEntry(RcorrJZBeemm,"observed","p");
526     legBpredc->AddEntry(Bpredem,"OFZP","l");
527     legBpredc->AddEntry(BpredSBem,"OFSB","l");
528     legBpredc->AddEntry(BpredSBeemm,"SFSB","l");
529     legBpredc->Draw();
530     CompleteSave(globalcanvas,"Bpred_Data_comparison");
531     }
532     else {
533     legBpredc->AddEntry(RcorrJZBeemm,"MC observed","p");
534     legBpredc->AddEntry(Bpredem,"MC OFZP","l");
535     legBpredc->AddEntry(BpredSBem,"MC OFSB","l");
536     legBpredc->AddEntry(BpredSBeemm,"MC SFSB","l");
537     legBpredc->Draw();
538     legBpredc->Draw();
539     CompleteSave(globalcanvas,"Bpred_MC_comparison_ttbar");
540     }
541     delete RcorrJZBeemm;
542     delete LcorrJZBeemm;
543     delete RcorrJZBem;
544     delete LcorrJZBem;
545     delete RcorrJZBSBem;
546     delete LcorrJZBSBem;
547     delete RcorrJZBSBeemm;
548     delete LcorrJZBSBeemm;
549     delete lm4RcorrJZBeemm;
550     }
551    
552     void do_new_prediction_plots(string mcjzb, string datajzb, float DataSigma, float MCSigma, bool overlay_signal ) {
553     TCanvas *globalcanvas = new TCanvas("globalcanvas","Prediction Canvas");
554     // do_new_prediction_plot(datajzb,globalcanvas,DataSigma,jzbHigh ,data,overlay_signal);
555     do_new_prediction_plot(mcjzb,globalcanvas,MCSigma,jzbHigh ,mc,overlay_signal);
556     }
557