ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/MetPlotting.C
Revision: 1.39
Committed: Tue Nov 13 13:44:03 2012 UTC (12 years, 5 months ago) by fronga
Content type: text/plain
Branch: MAIN
Changes since 1.38: +131 -58 lines
Log Message: 
Added MC expectation for CR. Some cleanup.

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