ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/MetPlotting.C
Revision: 1.26
Committed: Thu Sep 13 07:39:57 2012 UTC (12 years, 7 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.25: +19 -22 lines
Log Message:
Cleaned up code (less verbose now)

File Contents

# User Rev Content
1 buchmann 1.1 #include <iostream>
2    
3     using namespace std;
4    
5 buchmann 1.17 float Get_Met_Z_Prediction(TCut JetCut, float MetCut, int isdata, bool isDYonly);
6    
7 buchmann 1.24 namespace MetPlotsSpace {
8     float Zprediction_Uncertainty=0.3;
9     }
10 buchmann 1.17
11 buchmann 1.1 TGraphErrors* MakeErrorGraph(TH1F *histo) {
12    
13     float dx[histo->GetNbinsX()];
14     float dy[histo->GetNbinsX()];
15     float x[histo->GetNbinsX()];
16     float y[histo->GetNbinsX()];
17     for(int i=1;i<=histo->GetNbinsX();i++) {
18     x[i-1]=histo->GetBinCenter(i);
19     y[i-1]=histo->GetBinContent(i);
20     if(i>1) dx[i-1]=(histo->GetBinCenter(i)-histo->GetBinCenter(i-1))/2.0;
21     else dx[i-1]=(histo->GetBinCenter(i+1)-histo->GetBinCenter(i))/2.0;
22     dy[i-1]=histo->GetBinError(i);
23     }
24    
25     TGraphErrors *gr = new TGraphErrors(histo->GetNbinsX(),x,y,dx,dy);
26     gr->SetFillColor(TColor::GetColor("#2E9AFE"));
27     return gr;
28     }
29    
30 fronga 1.22 void ProduceMetPlotsWithCut(TCut cut, string name, float cutat, int njets, bool doMC = false, float ymax = 80 ) {
31 fronga 1.16
32 buchmann 1.1 TCanvas *tcan = new TCanvas("tcan","tcan");
33     cout << "Doing met plots" << endl;
34 buchmann 1.2 stringstream MetBaseCuts;
35 fronga 1.7 MetBaseCuts << "met[4]>" << cutat << "&&" << cut.GetTitle();
36     stringstream snjets;
37     snjets << njets;
38 buchmann 1.2 TCut MetBaseCut(MetBaseCuts.str().c_str());
39 fronga 1.10 TCut nJetsSignal(PlottingSetup::basicqualitycut&&("pfJetGoodNum40>="+snjets.str()).c_str());
40 fronga 1.22 TCut nJetsControl(PlottingSetup::basiccut&&("met[4]>100&&met[4]<150&&pfJetGoodID[0]!=0&&pfJetGoodNum40=="+snjets.str()+"-1").c_str()); // Only njets vary
41 fronga 1.16
42     if ( !PlottingSetup::openBox ) {
43     nJetsSignal += TCut("mll>120");
44     }
45 buchmann 1.1
46     //compute SF / OF rate in (CR1+CR2), should give 0.941 +/- 0.05
47 fronga 1.7
48     // Create histograms
49 fronga 1.9 //int nbins = 30;
50     int nbins = 60;
51 buchmann 1.17 float xmin=15., xmax = 315.;
52     TH1F *mllsigEE = allsamples.Draw("mllsigEE","mll",nbins,xmin,xmax,"m_{ee} [GeV]", "events",TCut(cutOSSF&&MetBaseCut&&nJetsSignal&&"id1==0"),data,PlottingSetup::luminosity);
53     TH1F *mllsigMM = allsamples.Draw("mllsigMM","mll",nbins,xmin,xmax,"m_{#mu#mu} [GeV]","events",TCut(cutOSSF&&MetBaseCut&&nJetsSignal&&"id1==1"),data,PlottingSetup::luminosity);
54 fronga 1.22 TH1F *mllscon = allsamples.Draw("mllscon","mll",nbins,xmin,xmax,"m_{ll} [GeV]", "events",TCut(cutOSSF&&cut&&nJetsControl),data,PlottingSetup::luminosity);
55 fronga 1.7 TH1F *mllOsig = allsamples.Draw("mllOsig", "mll",nbins,xmin,xmax,"m_{ll} [GeV]","events",TCut(cutOSOF&&MetBaseCut&&nJetsSignal),data,PlottingSetup::luminosity);
56 fronga 1.22 TH1F *mllOscon = allsamples.Draw("mllOscon","mll",nbins,xmin,xmax,"m_{ll} [GeV]","events",TCut(cutOSOF&&cut&&nJetsControl),data,PlottingSetup::luminosity);
57 fronga 1.7
58 fronga 1.8 TH1F* mllsig = (TH1F*)mllsigEE->Clone("mllsig");
59     mllsig->Add(mllsigMM);
60     mllsig->GetXaxis()->SetTitle("m_{ll} [GeV]");
61    
62     THStack *mcMllsig, *mcMllsigEE,*mcMllsigMM,*mcMllscon,*mcMllsconEE,*mcMllsconMM, *mcMllOsig, *mcMllOscon;
63 fronga 1.7 if ( doMC ) {
64     name += "_mc";
65 fronga 1.8 mcMllsig = new THStack(allsamples.DrawStack("mcMllsig","mll",nbins,xmin,xmax,"m_{ll} [GeV]","events",TCut(cutOSSF&&MetBaseCut&&nJetsSignal),mc,PlottingSetup::luminosity));
66     mcMllsigEE = new THStack(allsamples.DrawStack("mcMllsigEE","mll",nbins,xmin,xmax,"m_{ee} [GeV]","events",TCut(cutOSSF&&MetBaseCut&&nJetsSignal&&"id1==0"),mc,PlottingSetup::luminosity));
67     mcMllsigMM = new THStack(allsamples.DrawStack("mcMllsigMM","mll",nbins,xmin,xmax,"m_{#mu#mu} [GeV]","events",TCut(cutOSSF&&MetBaseCut&&nJetsSignal&&"id1==1"),mc,PlottingSetup::luminosity));
68 fronga 1.22 mcMllscon = new THStack(allsamples.DrawStack("mcMllscon","mll",nbins,xmin,xmax,"m_{ll} [GeV]","events",TCut(cutOSSF&&cut&&nJetsControl),mc,PlottingSetup::luminosity));
69 fronga 1.8 mcMllOsig = new THStack(allsamples.DrawStack("mcMllOsig","mll",nbins,xmin,xmax,"m_{ll} [GeV]","events",TCut(cutOSOF&&MetBaseCut&&nJetsSignal),mc,PlottingSetup::luminosity));
70 fronga 1.22 mcMllOscon= new THStack(allsamples.DrawStack("mcMllOscon","mll",nbins,xmin,xmax,"m_{ll} [GeV]","events",TCut(cutOSOF&&cut&&nJetsControl),mc,PlottingSetup::luminosity));
71 fronga 1.8
72 fronga 1.7 }
73 fronga 1.8
74 buchmann 1.1 mllOsig->SetLineColor(kRed);
75     mllOscon->SetLineColor(kRed);
76    
77 fronga 1.7 TH1F *zlineshape = allsamples.Draw("zlineshape","mll",nbins,xmin,xmax,"m_{ll} (GeV)","events",cutOSSF&&TCut("pfJetGoodNum40==1")&&cut,data,PlottingSetup::luminosity);
78 buchmann 1.13 TH1F *zlineshapeControl = (TH1F*)zlineshape->Clone("zlineshapeControl");
79 buchmann 1.21 // TH1F *zlineshapeFINE = allsamples.Draw("zlineshapeFINE","mll",50*nbins,xmin,xmax,"m_{ll} (GeV)","events",cutOSSF&&TCut("pfJetGoodNum40==1")&&cut,data,PlottingSetup::luminosity);
80     //
81     // float scalefactor = Get_Met_Z_Prediction(nJetsSignal,cutat, data, false) / (zlineshapeFINE->Integral(zlineshapeFINE->FindBin(91.1-20),zlineshapeFINE->FindBin(91.1+20)));
82     // float scalefactor_Control = Get_Met_Z_Prediction(nJetsControl,cutat, data, false) / (zlineshapeFINE->Integral(zlineshapeFINE->FindBin(91.1-20),zlineshapeFINE->FindBin(91.1+20)));
83     // delete zlineshapeFINE;
84    
85    
86     Int_t scaleBinLow = mllsig->FindBin(86);
87     Int_t scaleBinHigh = mllsig->FindBin(94);
88     float scalefactor = (mllsig->Integral(scaleBinLow,scaleBinHigh)-mllOsig->Integral(scaleBinLow,scaleBinHigh));
89     scalefactor /= zlineshape->Integral(scaleBinLow,scaleBinHigh);
90    
91     float scalefactor_Control = (mllscon->Integral(scaleBinLow,scaleBinHigh)-mllOscon->Integral(scaleBinLow,scaleBinHigh));
92     scalefactor_Control /= zlineshapeControl->Integral(scaleBinLow,scaleBinHigh);
93    
94     cout << "Bins for scaling : " << scaleBinLow << " : " << scaleBinHigh << endl;
95     cout << "Scale factors : " << scalefactor << " : " << scalefactor_Control << endl;
96 buchmann 1.17
97 buchmann 1.3 zlineshape->Scale(scalefactor);
98     zlineshape->SetLineColor(kBlue);
99     zlineshape->SetLineStyle(2);
100    
101 buchmann 1.13 zlineshapeControl->Scale(scalefactor_Control);
102     zlineshapeControl->SetLineColor(kBlue);
103     zlineshapeControl->SetLineStyle(2);
104    
105 buchmann 1.3 TH1F *subtracted = (TH1F*)mllsig->Clone("subtracted");
106 buchmann 1.13 TH1F *subtractedControl = (TH1F*)mllscon->Clone("subtractedControl");
107 buchmann 1.3 TH1F *baseline = (TH1F*)mllOsig->Clone("baseline");
108 buchmann 1.13 TH1F *baselineControl = (TH1F*)mllOscon->Clone("baselineControl");
109 buchmann 1.3 for(int i=1;i<=subtracted->GetNbinsX();i++) {
110     subtracted->SetBinContent(i,mllsig->GetBinContent(i)-mllOsig->GetBinContent(i));
111 buchmann 1.13 subtractedControl->SetBinContent(i,mllscon->GetBinContent(i)-mllOscon->GetBinContent(i));
112 buchmann 1.3 baseline->SetBinContent(i,0);
113 buchmann 1.13 baselineControl->SetBinContent(i,0);
114 buchmann 1.3 }
115    
116     TH1F *prediction = (TH1F*)mllOsig->Clone("prediction");
117     prediction->Add(zlineshape);
118     prediction->SetLineColor(TColor::GetColor("#CF35CA"));
119 fronga 1.7
120     TH1F *control_prediction = (TH1F*)mllOscon->Clone("control_prediction");
121     control_prediction->SetLineColor(TColor::GetColor("#CF35CA"));
122    
123     // FIX Y RANGE TO EASE COMPARISON
124 fronga 1.22 mllsig->SetMaximum(ymax);
125     mllsigEE->SetMaximum(ymax);
126     mllsigMM->SetMaximum(ymax);
127     mllOsig->SetMaximum(ymax);
128     mllOscon->SetMaximum(ymax);
129     subtracted->SetMaximum(60);
130     subtracted->SetMinimum(-30);
131     subtractedControl->SetMaximum(65);
132     subtractedControl->SetMinimum(-30);
133 fronga 1.7
134 buchmann 1.3
135 fronga 1.7 // 1.- Signal region comparison
136 fronga 1.14 TBox *srbox = new TBox(xmin,0,70,mllsig->GetMaximum());
137 buchmann 1.1 srbox->SetFillStyle(0);
138     srbox->SetLineColor(TColor::GetColor("#298A08"));
139     srbox->SetLineWidth(3);
140 buchmann 1.23
141 fronga 1.10
142 buchmann 1.2 stringstream MetHeader;
143 fronga 1.19 MetHeader << "N_{j}#geq" << snjets.str() << ", MET>" << cutat << " GeV";
144     stringstream MetHeaderCon;
145 fronga 1.22 MetHeaderCon << "N_{j}=" << njets-1 << ", 100<MET<150 GeV";
146 fronga 1.10 stringstream saveasSig;
147     saveasSig << "MetPlots/mll_sig" << cutat << "__" << name;
148    
149     TLegend* leg;
150     if ( !doMC ) {
151     tcan->cd();
152     //mllsig->GetYaxis()->SetRangeUser(0,mllsig->GetMaximum()*1.2);
153     mllsig->Draw();
154     TGraphErrors *stat3j = MakeErrorGraph(prediction);
155     stat3j->Draw("2,same");
156     mllOsig->Draw("histo,same");
157     zlineshape->Draw("histo,same");
158     prediction->Draw("histo,same");
159     mllsig->Draw("same");
160     DrawPrelim();
161     leg = make_legend();
162 fronga 1.19 leg->SetX1(0.52);
163     leg->SetHeader(MetHeader.str().c_str());
164 fronga 1.10 leg->AddEntry(mllsig,"Data","PL");
165     leg->AddEntry(prediction,"All bg prediction","L");
166 fronga 1.7 leg->AddEntry(mllOsig,"bg without Z","L");
167     leg->AddEntry(zlineshape,"Z lineshape","L");
168     leg->AddEntry(stat3j,"stat. uncert.","F");
169 fronga 1.10 leg->AddEntry(srbox,"SR","F");
170     leg->Draw();
171     srbox->Draw();
172     CompleteSave(tcan,saveasSig.str());
173     } else {
174 fronga 1.19
175 fronga 1.10 TPad* rcan = new TPad("rcan","rcan",0,0,1,1);
176     rcan->cd();
177     mllsig->Draw();
178     mcMllsig->Draw("same");
179     prediction->Draw("histo,same");
180     mllsig->Draw("same");
181     DrawPrelim();
182     leg = allsamples.allbglegend();
183 fronga 1.19 leg->SetHeader(MetHeader.str().c_str());
184     leg->SetX1(0.52);
185 fronga 1.10 leg->AddEntry(prediction,"All bg prediction","L");
186     leg->AddEntry(srbox,"SR","F");
187     leg->Draw();
188     srbox->Draw();
189     save_with_ratio( mllsig, *mcMllsig, rcan, saveasSig.str() );
190 fronga 1.7 }
191 buchmann 1.1
192    
193 fronga 1.8 // 1b. MC: split ee and mumu
194     if ( doMC ) {
195 fronga 1.10 TPad* rcan = new TPad("rcan","rcan",0,0,1,1);
196     rcan->cd();
197 fronga 1.8 mllsigEE->Draw();
198     mcMllsigEE->Draw("same");
199     mllsigEE->Draw("same");
200     DrawPrelim();
201     leg->Draw();
202     srbox->Draw();
203 fronga 1.10 save_with_ratio( mllsigEE, *mcMllsigEE,rcan->cd(),saveasSig.str()+"_ee" );
204 fronga 1.8
205 fronga 1.10 rcan = new TPad("rcan","rcan",0,0,1,1);
206     rcan->cd();
207 fronga 1.8 mllsigMM->Draw();
208     mcMllsigMM->Draw("same");
209     mllsigMM->Draw("same");
210     DrawPrelim();
211     leg->Draw();
212     srbox->Draw();
213 fronga 1.10 save_with_ratio( mllsigMM, *mcMllsigMM,rcan,saveasSig.str()+"_mm" );
214 fronga 1.8 }
215 fronga 1.14
216     // 1c. MC: compare of and sf
217     if ( doMC ) {
218     TH1F* hMcMllsig = CollapseStack( *mcMllsig);
219 fronga 1.15 leg = allsamples.allbglegend("");
220 fronga 1.19 leg->SetHeader(MetHeader.str().c_str());
221 fronga 1.15 // Change "Data" label by hand
222     ((TLegendEntry*)leg->GetListOfPrimitives()->At(0))->SetLabel("Same-flavor (MC)");
223 fronga 1.14 TPad* rcan = new TPad("rcan","rcan",0,0,1,1);
224     rcan->cd();
225 fronga 1.22 hMcMllsig->SetMaximum(ymax);
226 fronga 1.14 hMcMllsig->Draw("E");
227     mcMllOsig->Draw("same,hist");
228     hMcMllsig->Draw("same,E");
229     DrawMCPrelim();
230 fronga 1.19 leg->SetX1(0.52);
231 fronga 1.14 leg->AddEntry(srbox,"SR","F");
232     leg->Draw();
233     srbox->Draw();
234     save_with_ratio( hMcMllsig, *mcMllOsig, rcan, saveasSig.str()+"_mconly");
235    
236     }
237 buchmann 1.1
238 fronga 1.7 // 2.- Signal region comparison - LOG scale
239 fronga 1.10 if ( !doMC ) {
240     tcan->cd();
241     mllsig->SetMinimum(0.2); // FIX Y RANGE TO EASE COMPARISON
242     //mllsig->SetMaximum(mllsig->GetMaximum()*4.0);
243     srbox->SetY2(mllsig->GetMaximum());
244    
245     tcan->SetLogy(1);
246     stringstream saveasSig2;
247     saveasSig2 << "MetPlots/mll_sig_ZLINESHAPE_" << cutat << "__" << name;
248    
249     CompleteSave(tcan,saveasSig2.str());
250     tcan->SetLogy(0);
251     }
252 buchmann 1.1
253 fronga 1.7
254     // 3.- Signal region, background subtracted
255 fronga 1.8 if ( !doMC ) {
256 fronga 1.10 tcan->cd();
257 buchmann 1.13 for(int i=1;i<=subtracted->GetNbinsX();i++) {
258 fronga 1.8 subtracted->SetBinContent(i,subtracted->GetBinContent(i)-zlineshape->GetBinContent(i));
259 buchmann 1.13 subtractedControl->SetBinContent(i,subtractedControl->GetBinContent(i)-zlineshapeControl->GetBinContent(i));
260 fronga 1.8 }
261 buchmann 1.3
262 fronga 1.8 TGraphErrors *subtrerr = MakeErrorGraph(baseline);
263     subtracted->Draw();
264     subtrerr->Draw("2,same");
265     subtracted->Draw("same");
266     DrawPrelim();
267     TLegend *DiffLeg = make_legend();
268 fronga 1.22 DiffLeg->SetX1(0.4);
269 fronga 1.8 DiffLeg->SetFillStyle(0);
270 fronga 1.19 DiffLeg->SetHeader(MetHeader.str().c_str());
271 fronga 1.8 DiffLeg->AddEntry(subtracted,"observed - predicted","PL");
272 fronga 1.19 DiffLeg->AddEntry(subtrerr,"stat. uncert","F");
273 fronga 1.8 DiffLeg->AddEntry((TObject*)0,"","");
274     DiffLeg->AddEntry((TObject*)0,"","");
275     DiffLeg->Draw();
276    
277     stringstream saveasSigSub;
278     saveasSigSub << "MetPlots/mll_sig_SUBTRACTED_" << cutat << "__" << name;
279    
280 fronga 1.22 //CompleteSave(tcan,saveasSigSub.str());
281 buchmann 1.13
282     // 3a.- Control region, background subtracted
283     TGraphErrors *subtrerrControl = MakeErrorGraph(baselineControl);
284     subtractedControl->Draw();
285     subtrerrControl->Draw("2,same");
286     subtractedControl->Draw("same");
287     DrawPrelim();
288 fronga 1.19 DiffLeg->SetHeader(MetHeaderCon.str().c_str());
289 buchmann 1.13 DiffLeg->Draw();
290     saveasSigSub.str("");
291 fronga 1.16 saveasSigSub << "MetPlots/mll_con_SUBTRACTED_" << cutat << "__" << name;
292 fronga 1.22 //CompleteSave(tcan,saveasSigSub.str());
293 buchmann 1.13
294    
295    
296 fronga 1.8 // 4.- Signal region, background subtracted, errors added in quadrature
297     TGraphErrors *subtrerr2 = (TGraphErrors*)subtrerr->Clone("subtrerr2");
298 buchmann 1.13 for(int i=1;i<=subtrerr2->GetN();i++) {
299     subtrerr2->SetPoint(i-1,subtracted->GetBinCenter(i),subtracted->GetBinContent(i));
300     subtrerr2->SetPointError(i-1,subtrerr2->GetErrorX(i),TMath::Sqrt(subtrerr2->GetErrorY(i)*subtrerr2->GetErrorY(i)+subtracted->GetBinError(i)*subtracted->GetBinError(i)));
301 fronga 1.8 }
302     TLine* l = new TLine(subtracted->GetBinLowEdge(1),0.,subtracted->GetBinLowEdge(subtracted->GetNbinsX()-1)+subtracted->GetBinWidth(1),0.);
303     l->SetLineWidth(subtracted->GetLineWidth());
304     subtracted->Draw();
305     subtrerr2->Draw("2,same");
306     l->Draw("same");
307     subtracted->Draw("same");
308     DrawPrelim();
309     TLegend *DiffLeg2 = make_legend();
310 fronga 1.22 DiffLeg2->SetX1(0.4);
311 fronga 1.19 DiffLeg2->SetHeader(MetHeader.str().c_str());
312 fronga 1.8 DiffLeg2->SetFillStyle(0);
313     DiffLeg2->AddEntry(subtracted,"observed - predicted","PL");
314 fronga 1.19 DiffLeg2->AddEntry(subtrerr2,"stat. uncert","F");
315 fronga 1.8 DiffLeg2->AddEntry((TObject*)0,"","");
316     DiffLeg2->AddEntry((TObject*)0,"","");
317     DiffLeg2->Draw();
318    
319     stringstream saveasSigSub2;
320     saveasSigSub2 << "MetPlots/mll_sig_SUBTRACTED_quadr_" << cutat << "__" << name;
321 fronga 1.5
322 fronga 1.8 CompleteSave(tcan,saveasSigSub2.str());
323 fronga 1.5
324 buchmann 1.13
325    
326     //4a.- Control region, background subtracted, errors added in quadrature
327     TGraphErrors *subtrerr2Control = (TGraphErrors*)subtrerrControl->Clone("subtrerr2Control");
328     for(int i=1;i<=subtrerr2Control->GetN();i++) {
329     subtrerr2Control->SetPoint(i-1,subtractedControl->GetBinCenter(i),subtractedControl->GetBinContent(i));
330     float width=subtrerr2Control->GetErrorX(i);
331     if(i==subtrerr2Control->GetN()) width=subtrerr2Control->GetErrorX(i-1);
332     subtrerr2Control->SetPointError(i-1,width,TMath::Sqrt(subtrerr2Control->GetErrorY(i)*subtrerr2Control->GetErrorY(i)+subtractedControl->GetBinError(i)*subtractedControl->GetBinError(i)));
333     }
334     TLine* lControl = new TLine(subtractedControl->GetBinLowEdge(1),0.,subtractedControl->GetBinLowEdge(subtractedControl->GetNbinsX()-1)+subtractedControl->GetBinWidth(1),0.);
335     lControl->SetLineWidth(subtractedControl->GetLineWidth());
336     subtractedControl->Draw();
337     subtrerr2Control->Draw("2,same");
338     lControl->Draw("same");
339     subtractedControl->Draw("same");
340     DrawPrelim();
341 fronga 1.22 DiffLeg2->SetHeader(MetHeaderCon.str().c_str());
342 buchmann 1.13 DiffLeg2->Draw();
343    
344     saveasSigSub2.str("");
345 fronga 1.16 saveasSigSub2 << "MetPlots/mll_con_SUBTRACTED_quadr_" << cutat << "__" << name;
346 buchmann 1.13
347     CompleteSave(tcan,saveasSigSub2.str());
348    
349 fronga 1.8 delete DiffLeg;
350     delete DiffLeg2;
351 buchmann 1.13
352 fronga 1.8 } // !doMC
353 buchmann 1.3
354    
355 fronga 1.7 // 5.- Control region comparison
356 fronga 1.16 // scalefactor = (mllscon->Integral(scaleBinLow,scaleBinHigh)-mllOscon->Integral(scaleBinLow,scaleBinHigh));
357     // scalefactor /= zlineshape->Integral(scaleBinLow,scaleBinHigh);
358     // zlineshape->Scale(scalefactor);
359     control_prediction->Add(zlineshapeControl);
360    
361 fronga 1.22 control_prediction->SetMaximum(ymax); // FIX MAXIMUM TO EASE COMPARISON
362 fronga 1.7
363 fronga 1.14 TBox *cr1box = new TBox(xmin,0,70,control_prediction->GetMaximum());
364 buchmann 1.1 cr1box->SetFillStyle(0);
365     cr1box->SetLineColor(TColor::GetColor("#0404B4"));
366     cr1box->SetLineWidth(3);
367    
368 fronga 1.14 TBox *cr2box = new TBox(120,0,xmax,control_prediction->GetMaximum());
369 buchmann 1.1 cr2box->SetFillStyle(0);
370     cr2box->SetLineColor(TColor::GetColor("#0404B4"));
371     cr2box->SetLineWidth(3);
372     cr2box->SetLineStyle(2);
373    
374 fronga 1.10 stringstream saveasCon;
375     saveasCon << "MetPlots/mll_con" << cutat << "__" << name;
376    
377 fronga 1.7 TLegend *legc;
378 fronga 1.10 //control_prediction->GetYaxis()->SetRangeUser(0,control_prediction->GetMaximum()*1.3);
379     if ( !doMC ) {
380     tcan->cd();
381     control_prediction->Draw("hist");
382     TGraphErrors *stat2j = MakeErrorGraph(control_prediction);
383     stat2j->Draw("2,same");
384     mllOscon->Draw("same,hist");
385 fronga 1.16 zlineshapeControl->Draw("histo,same");
386 fronga 1.10 control_prediction->Draw("histo,same");
387     mllscon->Draw("same");
388     DrawPrelim();
389 fronga 1.7 legc = make_legend();
390 fronga 1.19 legc->SetX1(0.52);
391     legc->SetHeader(MetHeaderCon.str().c_str());
392 fronga 1.10 legc->AddEntry(mllscon,"Data","PL");
393     legc->AddEntry(control_prediction,"All bg","L");
394 fronga 1.16 legc->AddEntry(zlineshapeControl,"Z lineshape","L");
395 fronga 1.7 legc->AddEntry(mllOscon,"bg without Z","L");
396     legc->AddEntry(stat2j,"stat. uncert.","F");
397 fronga 1.10 legc->AddEntry(cr1box,"CR1","F");
398     legc->AddEntry(cr2box,"CR2","F");
399     legc->Draw();
400     cr1box->Draw();
401     cr2box->Draw();
402     CompleteSave(tcan,saveasCon.str());
403     } else {
404     TPad* rcan = new TPad("rcan","rcan",0,0,1,1);
405     rcan->cd();
406     control_prediction->Draw("hist");
407     mcMllscon->Draw("same");
408     control_prediction->Draw("histo,same");
409     mllscon->Draw("same");
410     DrawPrelim();
411     legc = allsamples.allbglegend();
412 fronga 1.19 legc->SetX1(0.52);
413     legc->SetHeader(MetHeaderCon.str().c_str());
414 fronga 1.10 legc->AddEntry(control_prediction,"All bg","L");
415     legc->AddEntry(cr1box,"CR1","F");
416     legc->AddEntry(cr2box,"CR2","F");
417     legc->Draw();
418     cr1box->Draw();
419     cr2box->Draw();
420     save_with_ratio( mllscon, *mcMllscon, rcan, saveasCon.str());
421     }
422 buchmann 1.1
423 fronga 1.7 // 6. - Opposite-flavour data/MC comparison
424     if ( doMC ) {
425 fronga 1.10 TPad* rcan = new TPad("rcan","rcan",0,0,1,1);
426     rcan->cd();
427 fronga 1.7 mllOsig->SetLineColor(kBlack);
428     mllOsig->Draw();
429     mcMllOsig->Draw("same");
430     mllOsig->Draw("same");
431     TLegend *legsdm = allsamples.allbglegend();
432 fronga 1.19 legsdm->SetHeader((MetHeader.str()+", OF").c_str());
433     legsdm->SetX1(0.52);
434 fronga 1.7 legsdm->Draw();
435     stringstream saveasSigOF;
436     saveasSigOF << "MetPlots/mll_sig_of_" << cutat << "__" << name;
437 fronga 1.10 save_with_ratio( mllOsig, *mcMllOsig, rcan, saveasSigOF.str());
438 fronga 1.7
439 fronga 1.10 rcan = new TPad("rcan","rcan",0,0,1,1);
440     rcan->cd();
441 fronga 1.7 mllOscon->SetLineColor(kBlack);
442     mllOscon->Draw();
443     mcMllOscon->Draw("same");
444     mllOscon->Draw("same");
445     TLegend *legcdm = allsamples.allbglegend();
446 fronga 1.19 legcdm->SetHeader((MetHeaderCon.str()+", OF").c_str());
447     legcdm->SetX1(0.52);
448 fronga 1.7 legcdm->Draw();
449     stringstream saveasConOF;
450     saveasConOF << "MetPlots/mll_con_of_" << cutat << "__" << name;
451 fronga 1.10 save_with_ratio( mllOscon, *mcMllOscon, rcan, saveasConOF.str());
452    
453 fronga 1.7 delete legsdm;
454     delete legcdm;
455 fronga 1.10 }
456 fronga 1.7
457 fronga 1.10 // Memory clean-up
458     if (doMC) {
459 fronga 1.7 delete mcMllscon;
460     delete mcMllOscon;
461     delete mcMllsig;
462 fronga 1.8 delete mcMllsigEE;
463     delete mcMllsigMM;
464 fronga 1.7 delete mcMllOsig;
465     }
466 buchmann 1.1
467     delete cr1box;
468     delete cr2box;
469     delete srbox;
470     delete legc;
471     delete leg;
472 fronga 1.7
473 buchmann 1.1 delete mllscon;
474     delete mllOscon;
475     delete mllsig;
476 fronga 1.8 delete mllsigEE;
477     delete mllsigMM;
478 buchmann 1.1 delete mllOsig;
479 fronga 1.7 delete zlineshape;
480 fronga 1.16 delete zlineshapeControl;
481 buchmann 1.6 delete tcan;
482 buchmann 1.1 }
483    
484 fronga 1.14 void DoMetPlots(string datajzb, string mcjzb) {
485 fronga 1.7 float metCuts[] = { 100., 150. };
486 fronga 1.22 float ymax[] = { 90., 80. };
487 fronga 1.7 int jetCuts[] = { 3, 2 };
488 fronga 1.11 //int jetCuts[] = { 3, 3 };
489 fronga 1.22 string leptCuts[] = { "pt1>20&&pt2>20", "pt1>20&&pt2>10" };
490 fronga 1.7 bool nomc(0),domc(1);
491     for ( int i=0; i<2; ++i ) {
492 fronga 1.22 ProduceMetPlotsWithCut(TCut(("mll>15&&"+leptCuts[i]).c_str()),"",metCuts[i],jetCuts[i],nomc,ymax[i]);
493     ProduceMetPlotsWithCut(TCut(("mll>15&&"+leptCuts[i]).c_str()),"",metCuts[i],jetCuts[i],domc,ymax[i]);
494     continue;
495 fronga 1.19 // stringstream jzbcut;
496     // jzbcut << "((is_data&&("<<datajzb<<")>0)||(!is_data&&("<<mcjzb<<")>0))";
497     // ProduceMetPlotsWithCut(TCut((jzbcut.str()+"&&mll>20&&"+leptCuts[i]).c_str()),"JZBpos",metCuts[i], jetCuts[i],domc);
498 fronga 1.11 //continue;
499 fronga 1.7 if (!PlottingSetup::is53reco) { // Old 5_2
500 fronga 1.22 ProduceMetPlotsWithCut(TCut(("mll>15&&pfJetGoodNumBtag==0&&"+leptCuts[i]).c_str()),"bTagVeto30",metCuts[i], jetCuts[i], nomc,ymax[i]);
501     ProduceMetPlotsWithCut(TCut(("mll>15&&pfJetGoodNumBtag>0&&"+leptCuts[i]).c_str()),"AtLeastOneBJet30",metCuts[i],jetCuts[i], nomc,ymax[i]);
502     ProduceMetPlotsWithCut(TCut(("mll>15&&pfJetGoodNumBtag==0&&"+leptCuts[i]).c_str()),"bTagVeto30",metCuts[i], jetCuts[i], domc,ymax[i]);
503     ProduceMetPlotsWithCut(TCut(("mll>15&&pfJetGoodNumBtag>0&&"+leptCuts[i]).c_str()),"AtLeastOneBJet30",metCuts[i],jetCuts[i], domc,ymax[i]);
504 fronga 1.7 } else {
505 fronga 1.22 ProduceMetPlotsWithCut(TCut(("mll>15&&pfJetGoodNumBtag30==0&&"+leptCuts[i]).c_str()),"bTagVeto30",metCuts[i], jetCuts[i],nomc,ymax[i]);
506     ProduceMetPlotsWithCut(TCut(("mll>15&&pfJetGoodNumBtag30>0&&"+leptCuts[i]).c_str()),"AtLeastOneBJet30",metCuts[i],jetCuts[i],nomc,ymax[i]);
507     ProduceMetPlotsWithCut(TCut(("mll>15&&pfJetGoodNumBtag30==0&&"+leptCuts[i]).c_str()),"bTagVeto30",metCuts[i], jetCuts[i],domc,ymax[i]);
508     ProduceMetPlotsWithCut(TCut(("mll>15&&pfJetGoodNumBtag30>0&&"+leptCuts[i]).c_str()),"AtLeastOneBJet30",metCuts[i], jetCuts[i],domc,ymax[i]);
509 fronga 1.7 }
510     }
511 buchmann 1.1 }
512 buchmann 1.12
513 buchmann 1.13 void compute_r_in_out() {
514     int nbins=10;
515     float xmin=10;
516     float xmax=10;
517     TCut InCut("mll>71&&mll<111&&pfJetGoodNum40>0");
518     TCut OutCut("mll>20&&mll<70&&pfJetGoodNum40>0");
519    
520     TH1F *mllIN =allsamples.Draw("mllIN", "mll",nbins,xmin,xmax,"m_{ee} [GeV]", "events",TCut(cutOSSF&&InCut),mc,PlottingSetup::luminosity,allsamples.FindSample("DYJetsToLLNoTau"));
521     TH1F *mllOUT =allsamples.Draw("mllOUT","mll",nbins,xmin,xmax,"m_{ee} [GeV]", "events",TCut(cutOSSF&&OutCut),mc,PlottingSetup::luminosity,allsamples.FindSample("DYJetsToLLNoTau"));
522    
523     cout << "IN: " << mllIN->Integral() << endl;
524     cout << "OUT: " << mllOUT->Integral() << endl;
525     cout << "r_oi: " <<((float) mllOUT->Integral()) / mllIN->Integral()<< endl;
526    
527     delete mllIN;
528     delete mllOUT;
529    
530     }
531    
532 buchmann 1.17
533    
534     void LabelHisto(TH1 *MET_ratio,string titlex, string titley) {
535     MET_ratio->GetXaxis()->SetTitle(titlex.c_str());
536     MET_ratio->GetXaxis()->CenterTitle();
537     MET_ratio->GetYaxis()->SetTitle(titley.c_str());
538     MET_ratio->GetYaxis()->CenterTitle();
539     }
540    
541 buchmann 1.23 TH1F* GetPredictedAndObservedMetShapes(TCut JetCut, string sPositiveCut,string sNegativeCut,string CorrectedMet,string ObservedMet, string JZBPosvar, string JZBNegvar, float MetCut, int is_data, bool isDYonly, bool isAachen) {
542 buchmann 1.17
543     //Steps:
544     // 1) Prepare samples and histo definition (with "optimal" binning for MET cut)
545     // 2) Fill MET histograms
546     // 3) Fill JZB histograms
547     // 4) Draw them and store them
548     // 5) return predicted MET distribution as is (i.e. not scaled by factor of 2!)
549    
550 buchmann 1.24 cout << "*************************************" << endl;
551 buchmann 1.26 // cout << "** SUMMARY BEFORE STARTING DRAWING **" << endl;
552     // cout << "MET variable: " << ObservedMet << endl;
553     // cout << "Corr. MET var:" << CorrectedMet << endl;
554     // cout << "JZB pos. var: " << JZBPosvar << endl;
555     // cout << "JZB neg. var: " << JZBNegvar << endl;
556     // cout << "JZB pos cut : " << sPositiveCut << endl;
557     // cout << "JZB neg cut : " << sNegativeCut << endl;
558 buchmann 1.17 //Step 1: Prepare samples and histo definition
559     vector<int> SelectedSamples;
560     if(is_data==mc&&isDYonly) {
561     SelectedSamples=allsamples.FindSample("Z_em_DYJetsToL");
562     if(SelectedSamples.size()==0) {
563     write_error(__FUNCTION__,"Cannot continue, there seems to be no DY sample without Taus - goodbye!");
564     assert(SelectedSamples.size()>0);
565     }
566     }
567    
568     float DisplayedBinSize=10.0; // this is the bin size that we use for plotting
569    
570 buchmann 1.21 float BinWidth=1.0;
571 buchmann 1.17 float xmin=0;
572 buchmann 1.23 float xmax=110;
573     if(isAachen) xmax=160;
574 buchmann 1.21 if(MetCut>=xmax) xmax=MetCut+10;
575 buchmann 1.17 int nbins=int((xmax-xmin)/BinWidth);
576 buchmann 1.23
577     float pt2cut=20;
578     if(isAachen)pt2cut=10;
579    
580 buchmann 1.17 stringstream basiccut;
581 buchmann 1.23 basiccut << (const char*) JetCut << "&&" << (const char*) Restrmasscut << "&&" << (const char*) leptoncut << "&&pt1>20&&pt2>" << pt2cut;
582 buchmann 1.17
583    
584     stringstream cMET_observed;
585     cMET_observed << "(" << basiccut.str() << "&&(" << sPositiveCut << ")&&" << (const char*) cutOSSF << ")";
586     stringstream cMET_ttbar_pred;
587     cMET_ttbar_pred << "(" << basiccut.str() << "&&(" << sPositiveCut << ")&&" << (const char*) cutOSOF << ")";
588     stringstream cMET_osof_pred;
589     cMET_osof_pred << "(" << basiccut.str() << "&&(" << sNegativeCut << ")&&" << (const char*) cutOSOF << ")";
590     stringstream cMET_ossf_pred;
591     cMET_ossf_pred << "(" << basiccut.str() << "&&(" << sNegativeCut << ")&&" << (const char*) cutOSSF << ")";
592    
593     write_warning(__FUNCTION__,"Once the rush is over you might want to define the potential sidebands ... ");
594     //Step 2: Fill Met histograms
595     TCanvas *PredictionCanvas = new TCanvas("PredictionCanvas","PredictionCanvas");
596     TH1F *MET_observed = allsamples.Draw("MET_observed",ObservedMet,nbins,xmin,xmax,"MET [GeV]","events",
597 buchmann 1.26 TCut(cMET_observed.str().c_str()),is_data,PlottingSetup::luminosity,SelectedSamples);
598 buchmann 1.17 TH1F *MET_ossf_pred = allsamples.Draw("MET_ossf_pred",CorrectedMet,nbins,xmin,xmax,"MET [GeV]","events",
599 buchmann 1.26 TCut(cMET_ossf_pred.str().c_str()),is_data,PlottingSetup::luminosity,SelectedSamples);
600 buchmann 1.17 TH1F *MET_osof_pred = allsamples.Draw("MET_osof_pred",CorrectedMet,nbins,xmin,xmax,"MET [GeV]","events",
601 buchmann 1.26 TCut(cMET_osof_pred.str().c_str()),is_data,PlottingSetup::luminosity,SelectedSamples);
602 buchmann 1.17 TH1F *MET_ttbar_pred= allsamples.Draw("MET_ttbar_pred",ObservedMet,nbins,xmin,xmax,"MET [GeV]","events",
603 buchmann 1.26 TCut(cMET_ttbar_pred.str().c_str()),is_data,PlottingSetup::luminosity,SelectedSamples);
604 buchmann 1.17
605 buchmann 1.25
606     if(isDYonly && is_data==mc) {
607     TH1F *MET_truth = allsamples.Draw("MET_truth",ObservedMet,1,MetCut,10000,"MET [GeV]","events",TCut(((string)"met[4]>"+any2string(MetCut)).c_str())&&cutOSSF&&TCut(basiccut.str().c_str()),is_data,PlottingSetup::luminosity,SelectedSamples);
608     write_info(__FUNCTION__,"DY Truth is : "+any2string(MET_truth->Integral()));
609     delete MET_truth;
610     }
611    
612    
613 buchmann 1.17 TH1F *MET_predicted=(TH1F*)MET_ossf_pred->Clone("MET_predicted");
614     MET_predicted->Add(MET_osof_pred,-1);
615     MET_predicted->Add(MET_ttbar_pred);
616     MET_predicted->SetLineColor(kRed);
617     MET_observed->SetLineColor(kBlack);
618    
619     TH1F *MET_Z_prediction=(TH1F*)MET_ossf_pred->Clone("MET_Z_prediction");
620     MET_Z_prediction->Add(MET_osof_pred,-1);
621     MET_Z_prediction->SetLineColor(kBlue);
622    
623     LabelHisto(MET_observed,"MET (GeV)","events");
624    
625     //Step 3: Fill JZB histograms
626 buchmann 1.25
627 buchmann 1.17 TH1F *JZB_observed = allsamples.Draw("JZB_observed",JZBPosvar,nbins,xmin,xmax,"JZB [GeV]","events",
628     TCut(cMET_observed.str().c_str()),is_data,PlottingSetup::luminosity,SelectedSamples);
629     TH1F *JZB_ossf_pred = allsamples.Draw("JZB_ossf_pred",JZBNegvar,nbins,xmin,xmax,"JZB [GeV]","events",
630     TCut(cMET_ossf_pred.str().c_str()),is_data,PlottingSetup::luminosity,SelectedSamples);
631     TH1F *JZB_osof_pred = allsamples.Draw("JZB_osof_pred",JZBNegvar,nbins,xmin,xmax,"JZB [GeV]","events",
632     TCut(cMET_osof_pred.str().c_str()),is_data,PlottingSetup::luminosity,SelectedSamples);
633     TH1F *JZB_ttbar_pred= allsamples.Draw("JZB_ttbar_pred",JZBPosvar,nbins,xmin,xmax,"JZB [GeV]","events",
634     TCut(cMET_ttbar_pred.str().c_str()),is_data,PlottingSetup::luminosity,SelectedSamples);
635    
636     TH1F *JZB_predicted=(TH1F*)JZB_ossf_pred->Clone("JZB_predicted");
637     JZB_predicted->Add(JZB_osof_pred,-1);
638     JZB_predicted->Add(JZB_ttbar_pred);
639     JZB_predicted->SetLineColor(kRed);
640     JZB_observed->SetLineColor(kBlack);
641    
642     TH1F *JZB_Z_prediction=(TH1F*)JZB_ossf_pred->Clone("JZB_Z_prediction");
643     JZB_Z_prediction->Add(JZB_osof_pred,-1);
644     MET_Z_prediction->SetLineColor(kBlue);
645    
646     LabelHisto(JZB_observed,"JZB (GeV)","events");
647    
648     // Step 4: Draw them and store them
649     PredictionCanvas->cd();
650     TPad *toppad = new TPad("toppad","toppad",0.0,0.3,1.0,1.0);
651     toppad->Draw();
652     PredictionCanvas->cd();
653     TPad *bottompad = new TPad("bottompad","bottompad",0.0,0.0,1.0,0.3);
654     bottompad->Draw();
655     PredictionCanvas->cd();
656    
657     TLegend *legend = new TLegend(0.6,0.6,0.89,0.89);
658    
659     MET_ttbar_pred->SetLineColor(TColor::GetColor("#005C00"));
660     JZB_ttbar_pred->SetLineColor(TColor::GetColor("#005C00"));
661    
662     legend->SetFillColor(kWhite);
663     legend->SetBorderSize(0);
664     legend->AddEntry(MET_predicted,"prediction","l");
665     legend->AddEntry(MET_observed,"observed","p");
666     legend->AddEntry(MET_Z_prediction,"predicted Z","l");
667     legend->AddEntry(MET_ttbar_pred,"OF-based prediction","l");
668    
669     if(is_data==mc) legend->SetHeader("Simulation:");
670     if(is_data==mc&&isDYonly) legend->SetHeader("DY #rightarrow ee,#mu#mu only:");
671     if(is_data==data) legend->SetHeader("Data:");
672    
673     stringstream SaveJZBname;
674     stringstream SaveMETname;
675     if(is_data==data) {
676     SaveJZBname << "MetPrediction/JZBdistribution_Data_METCutAt" << MetCut;
677     SaveMETname << "MetPrediction/METdistribution_Data_METCutAt" << MetCut;
678     }
679     if(is_data==mc&&!isDYonly) {
680     SaveJZBname << "MetPrediction/JZBdistribution_FullMC_METCutAt" << MetCut;
681     SaveMETname << "MetPrediction/METdistribution_FullMC_METCutAt" << MetCut;
682     }
683     if(is_data==mc&&isDYonly) {
684     SaveJZBname << "MetPrediction/JZBdistribution_DYMC_METCutAt" << MetCut;
685     SaveMETname << "MetPrediction/METdistribution_DYMC_METCutAt" << MetCut;
686     }
687    
688 buchmann 1.26 dout << "Shape summary (MET>50) for ";
689     if(is_data==data) cout << "data";
690     if(is_data==mc&&isDYonly) cout<< "DY ";
691     if(is_data==mc&&!isDYonly) cout << " Full MC";
692     cout << " : " << endl;
693    
694 buchmann 1.24 dout << " observed : " << MET_observed->Integral(MET_observed->FindBin(50),MET_observed->FindBin(xmax)) << endl;
695     dout << " predicted : " << MET_predicted->Integral(MET_predicted->FindBin(50),MET_predicted->FindBin(xmax)) << endl;
696     dout << " Z pred : " << MET_Z_prediction->Integral(MET_Z_prediction->FindBin(50),MET_Z_prediction->FindBin(xmax)) << endl;
697     dout << " ttbar : " << MET_ttbar_pred->Integral(MET_ttbar_pred->FindBin(50),MET_ttbar_pred->FindBin(xmax)) << endl;
698    
699    
700 buchmann 1.17 TH1F *ZpredClone = (TH1F*)MET_Z_prediction->Clone("ZpredClone");
701     ZpredClone->SetLineColor(kBlue);
702     MET_observed->Rebin(int(DisplayedBinSize/BinWidth));
703     ZpredClone->Rebin(int(DisplayedBinSize/BinWidth));
704     MET_predicted->Rebin(int(DisplayedBinSize/BinWidth));
705     MET_ttbar_pred->Rebin(int(DisplayedBinSize/BinWidth));
706    
707     TH1F *JZBZpredClone = (TH1F*)JZB_Z_prediction->Clone("ZpredClone");
708     JZBZpredClone->SetLineColor(kBlue);
709     JZB_observed->Rebin(int(DisplayedBinSize/BinWidth));
710     JZBZpredClone->Rebin(int(DisplayedBinSize/BinWidth));
711     JZB_predicted->Rebin(int(DisplayedBinSize/BinWidth));
712     JZB_ttbar_pred->Rebin(int(DisplayedBinSize/BinWidth));
713    
714     TH1F *JZB_ratio = (TH1F*)JZB_observed->Clone("JZB_ratio");
715     JZB_ratio->Divide(JZB_predicted);
716     LabelHisto(JZB_ratio,"JZB (GeV)","obs/pred");
717     TH1F *MET_ratio = (TH1F*)MET_observed->Clone("MET_ratio");
718     MET_ratio->Divide(MET_predicted);
719 buchmann 1.25 MET_observed->SetMinimum(0.1);
720     MET_observed->SetMaximum(15*MET_observed->GetMaximum());
721     JZB_observed->SetMinimum(0.1);
722     JZB_observed->SetMaximum(15*JZB_observed->GetMaximum());
723 buchmann 1.17 LabelHisto(MET_ratio,"MET (GeV)","obs/pred");
724 buchmann 1.24 TBox *sysenvelope = new TBox(xmin,1.0-MetPlotsSpace::Zprediction_Uncertainty,xmax,1.0+MetPlotsSpace::Zprediction_Uncertainty);
725 buchmann 1.17 sysenvelope->SetFillColor(TColor::GetColor("#82FA58")); // light green
726     sysenvelope->SetLineWidth(0);
727 buchmann 1.24 TBox *dsysenvelope = new TBox(xmin,1.0-2*MetPlotsSpace::Zprediction_Uncertainty,xmax,1.0+2*MetPlotsSpace::Zprediction_Uncertainty);
728     dsysenvelope->SetFillColor(TColor::GetColor("#F3F781")); // light yellow
729     dsysenvelope->SetLineWidth(0);
730 buchmann 1.17
731     toppad->cd();
732 buchmann 1.25 toppad->SetLogy(1);
733 buchmann 1.17 MET_observed->Draw("e1");
734     MET_ttbar_pred->Draw("histo,same");
735     ZpredClone->Draw("histo,same");
736     MET_predicted->Draw("histo,same");
737     MET_observed->Draw("e1,same");
738     legend->Draw();
739     if(is_data==data) DrawPrelim();
740     else DrawMCPrelim();
741    
742     bottompad->cd();
743     MET_ratio->GetYaxis()->SetRangeUser(0,2);
744     MET_ratio->Draw("e1");
745 buchmann 1.24 dsysenvelope->Draw();
746 buchmann 1.17 sysenvelope->Draw();
747     MET_ratio->Draw("AXIS,same");
748     MET_ratio->Draw("e1,same");
749     TLine *metl = new TLine(xmin,1.0,xmax,1.0);
750     metl->SetLineColor(kBlue);
751     metl->Draw();
752     CompleteSave(PredictionCanvas,SaveMETname.str());
753    
754     toppad->cd();
755     JZB_observed->Draw("e1");
756     JZB_ttbar_pred->Draw("histo,same");
757     JZBZpredClone->Draw("histo,same");
758     JZB_predicted->Draw("histo,same");
759     JZB_observed->Draw("e1,same");
760     legend->Draw();
761     if(is_data==data) DrawPrelim();
762     else DrawMCPrelim();
763    
764     bottompad->cd();
765     JZB_ratio->GetYaxis()->SetRangeUser(0,2);
766     JZB_ratio->Draw("e1");
767 buchmann 1.24 dsysenvelope->Draw();
768 buchmann 1.17 sysenvelope->Draw();
769     JZB_ratio->Draw("AXIS,same");
770     JZB_ratio->Draw("e1,same");
771     metl->Draw();
772    
773     CompleteSave(PredictionCanvas,SaveJZBname.str());
774    
775     delete PredictionCanvas;
776     delete MET_observed;
777     delete MET_predicted;
778     //do NOT delete MET_Z_prediction (it's the return value)
779     delete MET_osof_pred;
780     delete MET_ossf_pred;
781     delete MET_ttbar_pred;
782    
783     delete JZB_observed;
784     delete JZB_predicted;
785     delete JZB_osof_pred;
786     delete JZB_ossf_pred;
787     delete JZB_Z_prediction;
788     delete JZB_ttbar_pred;
789    
790     return MET_Z_prediction;
791     }
792    
793 buchmann 1.25 float extract_correction(string jzbvariable) {
794     int position = (int)jzbvariable.find("[1]");
795     if(position==-1) return 0.0;
796     string correction=jzbvariable.substr(position+3,jzbvariable.length()-position-3);
797     position = (int)correction.find("*");
798     if(position==-1) return 0.0;
799     correction=correction.substr(0,position);
800     float correctionvalue=atof(correction.c_str());
801     assert(correctionvalue<1&&correctionvalue>0);
802     return correctionvalue;
803     }
804    
805 buchmann 1.23 float Get_Met_Z_Prediction(TCut JetCut, float MetCut, int isdata, bool isDYonly, bool isAachen=false) {
806 buchmann 1.17 dout << "Going to compute Z region prediction for a MET cut at " << MetCut << " GeV" << endl;
807     // Steps:
808 buchmann 1.25 // 1) Get peak
809     // 2) use the peak and pt correction for sample splitting
810     // and for MET distribution shifting
811 buchmann 1.17 // 3) compute the estimate for MET>MetCut
812    
813     // do this for data if isdata==data, otherwise for MC (full closure if isDYonly==false, otherwise use only DY sample)
814    
815 buchmann 1.25 // Step 1) Get peak
816 buchmann 1.17 float MCPeakNoPtCorr=0,MCPeakErrorNoPtCorr=0,DataPeakNoPtCorr=0,DataPeakErrorNoPtCorr=0,MCSigma=0,DataSigma=0;
817     stringstream resultsNoPtCorr;
818     stringstream NoPtCorrdatajzb;
819     stringstream NoPtCorrmcjzb;
820    
821 buchmann 1.24 if(isAachen) {
822     //need to make sure that none of the typical basic cuts contain problematic selections!
823     string Sleptoncut = (const char*) leptoncut;
824     cout << "Lepton cut is " << Sleptoncut << endl;
825     if((int)Sleptoncut.find("pt2>20")>-1) {
826     write_error(__FUNCTION__,"You're trying to compute the Aachen estimate but are requiring pt2>20 ... please check your config.");
827     assert((int)Sleptoncut.find("pt2>20")==-1);
828     }
829     } else {
830     string Sleptoncut = (const char*) leptoncut;
831     cout << "Lepton cut is " << Sleptoncut << endl;
832     if((int)Sleptoncut.find("pt2>10")>-1) {
833     write_error(__FUNCTION__,"You're trying to compute the ETH estimate but are requiring pt2>10 ... please check your config.");
834     assert((int)Sleptoncut.find("pt2>10")==-1);
835     }
836     }
837    
838    
839 buchmann 1.25 float Ptcorrection=0.0;
840 buchmann 1.24
841 buchmann 1.25 if(isdata==data) Ptcorrection=extract_correction(PlottingSetup::jzbvariabledata);
842     else Ptcorrection=extract_correction(PlottingSetup::jzbvariablemc);
843 buchmann 1.24
844 buchmann 1.25 find_peaks(MCPeakNoPtCorr,MCPeakErrorNoPtCorr, DataPeakNoPtCorr,DataPeakErrorNoPtCorr,resultsNoPtCorr,true,NoPtCorrdatajzb,NoPtCorrmcjzb,(const char*) JetCut, true);
845 buchmann 1.17
846 buchmann 1.25 float PeakPosition=0.0;
847     string jzbvariable;
848 buchmann 1.17 if(isdata==data) {
849     PeakPosition=DataPeakNoPtCorr;
850 buchmann 1.25 jzbvariable=jzbvariabledata;
851 buchmann 1.17 dout << "Found peak in data at " << DataPeakNoPtCorr << " +/- " << DataPeakErrorNoPtCorr << " ; will use this result (" << PeakPosition << ")" << endl;
852     } else {
853     PeakPosition=MCPeakNoPtCorr;
854 buchmann 1.25 jzbvariable=jzbvariablemc;
855 buchmann 1.17 dout << "Found peak in mc at " << MCPeakNoPtCorr << " +/- " << MCPeakErrorNoPtCorr << " ; will use this result (" << PeakPosition << ")" << endl;
856     }
857    
858     // Step 2: Use peak for sample splitting and MET shifting
859 buchmann 1.25 string CorrectedMet="met[4]-"+any2string(Ptcorrection)+"*pt +"+any2string(abs(1.0*(PeakPosition)));
860     if(2*(PeakPosition)<0) CorrectedMet="met[4]-"+any2string(Ptcorrection)+"*pt -"+any2string(abs(1.0*(PeakPosition)));
861 buchmann 1.17
862     stringstream sPositiveCut;
863 buchmann 1.25 if(PeakPosition>0) sPositiveCut << "((" << jzbvariable << "-" << PeakPosition << ")>0)";
864     else sPositiveCut << "( " << jzbvariable << "+" << abs(PeakPosition) << ")>0)";
865 buchmann 1.17
866     stringstream sNegativeCut;
867 buchmann 1.25 if(PeakPosition<0) sNegativeCut << "((" << jzbvariable << "+" << abs(PeakPosition) << ")<0)";
868     else sNegativeCut << "(( " << jzbvariable << "-" << abs(PeakPosition) << ")<0)";
869 buchmann 1.17
870     string ObservedMet="met[4]";
871    
872     stringstream JZBPosvar;
873 buchmann 1.25 JZBPosvar<<jzbvariable;
874     if(PeakPosition>0) JZBPosvar << "-" << PeakPosition;
875     else JZBPosvar << "+" << abs(PeakPosition);
876    
877 buchmann 1.17 stringstream JZBNegvar;
878 buchmann 1.25 JZBNegvar<<"-(" << jzbvariable;
879     if(PeakPosition>0) JZBNegvar << "-" << PeakPosition << ")";
880     else JZBNegvar << "+" << abs(PeakPosition) << ")";
881    
882 buchmann 1.17
883     // Step 3: Compute estimate
884 buchmann 1.23 TH1F *predicted = GetPredictedAndObservedMetShapes(JetCut, sPositiveCut.str(),sNegativeCut.str(),CorrectedMet,ObservedMet,JZBPosvar.str(),JZBNegvar.str(), MetCut, isdata, isDYonly, isAachen);
885 buchmann 1.17 float ZregionZestimate=0;
886     for(int ibin=1;ibin<=(int)predicted->GetNbinsX();ibin++) {
887     if(predicted->GetBinLowEdge(ibin)+predicted->GetBinWidth(ibin)>MetCut) {
888     ZregionZestimate+=2*(predicted->GetBinContent(ibin));
889     }
890     }
891    
892 buchmann 1.25 cout << " Z region estimate in MET>" << MetCut << " for this sample: " << ZregionZestimate << endl;
893    
894    
895 buchmann 1.17 return ZregionZestimate;
896     }
897    
898 buchmann 1.12 void ExperimentalMetPrediction() {
899 buchmann 1.17
900 buchmann 1.23 bool isAachen=false;
901 buchmann 1.24
902     bool HighPurityMode=false; // High Purity = |mll-91.2|<10 GeV , else <20
903    
904     string restrmasscutbkp=(const char*) PlottingSetup::Restrmasscut;
905    
906     if(HighPurityMode) PlottingSetup::Restrmasscut=TCut("abs(mll-91.2)<10");
907     else PlottingSetup::Restrmasscut= TCut("abs(mll-91.2)<20");
908    
909    
910    
911 buchmann 1.23 cout << "Aachen mode (20/10, 2 jets) ? " << isAachen << endl;
912 buchmann 1.24 cout << "High Purity mode? " << HighPurityMode << endl;
913    
914 buchmann 1.23
915     if(isAachen) write_warning(__FUNCTION__,"Please don't forget to adapt the global lepton cut (to 20/10) for Aachen!");
916 buchmann 1.24 stringstream snjets;
917     if(isAachen) snjets << 2;
918     else snjets << 3;
919     float maxMET=100;
920     if(isAachen) maxMET=150;
921 buchmann 1.12
922 buchmann 1.17 TCut nJetsSignal(PlottingSetup::basicqualitycut&&("pfJetGoodNum40>="+snjets.str()).c_str());
923    
924     cout << " ***** TESTING Z PREDICTION ***** " << endl;
925     cout << "Notation (you can copy & paste this to evaluate it further)" << endl;
926     cout << "Cut;Data;MC;DY;" << endl;
927 buchmann 1.26 float DataEstimate = Get_Met_Z_Prediction(Restrmasscut&&nJetsSignal,maxMET, data, false, isAachen);
928 buchmann 1.25 float DYEstimate = Get_Met_Z_Prediction(Restrmasscut&&nJetsSignal,maxMET, mc, true, isAachen);
929 buchmann 1.26 float MCEstimate = Get_Met_Z_Prediction(Restrmasscut&&nJetsSignal,maxMET, mc, false, isAachen);
930 buchmann 1.25
931     cout << maxMET << ";" << DataEstimate << ";" << MCEstimate << ";" << DYEstimate << endl;
932 buchmann 1.24 dout << "Found estimate in data of " << DataEstimate << endl;
933    
934     float Diff=20.0;
935     if(HighPurityMode) Diff=10;
936     TCut cut("mll>20&&pt1>20&&pt2>20");
937     TH1F *zlineshape = allsamples.Draw("zlineshape","mll",int((91.2)*5),20,91.2+20,"m_{ll} (GeV)","events",cutOSSF&&TCut("pfJetGoodNum40==1")&&cut,data,PlottingSetup::luminosity); // bins of 0.2 GeV
938     float a = (zlineshape->Integral(zlineshape->FindBin(20),zlineshape->FindBin(70)));
939     float b = (zlineshape->Integral(zlineshape->FindBin(91.2-Diff),zlineshape->FindBin(91.2+Diff)));
940     float r = a/b;
941     float dr= (a/b)*TMath::Sqrt(1/a+1/b);
942    
943     float SysUncertainty = TMath::Sqrt(DataEstimate*DataEstimate*dr*dr + r*r*(DataEstimate*MetPlotsSpace::Zprediction_Uncertainty*DataEstimate*MetPlotsSpace::Zprediction_Uncertainty));
944     float StatUncertainty = TMath::Sqrt(DataEstimate);
945    
946     cout << "Z estimate in peak : " << DataEstimate << " +/- " << DataEstimate*MetPlotsSpace::Zprediction_Uncertainty << " (sys) +/- " << TMath::Sqrt(2*DataEstimate) << " (stat) " << endl;
947     cout << "Z ESTIMATE IN SR : " << DataEstimate*r << " +/- " << SysUncertainty << " (sys) +/- " << StatUncertainty << " (stat) " << endl;
948     cout << endl;
949     cout << "r = " << r << " +/- " << dr << endl;
950    
951    
952    
953     delete zlineshape;
954    
955     PlottingSetup::Restrmasscut=TCut(restrmasscutbkp.c_str());
956    
957 buchmann 1.12 }
958 buchmann 1.17