ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/MetPlotting.C
Revision: 1.29
Committed: Thu Sep 13 20:39:50 2012 UTC (12 years, 7 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.28: +3 -1 lines
Log Message: 
Only restore previous over/underflow bin status, don't automatically switch it on after peak determination

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