ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/MetPlotting.C
Revision: 1.38
Committed: Thu Oct 11 09:05:24 2012 UTC (12 years, 6 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.37: +0 -2 lines
Log Message:
Added y range for plots; added frederic's Z+Jets r_{in/out} ratio; added check for Aachen selection; made it more verbose; added pt spectrum for comparison

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