ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/MetPlotting.C
Revision: 1.32
Committed: Wed Sep 19 09:16:24 2012 UTC (12 years, 7 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.31: +32 -7 lines
Log Message:
Updated window to 10 GeV around 91 GeV

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