ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/MetPlotting.C
Revision: 1.27
Committed: Thu Sep 13 09:44:34 2012 UTC (12 years, 7 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.26: +2 -0 lines
Log Message: 
Turned off underflow bin globally. Switched on overflow bin for plots requiring it. Removed superfluous pred/obs ratio function which is in the general prediction function anyway; adapted the range of the prediction & observation histos so the overflow bin is shown

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 buchmann 1.27 switch_overunderflow(true);
486 fronga 1.7 float metCuts[] = { 100., 150. };
487 fronga 1.22 float ymax[] = { 90., 80. };
488 fronga 1.7 int jetCuts[] = { 3, 2 };
489 fronga 1.11 //int jetCuts[] = { 3, 3 };
490 fronga 1.22 string leptCuts[] = { "pt1>20&&pt2>20", "pt1>20&&pt2>10" };
491 fronga 1.7 bool nomc(0),domc(1);
492     for ( int i=0; i<2; ++i ) {
493 fronga 1.22 ProduceMetPlotsWithCut(TCut(("mll>15&&"+leptCuts[i]).c_str()),"",metCuts[i],jetCuts[i],nomc,ymax[i]);
494     ProduceMetPlotsWithCut(TCut(("mll>15&&"+leptCuts[i]).c_str()),"",metCuts[i],jetCuts[i],domc,ymax[i]);
495     continue;
496 fronga 1.19 // stringstream jzbcut;
497     // jzbcut << "((is_data&&("<<datajzb<<")>0)||(!is_data&&("<<mcjzb<<")>0))";
498     // ProduceMetPlotsWithCut(TCut((jzbcut.str()+"&&mll>20&&"+leptCuts[i]).c_str()),"JZBpos",metCuts[i], jetCuts[i],domc);
499 fronga 1.11 //continue;
500 fronga 1.7 if (!PlottingSetup::is53reco) { // Old 5_2
501 fronga 1.22 ProduceMetPlotsWithCut(TCut(("mll>15&&pfJetGoodNumBtag==0&&"+leptCuts[i]).c_str()),"bTagVeto30",metCuts[i], jetCuts[i], nomc,ymax[i]);
502     ProduceMetPlotsWithCut(TCut(("mll>15&&pfJetGoodNumBtag>0&&"+leptCuts[i]).c_str()),"AtLeastOneBJet30",metCuts[i],jetCuts[i], nomc,ymax[i]);
503     ProduceMetPlotsWithCut(TCut(("mll>15&&pfJetGoodNumBtag==0&&"+leptCuts[i]).c_str()),"bTagVeto30",metCuts[i], jetCuts[i], domc,ymax[i]);
504     ProduceMetPlotsWithCut(TCut(("mll>15&&pfJetGoodNumBtag>0&&"+leptCuts[i]).c_str()),"AtLeastOneBJet30",metCuts[i],jetCuts[i], domc,ymax[i]);
505 fronga 1.7 } else {
506 fronga 1.22 ProduceMetPlotsWithCut(TCut(("mll>15&&pfJetGoodNumBtag30==0&&"+leptCuts[i]).c_str()),"bTagVeto30",metCuts[i], jetCuts[i],nomc,ymax[i]);
507     ProduceMetPlotsWithCut(TCut(("mll>15&&pfJetGoodNumBtag30>0&&"+leptCuts[i]).c_str()),"AtLeastOneBJet30",metCuts[i],jetCuts[i],nomc,ymax[i]);
508     ProduceMetPlotsWithCut(TCut(("mll>15&&pfJetGoodNumBtag30==0&&"+leptCuts[i]).c_str()),"bTagVeto30",metCuts[i], jetCuts[i],domc,ymax[i]);
509     ProduceMetPlotsWithCut(TCut(("mll>15&&pfJetGoodNumBtag30>0&&"+leptCuts[i]).c_str()),"AtLeastOneBJet30",metCuts[i], jetCuts[i],domc,ymax[i]);
510 fronga 1.7 }
511     }
512 buchmann 1.27 switch_overunderflow(false);
513 buchmann 1.1 }
514 buchmann 1.12
515 buchmann 1.13 void compute_r_in_out() {
516     int nbins=10;
517     float xmin=10;
518     float xmax=10;
519     TCut InCut("mll>71&&mll<111&&pfJetGoodNum40>0");
520     TCut OutCut("mll>20&&mll<70&&pfJetGoodNum40>0");
521    
522     TH1F *mllIN =allsamples.Draw("mllIN", "mll",nbins,xmin,xmax,"m_{ee} [GeV]", "events",TCut(cutOSSF&&InCut),mc,PlottingSetup::luminosity,allsamples.FindSample("DYJetsToLLNoTau"));
523     TH1F *mllOUT =allsamples.Draw("mllOUT","mll",nbins,xmin,xmax,"m_{ee} [GeV]", "events",TCut(cutOSSF&&OutCut),mc,PlottingSetup::luminosity,allsamples.FindSample("DYJetsToLLNoTau"));
524    
525     cout << "IN: " << mllIN->Integral() << endl;
526     cout << "OUT: " << mllOUT->Integral() << endl;
527     cout << "r_oi: " <<((float) mllOUT->Integral()) / mllIN->Integral()<< endl;
528    
529     delete mllIN;
530     delete mllOUT;
531    
532     }
533    
534 buchmann 1.17
535    
536     void LabelHisto(TH1 *MET_ratio,string titlex, string titley) {
537     MET_ratio->GetXaxis()->SetTitle(titlex.c_str());
538     MET_ratio->GetXaxis()->CenterTitle();
539     MET_ratio->GetYaxis()->SetTitle(titley.c_str());
540     MET_ratio->GetYaxis()->CenterTitle();
541     }
542    
543 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) {
544 buchmann 1.17
545     //Steps:
546     // 1) Prepare samples and histo definition (with "optimal" binning for MET cut)
547     // 2) Fill MET histograms
548     // 3) Fill JZB histograms
549     // 4) Draw them and store them
550     // 5) return predicted MET distribution as is (i.e. not scaled by factor of 2!)
551    
552 buchmann 1.24 cout << "*************************************" << endl;
553 buchmann 1.26 // cout << "** SUMMARY BEFORE STARTING DRAWING **" << endl;
554     // cout << "MET variable: " << ObservedMet << endl;
555     // cout << "Corr. MET var:" << CorrectedMet << endl;
556     // cout << "JZB pos. var: " << JZBPosvar << endl;
557     // cout << "JZB neg. var: " << JZBNegvar << endl;
558     // cout << "JZB pos cut : " << sPositiveCut << endl;
559     // cout << "JZB neg cut : " << sNegativeCut << endl;
560 buchmann 1.17 //Step 1: Prepare samples and histo definition
561     vector<int> SelectedSamples;
562     if(is_data==mc&&isDYonly) {
563     SelectedSamples=allsamples.FindSample("Z_em_DYJetsToL");
564     if(SelectedSamples.size()==0) {
565     write_error(__FUNCTION__,"Cannot continue, there seems to be no DY sample without Taus - goodbye!");
566     assert(SelectedSamples.size()>0);
567     }
568     }
569    
570     float DisplayedBinSize=10.0; // this is the bin size that we use for plotting
571    
572 buchmann 1.21 float BinWidth=1.0;
573 buchmann 1.17 float xmin=0;
574 buchmann 1.23 float xmax=110;
575     if(isAachen) xmax=160;
576 buchmann 1.21 if(MetCut>=xmax) xmax=MetCut+10;
577 buchmann 1.17 int nbins=int((xmax-xmin)/BinWidth);
578 buchmann 1.23
579     float pt2cut=20;
580     if(isAachen)pt2cut=10;
581    
582 buchmann 1.17 stringstream basiccut;
583 buchmann 1.23 basiccut << (const char*) JetCut << "&&" << (const char*) Restrmasscut << "&&" << (const char*) leptoncut << "&&pt1>20&&pt2>" << pt2cut;
584 buchmann 1.17
585    
586     stringstream cMET_observed;
587     cMET_observed << "(" << basiccut.str() << "&&(" << sPositiveCut << ")&&" << (const char*) cutOSSF << ")";
588     stringstream cMET_ttbar_pred;
589     cMET_ttbar_pred << "(" << basiccut.str() << "&&(" << sPositiveCut << ")&&" << (const char*) cutOSOF << ")";
590     stringstream cMET_osof_pred;
591     cMET_osof_pred << "(" << basiccut.str() << "&&(" << sNegativeCut << ")&&" << (const char*) cutOSOF << ")";
592     stringstream cMET_ossf_pred;
593     cMET_ossf_pred << "(" << basiccut.str() << "&&(" << sNegativeCut << ")&&" << (const char*) cutOSSF << ")";
594    
595     write_warning(__FUNCTION__,"Once the rush is over you might want to define the potential sidebands ... ");
596     //Step 2: Fill Met histograms
597     TCanvas *PredictionCanvas = new TCanvas("PredictionCanvas","PredictionCanvas");
598     TH1F *MET_observed = allsamples.Draw("MET_observed",ObservedMet,nbins,xmin,xmax,"MET [GeV]","events",
599 buchmann 1.26 TCut(cMET_observed.str().c_str()),is_data,PlottingSetup::luminosity,SelectedSamples);
600 buchmann 1.17 TH1F *MET_ossf_pred = allsamples.Draw("MET_ossf_pred",CorrectedMet,nbins,xmin,xmax,"MET [GeV]","events",
601 buchmann 1.26 TCut(cMET_ossf_pred.str().c_str()),is_data,PlottingSetup::luminosity,SelectedSamples);
602 buchmann 1.17 TH1F *MET_osof_pred = allsamples.Draw("MET_osof_pred",CorrectedMet,nbins,xmin,xmax,"MET [GeV]","events",
603 buchmann 1.26 TCut(cMET_osof_pred.str().c_str()),is_data,PlottingSetup::luminosity,SelectedSamples);
604 buchmann 1.17 TH1F *MET_ttbar_pred= allsamples.Draw("MET_ttbar_pred",ObservedMet,nbins,xmin,xmax,"MET [GeV]","events",
605 buchmann 1.26 TCut(cMET_ttbar_pred.str().c_str()),is_data,PlottingSetup::luminosity,SelectedSamples);
606 buchmann 1.17
607 buchmann 1.25
608     if(isDYonly && is_data==mc) {
609     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);
610     write_info(__FUNCTION__,"DY Truth is : "+any2string(MET_truth->Integral()));
611     delete MET_truth;
612     }
613    
614    
615 buchmann 1.17 TH1F *MET_predicted=(TH1F*)MET_ossf_pred->Clone("MET_predicted");
616     MET_predicted->Add(MET_osof_pred,-1);
617     MET_predicted->Add(MET_ttbar_pred);
618     MET_predicted->SetLineColor(kRed);
619     MET_observed->SetLineColor(kBlack);
620    
621     TH1F *MET_Z_prediction=(TH1F*)MET_ossf_pred->Clone("MET_Z_prediction");
622     MET_Z_prediction->Add(MET_osof_pred,-1);
623     MET_Z_prediction->SetLineColor(kBlue);
624    
625     LabelHisto(MET_observed,"MET (GeV)","events");
626    
627     //Step 3: Fill JZB histograms
628 buchmann 1.25
629 buchmann 1.17 TH1F *JZB_observed = allsamples.Draw("JZB_observed",JZBPosvar,nbins,xmin,xmax,"JZB [GeV]","events",
630     TCut(cMET_observed.str().c_str()),is_data,PlottingSetup::luminosity,SelectedSamples);
631     TH1F *JZB_ossf_pred = allsamples.Draw("JZB_ossf_pred",JZBNegvar,nbins,xmin,xmax,"JZB [GeV]","events",
632     TCut(cMET_ossf_pred.str().c_str()),is_data,PlottingSetup::luminosity,SelectedSamples);
633     TH1F *JZB_osof_pred = allsamples.Draw("JZB_osof_pred",JZBNegvar,nbins,xmin,xmax,"JZB [GeV]","events",
634     TCut(cMET_osof_pred.str().c_str()),is_data,PlottingSetup::luminosity,SelectedSamples);
635     TH1F *JZB_ttbar_pred= allsamples.Draw("JZB_ttbar_pred",JZBPosvar,nbins,xmin,xmax,"JZB [GeV]","events",
636     TCut(cMET_ttbar_pred.str().c_str()),is_data,PlottingSetup::luminosity,SelectedSamples);
637    
638     TH1F *JZB_predicted=(TH1F*)JZB_ossf_pred->Clone("JZB_predicted");
639     JZB_predicted->Add(JZB_osof_pred,-1);
640     JZB_predicted->Add(JZB_ttbar_pred);
641     JZB_predicted->SetLineColor(kRed);
642     JZB_observed->SetLineColor(kBlack);
643    
644     TH1F *JZB_Z_prediction=(TH1F*)JZB_ossf_pred->Clone("JZB_Z_prediction");
645     JZB_Z_prediction->Add(JZB_osof_pred,-1);
646     MET_Z_prediction->SetLineColor(kBlue);
647    
648     LabelHisto(JZB_observed,"JZB (GeV)","events");
649    
650     // Step 4: Draw them and store them
651     PredictionCanvas->cd();
652     TPad *toppad = new TPad("toppad","toppad",0.0,0.3,1.0,1.0);
653     toppad->Draw();
654     PredictionCanvas->cd();
655     TPad *bottompad = new TPad("bottompad","bottompad",0.0,0.0,1.0,0.3);
656     bottompad->Draw();
657     PredictionCanvas->cd();
658    
659     TLegend *legend = new TLegend(0.6,0.6,0.89,0.89);
660    
661     MET_ttbar_pred->SetLineColor(TColor::GetColor("#005C00"));
662     JZB_ttbar_pred->SetLineColor(TColor::GetColor("#005C00"));
663    
664     legend->SetFillColor(kWhite);
665     legend->SetBorderSize(0);
666     legend->AddEntry(MET_predicted,"prediction","l");
667     legend->AddEntry(MET_observed,"observed","p");
668     legend->AddEntry(MET_Z_prediction,"predicted Z","l");
669     legend->AddEntry(MET_ttbar_pred,"OF-based prediction","l");
670    
671     if(is_data==mc) legend->SetHeader("Simulation:");
672     if(is_data==mc&&isDYonly) legend->SetHeader("DY #rightarrow ee,#mu#mu only:");
673     if(is_data==data) legend->SetHeader("Data:");
674    
675     stringstream SaveJZBname;
676     stringstream SaveMETname;
677     if(is_data==data) {
678     SaveJZBname << "MetPrediction/JZBdistribution_Data_METCutAt" << MetCut;
679     SaveMETname << "MetPrediction/METdistribution_Data_METCutAt" << MetCut;
680     }
681     if(is_data==mc&&!isDYonly) {
682     SaveJZBname << "MetPrediction/JZBdistribution_FullMC_METCutAt" << MetCut;
683     SaveMETname << "MetPrediction/METdistribution_FullMC_METCutAt" << MetCut;
684     }
685     if(is_data==mc&&isDYonly) {
686     SaveJZBname << "MetPrediction/JZBdistribution_DYMC_METCutAt" << MetCut;
687     SaveMETname << "MetPrediction/METdistribution_DYMC_METCutAt" << MetCut;
688     }
689    
690 buchmann 1.26 dout << "Shape summary (MET>50) for ";
691     if(is_data==data) cout << "data";
692     if(is_data==mc&&isDYonly) cout<< "DY ";
693     if(is_data==mc&&!isDYonly) cout << " Full MC";
694     cout << " : " << endl;
695    
696 buchmann 1.24 dout << " observed : " << MET_observed->Integral(MET_observed->FindBin(50),MET_observed->FindBin(xmax)) << endl;
697     dout << " predicted : " << MET_predicted->Integral(MET_predicted->FindBin(50),MET_predicted->FindBin(xmax)) << endl;
698     dout << " Z pred : " << MET_Z_prediction->Integral(MET_Z_prediction->FindBin(50),MET_Z_prediction->FindBin(xmax)) << endl;
699     dout << " ttbar : " << MET_ttbar_pred->Integral(MET_ttbar_pred->FindBin(50),MET_ttbar_pred->FindBin(xmax)) << endl;
700    
701    
702 buchmann 1.17 TH1F *ZpredClone = (TH1F*)MET_Z_prediction->Clone("ZpredClone");
703     ZpredClone->SetLineColor(kBlue);
704     MET_observed->Rebin(int(DisplayedBinSize/BinWidth));
705     ZpredClone->Rebin(int(DisplayedBinSize/BinWidth));
706     MET_predicted->Rebin(int(DisplayedBinSize/BinWidth));
707     MET_ttbar_pred->Rebin(int(DisplayedBinSize/BinWidth));
708    
709     TH1F *JZBZpredClone = (TH1F*)JZB_Z_prediction->Clone("ZpredClone");
710     JZBZpredClone->SetLineColor(kBlue);
711     JZB_observed->Rebin(int(DisplayedBinSize/BinWidth));
712     JZBZpredClone->Rebin(int(DisplayedBinSize/BinWidth));
713     JZB_predicted->Rebin(int(DisplayedBinSize/BinWidth));
714     JZB_ttbar_pred->Rebin(int(DisplayedBinSize/BinWidth));
715    
716     TH1F *JZB_ratio = (TH1F*)JZB_observed->Clone("JZB_ratio");
717     JZB_ratio->Divide(JZB_predicted);
718     LabelHisto(JZB_ratio,"JZB (GeV)","obs/pred");
719     TH1F *MET_ratio = (TH1F*)MET_observed->Clone("MET_ratio");
720     MET_ratio->Divide(MET_predicted);
721 buchmann 1.25 MET_observed->SetMinimum(0.1);
722     MET_observed->SetMaximum(15*MET_observed->GetMaximum());
723     JZB_observed->SetMinimum(0.1);
724     JZB_observed->SetMaximum(15*JZB_observed->GetMaximum());
725 buchmann 1.17 LabelHisto(MET_ratio,"MET (GeV)","obs/pred");
726 buchmann 1.24 TBox *sysenvelope = new TBox(xmin,1.0-MetPlotsSpace::Zprediction_Uncertainty,xmax,1.0+MetPlotsSpace::Zprediction_Uncertainty);
727 buchmann 1.17 sysenvelope->SetFillColor(TColor::GetColor("#82FA58")); // light green
728     sysenvelope->SetLineWidth(0);
729 buchmann 1.24 TBox *dsysenvelope = new TBox(xmin,1.0-2*MetPlotsSpace::Zprediction_Uncertainty,xmax,1.0+2*MetPlotsSpace::Zprediction_Uncertainty);
730     dsysenvelope->SetFillColor(TColor::GetColor("#F3F781")); // light yellow
731     dsysenvelope->SetLineWidth(0);
732 buchmann 1.17
733     toppad->cd();
734 buchmann 1.25 toppad->SetLogy(1);
735 buchmann 1.17 MET_observed->Draw("e1");
736     MET_ttbar_pred->Draw("histo,same");
737     ZpredClone->Draw("histo,same");
738     MET_predicted->Draw("histo,same");
739     MET_observed->Draw("e1,same");
740     legend->Draw();
741     if(is_data==data) DrawPrelim();
742     else DrawMCPrelim();
743    
744     bottompad->cd();
745     MET_ratio->GetYaxis()->SetRangeUser(0,2);
746     MET_ratio->Draw("e1");
747 buchmann 1.24 dsysenvelope->Draw();
748 buchmann 1.17 sysenvelope->Draw();
749     MET_ratio->Draw("AXIS,same");
750     MET_ratio->Draw("e1,same");
751     TLine *metl = new TLine(xmin,1.0,xmax,1.0);
752     metl->SetLineColor(kBlue);
753     metl->Draw();
754     CompleteSave(PredictionCanvas,SaveMETname.str());
755    
756     toppad->cd();
757     JZB_observed->Draw("e1");
758     JZB_ttbar_pred->Draw("histo,same");
759     JZBZpredClone->Draw("histo,same");
760     JZB_predicted->Draw("histo,same");
761     JZB_observed->Draw("e1,same");
762     legend->Draw();
763     if(is_data==data) DrawPrelim();
764     else DrawMCPrelim();
765    
766     bottompad->cd();
767     JZB_ratio->GetYaxis()->SetRangeUser(0,2);
768     JZB_ratio->Draw("e1");
769 buchmann 1.24 dsysenvelope->Draw();
770 buchmann 1.17 sysenvelope->Draw();
771     JZB_ratio->Draw("AXIS,same");
772     JZB_ratio->Draw("e1,same");
773     metl->Draw();
774    
775     CompleteSave(PredictionCanvas,SaveJZBname.str());
776    
777     delete PredictionCanvas;
778     delete MET_observed;
779     delete MET_predicted;
780     //do NOT delete MET_Z_prediction (it's the return value)
781     delete MET_osof_pred;
782     delete MET_ossf_pred;
783     delete MET_ttbar_pred;
784    
785     delete JZB_observed;
786     delete JZB_predicted;
787     delete JZB_osof_pred;
788     delete JZB_ossf_pred;
789     delete JZB_Z_prediction;
790     delete JZB_ttbar_pred;
791    
792     return MET_Z_prediction;
793     }
794    
795 buchmann 1.25 float extract_correction(string jzbvariable) {
796     int position = (int)jzbvariable.find("[1]");
797     if(position==-1) return 0.0;
798     string correction=jzbvariable.substr(position+3,jzbvariable.length()-position-3);
799     position = (int)correction.find("*");
800     if(position==-1) return 0.0;
801     correction=correction.substr(0,position);
802     float correctionvalue=atof(correction.c_str());
803     assert(correctionvalue<1&&correctionvalue>0);
804     return correctionvalue;
805     }
806    
807 buchmann 1.23 float Get_Met_Z_Prediction(TCut JetCut, float MetCut, int isdata, bool isDYonly, bool isAachen=false) {
808 buchmann 1.17 dout << "Going to compute Z region prediction for a MET cut at " << MetCut << " GeV" << endl;
809     // Steps:
810 buchmann 1.25 // 1) Get peak
811     // 2) use the peak and pt correction for sample splitting
812     // and for MET distribution shifting
813 buchmann 1.17 // 3) compute the estimate for MET>MetCut
814    
815     // do this for data if isdata==data, otherwise for MC (full closure if isDYonly==false, otherwise use only DY sample)
816    
817 buchmann 1.25 // Step 1) Get peak
818 buchmann 1.17 float MCPeakNoPtCorr=0,MCPeakErrorNoPtCorr=0,DataPeakNoPtCorr=0,DataPeakErrorNoPtCorr=0,MCSigma=0,DataSigma=0;
819     stringstream resultsNoPtCorr;
820     stringstream NoPtCorrdatajzb;
821     stringstream NoPtCorrmcjzb;
822    
823 buchmann 1.24 if(isAachen) {
824     //need to make sure that none of the typical basic cuts contain problematic selections!
825     string Sleptoncut = (const char*) leptoncut;
826     cout << "Lepton cut is " << Sleptoncut << endl;
827     if((int)Sleptoncut.find("pt2>20")>-1) {
828     write_error(__FUNCTION__,"You're trying to compute the Aachen estimate but are requiring pt2>20 ... please check your config.");
829     assert((int)Sleptoncut.find("pt2>20")==-1);
830     }
831     } else {
832     string Sleptoncut = (const char*) leptoncut;
833     cout << "Lepton cut is " << Sleptoncut << endl;
834     if((int)Sleptoncut.find("pt2>10")>-1) {
835     write_error(__FUNCTION__,"You're trying to compute the ETH estimate but are requiring pt2>10 ... please check your config.");
836     assert((int)Sleptoncut.find("pt2>10")==-1);
837     }
838     }
839    
840    
841 buchmann 1.25 float Ptcorrection=0.0;
842 buchmann 1.24
843 buchmann 1.25 if(isdata==data) Ptcorrection=extract_correction(PlottingSetup::jzbvariabledata);
844     else Ptcorrection=extract_correction(PlottingSetup::jzbvariablemc);
845 buchmann 1.24
846 buchmann 1.25 find_peaks(MCPeakNoPtCorr,MCPeakErrorNoPtCorr, DataPeakNoPtCorr,DataPeakErrorNoPtCorr,resultsNoPtCorr,true,NoPtCorrdatajzb,NoPtCorrmcjzb,(const char*) JetCut, true);
847 buchmann 1.17
848 buchmann 1.25 float PeakPosition=0.0;
849     string jzbvariable;
850 buchmann 1.17 if(isdata==data) {
851     PeakPosition=DataPeakNoPtCorr;
852 buchmann 1.25 jzbvariable=jzbvariabledata;
853 buchmann 1.17 dout << "Found peak in data at " << DataPeakNoPtCorr << " +/- " << DataPeakErrorNoPtCorr << " ; will use this result (" << PeakPosition << ")" << endl;
854     } else {
855     PeakPosition=MCPeakNoPtCorr;
856 buchmann 1.25 jzbvariable=jzbvariablemc;
857 buchmann 1.17 dout << "Found peak in mc at " << MCPeakNoPtCorr << " +/- " << MCPeakErrorNoPtCorr << " ; will use this result (" << PeakPosition << ")" << endl;
858     }
859    
860     // Step 2: Use peak for sample splitting and MET shifting
861 buchmann 1.25 string CorrectedMet="met[4]-"+any2string(Ptcorrection)+"*pt +"+any2string(abs(1.0*(PeakPosition)));
862     if(2*(PeakPosition)<0) CorrectedMet="met[4]-"+any2string(Ptcorrection)+"*pt -"+any2string(abs(1.0*(PeakPosition)));
863 buchmann 1.17
864     stringstream sPositiveCut;
865 buchmann 1.25 if(PeakPosition>0) sPositiveCut << "((" << jzbvariable << "-" << PeakPosition << ")>0)";
866     else sPositiveCut << "( " << jzbvariable << "+" << abs(PeakPosition) << ")>0)";
867 buchmann 1.17
868     stringstream sNegativeCut;
869 buchmann 1.25 if(PeakPosition<0) sNegativeCut << "((" << jzbvariable << "+" << abs(PeakPosition) << ")<0)";
870     else sNegativeCut << "(( " << jzbvariable << "-" << abs(PeakPosition) << ")<0)";
871 buchmann 1.17
872     string ObservedMet="met[4]";
873    
874     stringstream JZBPosvar;
875 buchmann 1.25 JZBPosvar<<jzbvariable;
876     if(PeakPosition>0) JZBPosvar << "-" << PeakPosition;
877     else JZBPosvar << "+" << abs(PeakPosition);
878    
879 buchmann 1.17 stringstream JZBNegvar;
880 buchmann 1.25 JZBNegvar<<"-(" << jzbvariable;
881     if(PeakPosition>0) JZBNegvar << "-" << PeakPosition << ")";
882     else JZBNegvar << "+" << abs(PeakPosition) << ")";
883    
884 buchmann 1.17
885     // Step 3: Compute estimate
886 buchmann 1.23 TH1F *predicted = GetPredictedAndObservedMetShapes(JetCut, sPositiveCut.str(),sNegativeCut.str(),CorrectedMet,ObservedMet,JZBPosvar.str(),JZBNegvar.str(), MetCut, isdata, isDYonly, isAachen);
887 buchmann 1.17 float ZregionZestimate=0;
888     for(int ibin=1;ibin<=(int)predicted->GetNbinsX();ibin++) {
889     if(predicted->GetBinLowEdge(ibin)+predicted->GetBinWidth(ibin)>MetCut) {
890     ZregionZestimate+=2*(predicted->GetBinContent(ibin));
891     }
892     }
893    
894 buchmann 1.25 cout << " Z region estimate in MET>" << MetCut << " for this sample: " << ZregionZestimate << endl;
895    
896    
897 buchmann 1.17 return ZregionZestimate;
898     }
899    
900 buchmann 1.12 void ExperimentalMetPrediction() {
901 buchmann 1.17
902 buchmann 1.23 bool isAachen=false;
903 buchmann 1.24
904     bool HighPurityMode=false; // High Purity = |mll-91.2|<10 GeV , else <20
905    
906     string restrmasscutbkp=(const char*) PlottingSetup::Restrmasscut;
907    
908     if(HighPurityMode) PlottingSetup::Restrmasscut=TCut("abs(mll-91.2)<10");
909     else PlottingSetup::Restrmasscut= TCut("abs(mll-91.2)<20");
910    
911    
912    
913 buchmann 1.23 cout << "Aachen mode (20/10, 2 jets) ? " << isAachen << endl;
914 buchmann 1.24 cout << "High Purity mode? " << HighPurityMode << endl;
915    
916 buchmann 1.23
917     if(isAachen) write_warning(__FUNCTION__,"Please don't forget to adapt the global lepton cut (to 20/10) for Aachen!");
918 buchmann 1.24 stringstream snjets;
919     if(isAachen) snjets << 2;
920     else snjets << 3;
921     float maxMET=100;
922     if(isAachen) maxMET=150;
923 buchmann 1.12
924 buchmann 1.17 TCut nJetsSignal(PlottingSetup::basicqualitycut&&("pfJetGoodNum40>="+snjets.str()).c_str());
925    
926     cout << " ***** TESTING Z PREDICTION ***** " << endl;
927     cout << "Notation (you can copy & paste this to evaluate it further)" << endl;
928     cout << "Cut;Data;MC;DY;" << endl;
929 buchmann 1.26 float DataEstimate = Get_Met_Z_Prediction(Restrmasscut&&nJetsSignal,maxMET, data, false, isAachen);
930 buchmann 1.25 float DYEstimate = Get_Met_Z_Prediction(Restrmasscut&&nJetsSignal,maxMET, mc, true, isAachen);
931 buchmann 1.26 float MCEstimate = Get_Met_Z_Prediction(Restrmasscut&&nJetsSignal,maxMET, mc, false, isAachen);
932 buchmann 1.25
933     cout << maxMET << ";" << DataEstimate << ";" << MCEstimate << ";" << DYEstimate << endl;
934 buchmann 1.24 dout << "Found estimate in data of " << DataEstimate << endl;
935    
936     float Diff=20.0;
937     if(HighPurityMode) Diff=10;
938     TCut cut("mll>20&&pt1>20&&pt2>20");
939     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
940     float a = (zlineshape->Integral(zlineshape->FindBin(20),zlineshape->FindBin(70)));
941     float b = (zlineshape->Integral(zlineshape->FindBin(91.2-Diff),zlineshape->FindBin(91.2+Diff)));
942     float r = a/b;
943     float dr= (a/b)*TMath::Sqrt(1/a+1/b);
944    
945     float SysUncertainty = TMath::Sqrt(DataEstimate*DataEstimate*dr*dr + r*r*(DataEstimate*MetPlotsSpace::Zprediction_Uncertainty*DataEstimate*MetPlotsSpace::Zprediction_Uncertainty));
946     float StatUncertainty = TMath::Sqrt(DataEstimate);
947    
948     cout << "Z estimate in peak : " << DataEstimate << " +/- " << DataEstimate*MetPlotsSpace::Zprediction_Uncertainty << " (sys) +/- " << TMath::Sqrt(2*DataEstimate) << " (stat) " << endl;
949     cout << "Z ESTIMATE IN SR : " << DataEstimate*r << " +/- " << SysUncertainty << " (sys) +/- " << StatUncertainty << " (stat) " << endl;
950     cout << endl;
951     cout << "r = " << r << " +/- " << dr << endl;
952    
953    
954    
955     delete zlineshape;
956    
957     PlottingSetup::Restrmasscut=TCut(restrmasscutbkp.c_str());
958    
959 buchmann 1.12 }
960 buchmann 1.17