ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/ZbTools.C
Revision: 1.12
Committed: Mon Nov 19 16:21:52 2012 UTC (12 years, 5 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.11: +53 -13 lines
Log Message: 
Added inclusive selection so we can compute results also for C_{abs}^{inc}; note that mistag and efficiency corrections are switched off in this case

File Contents

# User Rev Content
1 buchmann 1.1 #include <iostream>
2     #include <vector>
3     #include <sys/stat.h>
4    
5     #include <TCut.h>
6     #include <TROOT.h>
7     #include <TCanvas.h>
8     #include <TMath.h>
9     #include <TColor.h>
10     #include <TPaveText.h>
11     #include <TRandom.h>
12     #include <TH1.h>
13     #include <TH2.h>
14     #include <TF1.h>
15     #include <TSQLResult.h>
16     #include <TProfile.h>
17     #include <TLegendEntry.h>
18    
19     using namespace std;
20    
21     using namespace PlottingSetup;
22    
23     namespace ZbData {
24     vector<float> data_over_mc;
25    
26 buchmann 1.10 TCut ZplusBsel("pt1>20&&pt2>20&&mll>60&&mll<120&&id1==id2");
27 buchmann 1.9 TCut LeadingB("Zb3010_bTagProbCSVBP[0]>0.898");
28     TCut EtaB("abs(Zb3010_pfJetGoodEta[0])<1.3");
29 buchmann 1.10 TCut PhiZcut("abs(Zb3010_pfJetDphiZ[0])>2.7");
30 buchmann 1.1
31     }
32    
33     using namespace ZbData;
34    
35 buchmann 1.7
36     //*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/* DELETE EVERYTHING BETWEEN THIS LINE /*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*//
37     void SpecialCutFlow(string identifier);
38    
39     //*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/* AND THIS LINE /*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*//
40    
41    
42 buchmann 1.1 void print_yield(TCut cut) {
43     TH1F *data = allsamples.Draw("data", "mll",1,0,10000000, "m_{ll} [GeV]", "events", cut,data,luminosity);
44 buchmann 1.7 dout << " Yields for " << (const char*) cut << endl;
45     dout << " data: " << data->Integral() << endl;
46 buchmann 1.9 THStack mcm = allsamples.DrawStack("mc", "mll",1,0,10000000, "m_{ll} [GeV]", "events", cut,mc,luminosity);
47 buchmann 1.1 int nhists = mcm.GetHists()->GetEntries();
48 buchmann 1.7 dout << "Going to loop over " << nhists << " histograms " << endl;
49 buchmann 1.1 TFile *test = new TFile("test.root","RECREATE");
50     mcm.Write();
51     test->Close();
52     for (int i = 0; i < nhists; i++) {
53     TH1* hist = static_cast<TH1*>(mcm.GetHists()->At(i));
54     cout << " " << hist->GetName() << ": " << hist->Integral() << endl;
55     }
56 buchmann 1.9
57 buchmann 1.1 cout << endl << endl;
58 buchmann 1.9
59 buchmann 1.1 delete data;
60     }
61    
62 buchmann 1.2 void draw_kin_variable(string variable, TCut cut, int nbins, float min, float max, string xlabel, string saveas, bool dologscale, string speciallabel="") {
63 buchmann 1.1
64     TCanvas *ckin = new TCanvas("ckin","kin variable canvas");
65     ckin->SetLogy(dologscale);
66     TH1F *datahisto = allsamples.Draw("datahisto",variable,nbins,min,max, xlabel, "events",cut,data,luminosity);
67     datahisto->SetMarkerSize(DataMarkerSize);
68     THStack mcstack = allsamples.DrawStack("mcstack",variable,nbins,min,max, xlabel, "events",cut,mc,luminosity);
69 buchmann 1.2 if(dologscale) datahisto->SetMaximum(3*datahisto->GetMaximum());
70     else datahisto->SetMaximum(1.5*datahisto->GetMaximum());
71     TText *speciall;
72     if(speciallabel!="") {
73     speciall = new TText(0.2,0.8,speciallabel.c_str());
74     speciall->Draw();
75     }
76 buchmann 1.1 datahisto->Draw("e1");
77     ckin->Update();
78     mcstack.Draw("same");
79     datahisto->Draw("same,e1");
80     TLegend *kinleg = allsamples.allbglegend();
81     kinleg->Draw();
82    
83     TPad *kinpad = new TPad("kinpad","kinpad",0,0,1,1);
84     kinpad->cd();
85     kinpad->SetLogy(dologscale);
86     datahisto->Draw("e1");
87     mcstack.Draw("same");
88     datahisto->Draw("same,e1");
89     datahisto->Draw("same,axis");
90     kinleg->Draw();
91     DrawPrelim();
92     saveas="kin/"+saveas;
93     save_with_ratio(datahisto,mcstack,kinpad->cd(),saveas);
94    
95    
96    
97    
98     ZbData::data_over_mc.push_back(datahisto->Integral()/CollapseStack(mcstack)->Integral());
99    
100     datahisto->Delete();
101     delete ckin;
102     }
103    
104     Value getfrom2Dmap(TH2F *map, int ixbin, int iybin) {
105     Value sum(0,0);
106     for(int iy=1;iy<=iybin;iy++) {
107     sum=sum+Value(map->GetBinContent(ixbin,iy),map->GetBinError(ixbin,iy));
108     }
109     return sum;
110     }
111    
112     void print_all_b_yields() {
113     cout << "Basic selection with a b jet" << endl;
114 buchmann 1.9 print_yield(ZplusBsel&&"Zb3010_pfJetGoodNumBtag>0");
115 buchmann 1.2 cout << "Leading jet is a b-jet " << endl;
116 buchmann 1.1 print_yield(ZplusBsel&&LeadingB);
117     cout << "|#eta_{b}|<1.3 " << endl;
118     print_yield(ZplusBsel&&LeadingB&&EtaB);
119     cout << "|#delta#phi(b<Z)|>2.7" << endl;
120     print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut);
121 buchmann 1.2 cout << "10<ptZ<1000 GeV" << endl;
122 buchmann 1.1 print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000");
123     cout << "#alpha < 0.3" << endl;
124 buchmann 1.9 print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.3");
125 buchmann 1.1 cout << "#alpha < 0.2" << endl;
126 buchmann 1.9 print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.2");
127 buchmann 1.1 cout << "#alpha < 0.15" << endl;
128 buchmann 1.9 print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.15");
129 buchmann 1.1 cout << "#alpha < 0.1" << endl;
130 buchmann 1.9 print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.1");
131     // cout << "#alpha < 0.05" << endl;
132     // print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.05");
133 buchmann 1.1 }
134    
135 buchmann 1.7 void DrawEvilCutFlow() {
136     stringstream MegaCut;
137    
138 buchmann 1.9 MegaCut << " 1*(" << (const char*) (ZplusBsel&&TCut("Zb3010_pfJetGoodNumBtag>0")) << ")";
139 buchmann 1.7 MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB) << ")";
140     MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB) << ")";
141     MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut) << ")";
142     MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>10&&pt<1000")) << ")";
143 buchmann 1.9 MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>10&&pt<1000&&Zb3010_alpha<0.3")) << ")";
144     MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>10&&pt<1000&&Zb3010_alpha<0.2")) << ")";
145     MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>10&&pt<1000&&Zb3010_alpha<0.15")) << ")";
146     MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>10&&pt<1000&&Zb3010_alpha<0.1")) << ")";
147     MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>10&&pt<1000&&Zb3010_alpha<0.05")) << ")";
148 buchmann 1.7 MegaCut << " ";
149    
150     TCanvas *can = new TCanvas("can","can");
151     can->SetLogy(1);
152     TCut basecut("mll>6||mll<10");
153    
154    
155    
156     TLegend *leg = allsamples.allbglegend();
157    
158    
159     TH1F *data = allsamples.Draw("data", MegaCut.str(),11,-0.5,10.5, "", "events", basecut,data,luminosity);
160     THStack mcm = allsamples.DrawStack("mc", MegaCut.str(),11,-0.5,10.5, "", "events", basecut,mc,luminosity);
161    
162     float runningsum=0;
163     for(int i=data->GetNbinsX();i>0;i--) {
164     runningsum+=data->GetBinContent(i);
165     data->SetBinContent(i,runningsum);
166     }
167    
168     float minimum=data->GetMaximum();
169    
170     TH1F* h;
171     TIter nextOF(mcm.GetHists());
172     while ( h = (TH1F*)nextOF() ) {
173     float runningsum=0;
174     minimum=data->GetMaximum();//done deliberately at each step so get the minimum of the last (!) contribution
175     for(int i=h->GetNbinsX();i>0;i--) {
176     runningsum+=h->GetBinContent(i);
177     h->SetBinContent(i,runningsum);
178     if(runningsum<minimum&&runningsum>0) minimum=runningsum;
179     }
180     }
181    
182    
183    
184     data->SetMinimum(0.05*minimum);
185     can->cd(1)->SetBottomMargin(0.25);
186     data->GetXaxis()->SetBinLabel(1,"Pre-selection");
187     data->GetXaxis()->SetBinLabel(2,"Contains b jet");
188     data->GetXaxis()->SetBinLabel(3,"Leading jet is b-jet");
189     data->GetXaxis()->SetBinLabel(4,"|#eta_{b}|<1.3 ");
190     data->GetXaxis()->SetBinLabel(5,"|#delta#phi(b,Z)|>2.7");
191     data->GetXaxis()->SetBinLabel(6,"10<p_{T}^{Z}<1000 GeV");
192     data->GetXaxis()->SetBinLabel(7,"#alpha < 0.3");
193     data->GetXaxis()->SetBinLabel(8,"#alpha < 0.2");
194     data->GetXaxis()->SetBinLabel(9,"#alpha < 0.15");
195     data->GetXaxis()->SetBinLabel(10,"#alpha < 0.1");
196     data->GetXaxis()->SetBinLabel(11,"#alpha < 0.05");
197     data->GetXaxis()->LabelsOption("v");
198    
199     data->Draw("e1");
200     mcm.Draw("same");
201     data->Draw("same,e1,axis");
202     data->Draw("same,e1");
203     // data->Draw("same,TEXT");
204    
205     leg->Draw();
206    
207     DrawPrelim();
208    
209     CompleteSave(can,"CutFlow");
210     }
211    
212 buchmann 1.1 void draw_Zb_kin_vars() {
213 buchmann 1.6 draw_kin_variable("pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",25,0,200,"Z p_{T}","Official/Zpt",1);
214 buchmann 1.9 draw_kin_variable("Zb3010_pfJetGoodPt[1]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",25,0,200,"Sub-Leading Jet Pt","Official/Jet2Pt",1);
215     draw_kin_variable("Zb3010_pfJetGoodPt[0]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",25,0,200,"Leading Jet Pt (B)","Official/LeadingJetPt",1);
216     draw_kin_variable("Zb3010_pfJetGoodPt[1]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,0,2,"p_{T}^{J2} / p_{T}^{Z}","Official/SubLeadingJetPt_Over_ZPt",1);
217     draw_kin_variable("Zb3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<50", 40,0,2,"#alpha","Official/alpha_pt_30_to_50",1);
218     draw_kin_variable("Zb3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>50&&pt<75", 40,0,2,"#alpha","Official/alpha_50_to_75",1);
219     draw_kin_variable("Zb3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>75&&pt<125", 40,0,2,"#alpha","Official/alpha_pt_75_to_125",1);
220     draw_kin_variable("Zb3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>125&&pt<1000",40,0,2,"#alpha","Official/alpha_pt_125_to_1000",1);
221     draw_kin_variable("Zb3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000",40,0,2,"#alpha","Official/alpha_pt_30_to_1000",1);
222 buchmann 1.6
223 buchmann 1.10 draw_kin_variable("Zb3010_pfBJetDphiZ[0]",ZplusBsel&&LeadingB&&"pt>10&&pt<1000&&pfBJetDphiZ[0]>0",30,0,3.2,"#delta#phi (Z,b lead)","DeltaPhiZBlead",0);
224 buchmann 1.6
225     /*
226     draw_kin_variable("pt1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,0,200,"p_{T}^{l1}","Lep/pt1",1);
227     draw_kin_variable("pt1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id1==0",20,0,200,"p_{T}^{e1}","Lep/pt1_e",1);
228     draw_kin_variable("pt1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id1==1",20,0,200,"p_{T}^{#mu1}","Lep/pt1_m",1);
229     draw_kin_variable("pt2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,0,200,"p_{T}^{l2}","Lep/pt2",1);
230     draw_kin_variable("pt2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id2==0",20,0,200,"p_{T}^{e2}","Lep/pt2_e",1);
231     draw_kin_variable("pt2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id2==1",20,0,200,"p_{T}^{#mu2}","Lep/pt2_m",1);
232     draw_kin_variable("eta1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,-3.1,3.1,"#eta^{l1}","Lep/eta1",1);
233     draw_kin_variable("eta1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id1==0",20,-3.1,3.1,"#eta^{e1}","Lep/eta1_e",1);
234     draw_kin_variable("eta1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id1==1",20,-3.1,3.1,"#eta^{m1}","Lep/eta1_m",1);
235     draw_kin_variable("eta2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,-3.1,3.1,"#eta^{l1}","Lep/eta2",1);
236     draw_kin_variable("eta2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id2==0",20,-3.1,3.1,"#eta^{e2}","Lep/eta2_e",1);
237     draw_kin_variable("eta2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id2==1",20,-3.1,3.1,"#eta^{#mu2}","Lep/eta2_m",1);
238     draw_kin_variable("numVtx",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",30,0,30,"N(Vtx)","Lep/nVtx",1);
239     */
240 buchmann 1.1 }
241    
242     void draw_mpf_vars() {
243     ZbData::data_over_mc.clear();
244 buchmann 1.9 draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.3",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p3",0,"#alpha<0.3, 10 GeV < p_{T}^{Z} < 1000 GeV");
245     draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.2",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p2",0,"#alpha<0.2, 10 GeV < p_{T}^{Z} < 1000 GeV");
246     draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.15",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p15",0,"#alpha<0.15, 10 GeV < p_{T}^{Z} < 1000 GeV");
247     draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.1",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p1",0,"#alpha<0.1, 10 GeV < p_{T}^{Z} < 1000 GeV");
248     draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.05",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p05",0,"#alpha<0.05, 10 GeV < p_{T}^{Z} < 1000 GeV");
249    
250     draw_kin_variable("Zb3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.3",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p3",0,"#alpha<0.3, 10 GeV < p_{T}^{Z} < 1000 GeV");
251     draw_kin_variable("Zb3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.2",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p2",0,"#alpha<0.2, 10 GeV < p_{T}^{Z} < 1000 GeV");
252     draw_kin_variable("Zb3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.15",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p15",0,"#alpha<0.15, 10 GeV < p_{T}^{Z} < 1000 GeV");
253     draw_kin_variable("Zb3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.1",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p1",0,"#alpha<0.1, 10 GeV < p_{T}^{Z} < 1000 GeV");
254     draw_kin_variable("Zb3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"Zb3010_alpha<0.05",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz__mpf_alpha_smaller_0p05",0,"#alpha<0.05, 10 GeV < p_{T}^{Z} < 1000 GeV");
255 buchmann 1.1 }
256    
257 buchmann 1.12 /*double GetMedian(TH1F *histo) {
258 buchmann 1.11 int numBins = histo->GetXaxis()->GetNbins();
259 buchmann 1.12 Double_t x[numBins];
260     Double_t y[numBins];
261 buchmann 1.11 for (int i = 1; i <= numBins; i++) {
262     x[i] = histo->GetBinCenter(i);
263     y[i] = histo->GetBinContent(i);
264     }
265 buchmann 1.12 return TMath::Median(numBins, &x, &y);
266     }*/
267 buchmann 1.11
268 buchmann 1.12 Value get_Zb_data_over_mc(string variable, string variable2, TCut cut, string saveas) {
269 buchmann 1.2 // write_warning(__FUNCTION__,"Debugging this function, therefore always returning 3.1415+/-1.010101"); return Value(3.1415,1.010101);
270 buchmann 1.12 TCanvas *cn = new TCanvas("cn","cn");
271     string varname="MPF";
272     if(Contains(variable,"pfJetGood")) varname="R_{abs}";
273     TH1F *hdata = allsamples.Draw("hdata",variable,1200,0,30, varname, "events", cut,data,luminosity);
274     TH1F *hmc = allsamples.Draw("hmc" , variable,1200,0,30, varname, "events", cut, mc, luminosity);
275 buchmann 1.2
276 buchmann 1.10 float a=hdata->GetMean();
277 buchmann 1.2 float b=hmc->GetMean();
278     float da=hdata->GetMeanError();
279 buchmann 1.10 float db=hmc->GetMeanError();
280 buchmann 1.2 float factor = a / b;
281     float error = TMath::Sqrt( (1/(b*b)) * (da*da) + ((a*a)/(b*b))*db*db);
282 buchmann 1.12
283     cn->cd();
284     hdata->GetYaxis()->SetTitle("events / 0.1");
285     hdata->Rebin(4);
286     hmc->Rebin(4);
287     hdata->GetXaxis()->SetRangeUser(0,2);
288     hdata->Draw("e1");
289     hmc->Draw("histo,same");
290     hdata->Draw("e1,same");
291     stringstream summary;
292     summary << "#splitline{" << varname << "}{#splitline{data: " << std::setprecision(4) << a << " #pm " << std::setprecision(4) << da << "}{#splitline{mc: " << std::setprecision(4) << b << " #pm " << std::setprecision(4) << db << "}{ratio: " << factor;
293     summary << " #pm " << error << "}}}";
294     TText *infobox = write_title(summary.str());
295     infobox->SetX(0.75);
296     infobox->SetTextSize(0.03);
297     infobox->SetY(0.75);
298     infobox->Draw();
299    
300    
301     DrawPrelim();
302     CompleteSave(cn,"ResponsePlots/"+saveas);
303    
304     cout << "Data;" << a << ";"<<da<<";MC"<<b<<";"<<db<<endl;
305     delete cn;
306 buchmann 1.2 delete hdata;
307     delete hmc;
308 buchmann 1.9
309     // cout << a << "+/-" << da << " " << b << "=/-" << db << endl;
310    
311 buchmann 1.2 return Value(factor,error);
312     }
313    
314    
315    
316     void new_data_mc_agreement_2d() {
317 buchmann 1.3 gStyle->SetOptFit(0);
318 buchmann 1.2 TCut basecut(ZplusBsel&&LeadingB&&EtaB&&PhiZcut);
319 buchmann 1.9 const int nalphacuts=5;
320     float alphacuts[nalphacuts] = {0.1,0.15,0.2,0.25,0.3};
321 buchmann 1.5 int PointThree=-1;
322     for(int i=0;i<nalphacuts;i++) {
323 buchmann 1.9 if((alphacuts[i]-0.3)<0.01) PointThree=i-1;
324 buchmann 1.5 }
325     if(PointThree<0) {
326     write_warning(__FUNCTION__,"alpha=0.3 MUST BE ONE OF THE ALPHA VALUES!");
327     assert(PointThree<0);
328     }
329    
330 buchmann 1.2 const int nptcuts=5;
331 buchmann 1.9 float ptcuts[nptcuts]={30,50,75,125,1000};
332 buchmann 1.2
333     float MPF_Results[nalphacuts][nptcuts];
334     float MPF_Errors[nalphacuts][nptcuts];
335     float RABS_Results[nalphacuts][nptcuts];
336     float RABS_Errors[nalphacuts][nptcuts];
337    
338    
339 buchmann 1.9 for(int ia=0;ia<nalphacuts;ia++) {
340 buchmann 1.2 for(int ipt=0;ipt<nptcuts-1;ipt++) {
341     stringstream specialcut;
342 buchmann 1.9 specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (Zb3010_alpha<" << alphacuts[ia] << "))";
343     cout << specialcut.str() << endl;
344 buchmann 1.12 Value MPF_data_over_mc = get_Zb_data_over_mc("mpf","",TCut(basecut && specialcut.str().c_str()),"MPF___pt_"+any2string(ptcuts[ipt])+"_"+any2string(ptcuts[ipt+1])+"__alpha_"+any2string(alphacuts[ia]));
345     Value RABS_data_over_mc = get_Zb_data_over_mc("Zb3010_pfJetGoodPt[0]/pt","",TCut(basecut && specialcut.str().c_str()),"Rabs___pt_"+any2string(ptcuts[ipt])+"_"+any2string(ptcuts[ipt+1])+"__alpha_"+any2string(alphacuts[ia]));
346 buchmann 1.2
347     MPF_Results[ia][ipt]=MPF_data_over_mc.getValue();
348     MPF_Errors[ia][ipt]=MPF_data_over_mc.getError();
349    
350     RABS_Results[ia][ipt]=RABS_data_over_mc.getValue();
351     RABS_Errors[ia][ipt]=RABS_data_over_mc.getError();
352    
353 buchmann 1.9 cout << "alpha cut at " << alphacuts[ia+1] << " and pt range " << ptcuts[ipt] << " , " << ptcuts[ipt+1] << " \t : \t" << specialcut.str() << " \t \t --> " << MPF_data_over_mc << " (MPF) " << RABS_data_over_mc << " (RABS)" << endl;
354 buchmann 1.2 }
355     }
356    
357     float MPF_ExtraPolatedResults[nptcuts];
358     float RABS_ExtraPolatedResults[nptcuts];
359    
360     TCanvas *can = new TCanvas("can");
361    
362     TGraphErrors *MPF_FinalGraph = new TGraphErrors();
363     TGraphErrors *RABS_FinalGraph = new TGraphErrors();
364 buchmann 1.5 TGraphErrors *FaceValueAtPoint3 = new TGraphErrors();
365 buchmann 1.2
366     for(int ipt=0;ipt<nptcuts-1;ipt++) {
367    
368     stringstream filename;
369     filename << "Extrapolation/Zplusb_data_over_mc___" << ptcuts[ipt] << "_to_" << ptcuts[ipt+1];
370     cout << "Going to create histo for " << filename.str() << endl;
371     TGraphErrors *mpf_gr =new TGraphErrors();
372     TGraphErrors *rabs_gr =new TGraphErrors();
373     for(int ia=0;ia<nalphacuts;ia++) {
374 buchmann 1.9 cout << "Setting a point for an alpha cut at " << alphacuts[ia] << " with value " << MPF_Results[ia][ipt] << " (MPF) and " << RABS_Results[ia][ipt] << " (RABS)" << endl;
375     mpf_gr->SetPoint(ia,alphacuts[ia],MPF_Results[ia][ipt]); // avoiding empty bins (nan's)
376     mpf_gr->SetPointError(ia,0,MPF_Errors[ia][ipt]); // avoiding empty bins (nan's)
377     rabs_gr->SetPoint(ia,alphacuts[ia],RABS_Results[ia][ipt]); // avoiding empty bins (nan's)
378     rabs_gr->SetPointError(ia,0,RABS_Errors[ia][ipt]); // avoiding empty bins (nan's)
379 buchmann 1.2 }
380 buchmann 1.9 cout << "About to add a point in pt (" << ptcuts[ipt]+ptcuts[ipt+1] << " ) to face value : " << MPF_Results[PointThree][ipt] << endl;
381 buchmann 1.5 FaceValueAtPoint3->SetPoint(ipt,(ptcuts[ipt]+ptcuts[ipt+1]),MPF_Results[PointThree][ipt]);
382     FaceValueAtPoint3->SetPointError(ipt,0,MPF_Errors[PointThree][ipt]);
383 buchmann 1.2
384     mpf_gr->SetMarkerStyle(21);
385     mpf_gr->SetMarkerColor(kBlue);
386     mpf_gr->Draw("AP");
387 buchmann 1.3 mpf_gr->Fit("pol1","QE"); // quiet, use minos
388     mpf_gr->GetXaxis()->SetTitle("#alpha");
389     mpf_gr->GetXaxis()->CenterTitle();
390 buchmann 1.9 mpf_gr->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
391 buchmann 1.3 mpf_gr->GetYaxis()->CenterTitle();
392 buchmann 1.2 TF1 *mpf_pol=(TF1*)mpf_gr->GetFunction("pol1");
393     float mpf_a=mpf_pol->GetParameter(0);
394     float mpf_b=mpf_pol->GetParameter(1);
395     MPF_FinalGraph->SetPoint(ipt,0.5*(ptcuts[ipt]+ptcuts[ipt+1]),mpf_a);
396     MPF_FinalGraph->SetPointError(ipt,0,mpf_pol->GetParError(0));
397    
398     MPF_ExtraPolatedResults[ipt]=mpf_a;
399    
400     stringstream mpf_info;
401 buchmann 1.3 mpf_info << "#splitline{#splitline{#splitline{" << ptcuts[ipt] << " GeV < p_{T}^{Z} < " << ptcuts[ipt+1] << " GeV }{ (MPF) }}{Extrapolated value: " << std::setprecision(3) << mpf_a << "}}{#chi^{2}/NDF : " << mpf_pol->GetChisquare() << " / " << mpf_pol->GetNDF() << "}";
402    
403 buchmann 1.2 TText *mpf_mark = write_title(mpf_info.str());
404     mpf_mark->SetX(0.75);
405     mpf_mark->SetTextSize(0.03);
406     mpf_mark->SetY(0.75);
407     mpf_mark->Draw();
408    
409     string filenamebkp=filename.str();
410     filename << "__MPF";
411 buchmann 1.3 DrawPrelim();
412 buchmann 1.2 CompleteSave(can,filename.str());
413     cout << "MPF : " << mpf_a << " + " << mpf_b << " * alpha " << endl;
414    
415     filename.str("");
416    
417     rabs_gr->SetMarkerStyle(21);
418     rabs_gr->SetMarkerColor(kBlue);
419 buchmann 1.3 rabs_gr->Fit("pol1","QE"); // quiet, use minos
420 buchmann 1.2 rabs_gr->Draw("AP");
421 buchmann 1.3 rabs_gr->GetXaxis()->SetTitle("#alpha");
422     rabs_gr->GetXaxis()->CenterTitle();
423 buchmann 1.9 rabs_gr->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
424 buchmann 1.3 rabs_gr->GetYaxis()->CenterTitle();
425    
426 buchmann 1.2 TF1 *rabs_pol=(TF1*)rabs_gr->GetFunction("pol1");
427     float rabs_a=rabs_pol->GetParameter(0);
428     float rabs_b=rabs_pol->GetParameter(1);
429     RABS_FinalGraph->SetPoint(ipt,0.5*(ptcuts[ipt]+ptcuts[ipt+1]),rabs_a);
430     RABS_FinalGraph->SetPointError(ipt,0,rabs_pol->GetParError(0));
431 buchmann 1.3
432 buchmann 1.2 cout << "!+!+!+!+!+!+!+!+!+ Just added a point to the final plot : " << rabs_a << " +/- " << rabs_pol->GetParError(0) << endl;
433    
434    
435     stringstream rabs_info;
436 buchmann 1.3 rabs_info << "#splitline{#splitline{#splitline{" << ptcuts[ipt] << " GeV < p_{T}^{Z} < " << ptcuts[ipt+1] << " GeV }{ (R_{abs}) }}{Extrapolated value: " << std::setprecision(3) << rabs_a << "}}{#chi^{2}/NDF : " << rabs_pol->GetChisquare() << " / " << rabs_pol->GetNDF() << "}";
437 buchmann 1.2 TText *rabs_mark = write_title(rabs_info.str());
438     rabs_mark->SetX(0.75);
439     rabs_mark->SetTextSize(0.03);
440     rabs_mark->SetY(0.75);
441     rabs_mark->Draw();
442 buchmann 1.3
443    
444 buchmann 1.2
445     RABS_ExtraPolatedResults[ipt]=rabs_a;
446     filename << filenamebkp << "__RABS";
447 buchmann 1.3 DrawPrelim();
448 buchmann 1.2 CompleteSave(can,filename.str());
449     cout << "RABS : " << rabs_a << " + " << rabs_b << " * alpha " << endl;
450     }
451    
452     MPF_FinalGraph->SetMarkerStyle(21);
453     MPF_FinalGraph->SetMarkerColor(kRed);
454 buchmann 1.3 MPF_FinalGraph->Fit("pol0","QE"); // quiet, use minos
455    
456 buchmann 1.2 TF1 *mpf_pol0=(TF1*)MPF_FinalGraph->GetFunction("pol0");
457     MPF_FinalGraph->Draw("AP");
458 buchmann 1.3 MPF_FinalGraph->GetXaxis()->SetTitle("p_{T}^{Z}");
459     MPF_FinalGraph->GetXaxis()->CenterTitle();
460 buchmann 1.9 MPF_FinalGraph->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
461 buchmann 1.3 MPF_FinalGraph->GetYaxis()->CenterTitle();
462     MPF_FinalGraph->Draw("AP");
463    
464 buchmann 1.2 stringstream mpf_result;
465 buchmann 1.3 mpf_result << "#splitline{C_{abs}= " << std::setprecision(5) << 1/mpf_pol0->GetParameter(0) << " #pm " << std::setprecision(2) << mpf_pol0->GetParError(0)/(mpf_pol0->GetParameter(0)*mpf_pol0->GetParameter(0)) << "}{ (MPF)}";
466 buchmann 1.2 TText *rmpf_final_mark = write_title(mpf_result.str());
467     rmpf_final_mark->SetX(0.75);
468     rmpf_final_mark->SetY(0.75);
469     rmpf_final_mark->SetTextSize(0.03);
470     rmpf_final_mark->Draw();
471 buchmann 1.3 DrawPrelim();
472 buchmann 1.5 float MPFResult=1/mpf_pol0->GetParameter(0);
473     float MPFResultError=mpf_pol0->GetParError(0)/(mpf_pol0->GetParameter(0)*mpf_pol0->GetParameter(0));
474 buchmann 1.2
475     stringstream filename2;
476     filename2 << "Extrapolation/FINAL_RESULT_MPF";
477     CompleteSave(can,filename2.str());
478    
479    
480     RABS_FinalGraph->SetMarkerStyle(21);
481     RABS_FinalGraph->SetMarkerStyle(21);
482     RABS_FinalGraph->SetMarkerColor(kRed);
483 buchmann 1.3 RABS_FinalGraph->Fit("pol0","QE"); // quiet, use minos
484     RABS_FinalGraph->Draw("AP");
485     RABS_FinalGraph->GetXaxis()->SetTitle("p_{T}^{Z}");
486     RABS_FinalGraph->GetXaxis()->CenterTitle();
487 buchmann 1.9 RABS_FinalGraph->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
488 buchmann 1.3 RABS_FinalGraph->GetYaxis()->CenterTitle();
489 buchmann 1.2 RABS_FinalGraph->Draw("AP");
490     TF1 *rabs_pol0=(TF1*)RABS_FinalGraph->GetFunction("pol0");
491     stringstream rabs_result;
492 buchmann 1.3 rabs_result << "#splitline{C_{abs}= " << std::setprecision(5) << 1/rabs_pol0->GetParameter(0) << " #pm " << std::setprecision(2) << rabs_pol0->GetParError(0)/(rabs_pol0->GetParameter(0)*rabs_pol0->GetParameter(0)) << "}{ (R_{abs})}";;
493 buchmann 1.2 TText *rabs_final_mark = write_title(rabs_result.str());
494     rabs_final_mark->SetX(0.75);
495     rabs_final_mark->SetY(0.75);
496     rabs_final_mark->SetTextSize(0.03);
497     rabs_final_mark->Draw();
498 buchmann 1.3 DrawPrelim();
499 buchmann 1.2
500 buchmann 1.5 float RabsResult=1/rabs_pol0->GetParameter(0);
501     float RabsResultError=rabs_pol0->GetParError(0)/(rabs_pol0->GetParameter(0)*rabs_pol0->GetParameter(0));
502    
503     filename2.str("");
504     filename2 << "Extrapolation/FINAL_RESULT_RABS";
505     CompleteSave(can,filename2.str());
506    
507 buchmann 1.9 can->cd();
508 buchmann 1.5 TCanvas *FaceValueCanvas = new TCanvas("FaceValueCanvas","FaceValueCanvas");
509     FaceValueAtPoint3->SetMarkerStyle(21);
510     FaceValueAtPoint3->SetMarkerStyle(21);
511     FaceValueAtPoint3->Fit("pol0","QE"); // quiet, use minos
512     FaceValueAtPoint3->Draw("AP");
513     FaceValueAtPoint3->GetXaxis()->SetTitle("p_{T}^{Z}");
514     FaceValueAtPoint3->GetXaxis()->CenterTitle();
515     FaceValueAtPoint3->GetYaxis()->SetTitle("data/mc for #alpha<0.3");
516     FaceValueAtPoint3->GetYaxis()->CenterTitle();
517     FaceValueAtPoint3->Draw("AP");
518     TF1 *faceval_pol0=(TF1*)FaceValueAtPoint3->GetFunction("pol0");
519     stringstream faceval_result;
520     faceval_result << "#splitline{C_{abs}= " << std::setprecision(5) << 1/faceval_pol0->GetParameter(0) << " #pm " << std::setprecision(2) << faceval_pol0->GetParError(0)/(faceval_pol0->GetParameter(0)*faceval_pol0->GetParameter(0)) << "}{ (MPF at #alpha<0.3)}";;
521     TText *faceval_final_mark = write_title(faceval_result.str());
522     faceval_final_mark->SetX(0.75);
523     faceval_final_mark->SetY(0.75);
524     faceval_final_mark->SetTextSize(0.03);
525     faceval_final_mark->Draw();
526     DrawPrelim();
527    
528     float FaceValResult=1/faceval_pol0->GetParameter(0);
529     float FaceValResultError=faceval_pol0->GetParError(0)/(faceval_pol0->GetParameter(0)*faceval_pol0->GetParameter(0));
530    
531 buchmann 1.2 filename2.str("");
532 buchmann 1.9 filename2 << "Extrapolation/FINAL_RESULT_MPF_FaceValp3";
533 buchmann 1.2 CompleteSave(can,filename2.str());
534    
535 buchmann 1.5 cout << "FINAL RESULTS: " << endl;
536     cout << "MPF : " << MPFResult << " +/- " << MPFResultError << endl;
537     cout << "Rabs: " << RabsResult << " +/- " << RabsResultError << endl;
538     cout << "Face value at 0.3 :" << FaceValResult << " +/- " << FaceValResultError << endl;
539 buchmann 1.2
540     delete can;
541    
542    
543     }
544    
545    
546 buchmann 1.9 void DoPUStudy(string identifier) {
547    
548     float numVtxcuts[7]={0,10,15,20};
549    
550     TCanvas *can = new TCanvas("can","can");
551     TH1F *malpha[8];
552     TH1F *dalpha[8];
553     cout << identifier << ";";
554     cout << endl;
555     for(int i=1;i<5;i++) {
556     stringstream sSpecialCut;
557     if(i<4) {
558     sSpecialCut << "numVtx>" << numVtxcuts[i-1] << " && numVtx<" << numVtxcuts[i];
559     cout << numVtxcuts[i-1] << " < nVtx < " << numVtxcuts[i] << ";";
560     } else {
561     sSpecialCut << "numVtx>" << numVtxcuts[i-1];
562     cout << numVtxcuts[i-1] << ";";
563     }
564    
565     TCut SpecialCut(sSpecialCut.str().c_str());
566    
567    
568     malpha[i] = allsamples.Draw("malpha"+any2string(i),"mpf",1000,0,10,"MPF","events",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&SpecialCut,mc,luminosity);
569     dalpha[i] = allsamples.Draw("dalpha"+any2string(i),"mpf",1000,0,10,"MPF","events",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&SpecialCut,data,luminosity);
570    
571    
572     float a=dalpha[i]->GetMean();
573     float b=malpha[i]->GetMean();
574     float da=dalpha[i]->GetMeanError();
575     float db=malpha[i]->GetMeanError();
576     float factor = a / b;
577     float error = TMath::Sqrt( (1/(b*b)) * (da*da) + ((a*a)/(b*b))*db*db);
578    
579     cout << malpha[i]->GetMean() << ";" << malpha[i]->GetMeanError() << ";" << dalpha[i]->GetMean() << ";" << dalpha[i]->GetMeanError() << ";" << malpha[i]->Integral()<<";"<<dalpha[i]->Integral() << ";"<<factor << ";" << error << endl;
580     malpha[i]->SetLineColor(i+1);
581     if(i==0) malpha[i]->Draw("histo");
582     else malpha[i]->Draw("histo,same");
583     }
584    
585     CompleteSave(can,"AlphaOverview_"+identifier);
586     delete can;
587    
588     /* for(int i=1;i<8;i++) {
589     delete malpha[i];
590     delete dalpha[i];
591     }
592     */
593    
594     }
595    
596     void ScenarioComparison() {
597    
598     TGraphAsymmErrors *gr[5];
599     TF1 *fit[5];
600     bool dofit=false;
601    
602     TCanvas *can = new TCanvas("can","can",500,500);
603    
604     TLegend *leg = new TLegend(0.2,0.2,0.5,0.45);
605     leg->SetBorderSize(0);
606     leg->SetFillColor(kWhite);
607    
608     for(int i=0;i<5;i++) {
609     gr[i] = new TGraphAsymmErrors();
610     double x[4] = {5,12.5,17.5,20};
611     double dxl[4] = {5,2.5,2.5,2.5};
612     double dxh[4] = {5,2.5,2.5,10};
613    
614     string name="";
615     string color="";
616    
617     if(i==0) {
618     double y[4] = {1.008,1.013,1.003,1.007};
619     double dyl[4] = {0.011,0.009,0.013,0.020};
620     double dyh[4] = {0.011,0.009,0.013,0.020};
621     name="30/30";
622     gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
623     color="#a01d99";
624     }
625     if(i==1) {
626     double y[4] = {1.008+0.001,1.013+0.001,1.003+0.001,1.007+0.001};
627     double dyl[4] = {0.011,0.009,0.013,0.020};
628     double dyh[4] = {0.011,0.009,0.013,0.020};
629     gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
630     name="30/10";
631     color="#2464bc";
632     }
633     if(i==2) {
634     double y[4] = {1.014,1.026,1.014,1.059};
635     double dyl[4] = {0.011,0.010,0.014,0.023};
636     double dyh[4] = {0.011,0.010,0.014,0.023};
637     gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
638     name="15/15";
639     color="#f49230";
640     }
641     if(i==3) {
642     double y[4] = {1.015,1.020,1.008,1.050};
643     double dyl[4] = {0.011,0.010,0.014,0.022};
644     double dyh[4] = {0.011,0.010,0.014,0.022};
645     gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
646     name="20/15";
647     color="#c72f3c";
648     }
649     if(i==4) {
650     double y[4] = {1.015+0.001,1.020+0.001,1.008+0.001,1.050+0.001};
651     double dyl[4] = {0.011,0.010,0.014,0.022};
652     double dyh[4] = {0.011,0.010,0.014,0.022};
653     gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
654     name="20/20";
655     color="#7cb433";
656     }
657    
658    
659     gr[i]->SetTitle();
660     gr[i]->GetXaxis()->SetTitle("N(Vtx)");
661     gr[i]->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
662     gr[i]->GetXaxis()->CenterTitle();
663     gr[i]->GetYaxis()->CenterTitle();
664     gr[i]->SetName(name.c_str());
665     gr[i]->SetMarkerSize(0.00001);
666     gr[i]->GetYaxis()->SetRangeUser(0.90,1.10);
667     gr[i]->SetLineColor(TColor::GetColor(color.c_str()));
668     gr[i]->SetMarkerColor(TColor::GetColor(color.c_str()));
669     if(dofit) {
670     gr[i]->Fit("pol1");
671     fit[i] = (TF1*)gr[i]->GetFunction("pol1");
672     fit[i]->SetLineColor(TColor::GetColor(color.c_str()));
673     }
674     leg->AddEntry(gr[i],name.c_str(),"l");
675     if(i==0) gr[i]->Draw("AP*");
676     else gr[i]->Draw("P");
677     }
678     leg->Draw();
679     DrawPrelim();
680    
681     CompleteSave(can,"ScenarioComparison");
682     delete can;
683    
684     }
685    
686     void ScenarioComparisonInclusive() {
687    
688     TGraphAsymmErrors *gr[5];
689     TF1 *fit[5];
690     bool dofit=true;
691    
692     TCanvas *can = new TCanvas("can","can",500,500);
693    
694     TLegend *leg = new TLegend(0.2,0.2,0.5,0.45);
695     leg->SetBorderSize(0);
696     leg->SetFillColor(kWhite);
697    
698     for(int i=0;i<5;i++) {
699     gr[i] = new TGraphAsymmErrors();
700     double x[4] = {5,12.5,17.5,20};
701     double dxl[4] = {5,2.5,2.5,2.5};
702     double dxh[4] = {5,2.5,2.5,10};
703    
704     string name="";
705     string color="";
706    
707     if(i==0) {
708     double y[4] = {1.01103,1.01496,1.01486,1.02785};
709     double dyl[4] = {0.00235044,0.00225264,0.00313574,0.00541047};
710     double dyh[4] = {0.00235044,0.00225264,0.00313574,0.00541047};
711     name="30/30";
712     gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
713     color="#a01d99";
714     }
715     if(i==1) {
716     double y[4] = {1.01103+0.001,1.01496+0.001,1.01486+0.001,1.02785+0.001};
717     double dyl[4] = {0.00235044,0.00225264,0.00313574,0.00541047};
718     double dyh[4] = {0.00235044,0.00225264,0.00313574,0.00541047};
719     gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
720     name="30/10";
721     color="#2464bc";
722     }
723     if(i==2) {
724     double y[4] = {1.01575,1.02879,1.0325,1.05456};
725     double dyl[4] = {0.00288281,0.00278536,0.00367451,0.00590421};
726     double dyh[4] = {0.00288281,0.00278536,0.00367451,0.00590421};
727     gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
728     name="15/15";
729     color="#f49230";
730     }
731     if(i==3) {
732     double y[4] = {1.01221,1.01922,1.02396,1.04677};
733     double dyl[4] = {0.00263821,0.00260853,0.00353673,0.00584728};
734     double dyh[4] = {0.00263821,0.00260853,0.00353673,0.00584728};
735     gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
736     name="20/15";
737     color="#c72f3c";
738     }
739     if(i==4) {
740     double y[4] = {1.01221+0.001,1.01922+0.001,1.02396+0.001,1.04677+0.001};
741     double dyl[4] = {0.00263821,0.00260853,0.00353673,0.00584728};
742     double dyh[4] = {0.00263821,0.00260853,0.00353673,0.00584728};
743     gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
744     name="20/20";
745     color="#7cb433";
746     }
747    
748    
749     gr[i]->SetTitle();
750     gr[i]->GetXaxis()->SetTitle("N(Vtx)");
751     gr[i]->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
752     gr[i]->GetXaxis()->CenterTitle();
753     gr[i]->GetYaxis()->CenterTitle();
754     gr[i]->SetName(name.c_str());
755     gr[i]->SetMarkerSize(0.00001);
756     gr[i]->GetYaxis()->SetRangeUser(0.90,1.10);
757     gr[i]->SetLineColor(TColor::GetColor(color.c_str()));
758     gr[i]->SetMarkerColor(TColor::GetColor(color.c_str()));
759     if(dofit) {
760     gr[i]->Fit("pol1");
761     fit[i] = (TF1*)gr[i]->GetFunction("pol1");
762     fit[i]->SetLineColor(TColor::GetColor(color.c_str()));
763     }
764     leg->AddEntry(gr[i],name.c_str(),"l");
765     if(i==0) gr[i]->Draw("AP*");
766     else gr[i]->Draw("P");
767     }
768     leg->Draw();
769     DrawPrelim();
770    
771     CompleteSave(can,"ScenarioComparisonInclusive");
772     delete can;
773    
774     }
775 buchmann 1.2
776 buchmann 1.7 void compare_selection(string identifier) {
777     bool recognized_scenario=false;
778    
779     cout << "Running with identifier " << identifier << endl;
780 buchmann 1.9 if(identifier=="Zb3010") {
781     LeadingB=TCut ("Zb3010_bTagProbCSVBP[0]>0.898");
782     EtaB=TCut("abs(Zb3010_pfJetGoodEta[0])<1.3");
783     PhiZcut=TCut("abs(Zb3010_pfJetDphiZ[0])>2.7");
784     recognized_scenario=true;
785     }
786    
787     if(identifier=="Zb3010") {
788     LeadingB=TCut ("Zb3010_bTagProbCSVBP[0]>0.898");
789     EtaB=TCut("abs(Zb3010_pfJetGoodEta[0])<1.3");
790     PhiZcut=TCut("abs(Zb3010_pfJetDphiZ[0])>2.7");
791 buchmann 1.7 recognized_scenario=true;
792     }
793    
794 buchmann 1.9 if(identifier=="Zb40") {
795     LeadingB=TCut ("Zb40_bTagProbCSVBP[0]>0.898");
796     EtaB=TCut("abs(Zb40_pfJetGoodEta[0])<1.3");
797     PhiZcut=TCut("abs(Zb40_pfJetDphiZ[0])>2.7");
798     recognized_scenario=true;
799     }
800    
801     if(identifier=="Zb301030") {
802     LeadingB=TCut ("Zb301030_bTagProbCSVBP[0]>0.898");
803     EtaB=TCut("abs(Zb301030_pfJetGoodEta[0])<1.3");
804     PhiZcut=TCut("abs(Zb301030_pfJetDphiZ[0])>2.7");
805 buchmann 1.7 recognized_scenario=true;
806     }
807    
808     if(identifier=="Zb1530") {
809 buchmann 1.9 LeadingB=TCut ("Zb1530_bTagProbCSVBP[0]>0.898");
810 buchmann 1.7 EtaB=TCut("abs(Zb1530_pfJetGoodEta[0])<1.3");
811     PhiZcut=TCut("abs(Zb1530_pfJetDphiZ[0])>2.7");
812     recognized_scenario=true;
813     }
814    
815 buchmann 1.9 if(identifier=="Zb301010") {
816     LeadingB=TCut ("Zb301010_bTagProbCSVBP[0]>0.898");
817     EtaB=TCut("abs(Zb301010_pfJetGoodEta[0])<1.3");
818     PhiZcut=TCut("abs(Zb301010_pfJetDphiZ[0])>2.7");
819     recognized_scenario=true;
820     }
821    
822     if(identifier=="Zb1510") {
823     LeadingB=TCut ("Zb1510_bTagProbCSVBP[0]>0.898");
824     EtaB=TCut("abs(Zb1510_pfJetGoodEta[0])<1.3");
825     PhiZcut=TCut("abs(Zb1510_pfJetDphiZ[0])>2.7");
826     recognized_scenario=true;
827     }
828    
829     if(identifier=="Zb301010") {
830     LeadingB=TCut ("Zb301010_bTagProbCSVBP[0]>0.898");
831     EtaB=TCut("abs(Zb301010_pfJetGoodEta[0])<1.3");
832     PhiZcut=TCut("abs(Zb301010_pfJetDphiZ[0])>2.7");
833     recognized_scenario=true;
834     }
835    
836     if(identifier=="ZbMikko") {
837     LeadingB=TCut ("ZbMikko_bTagProbCSVBP[0]>0.898 && ZbMikko__IsClean");
838     EtaB=TCut("abs(ZbMikko_pfJetGoodEta[0])<1.3 && ZbMikko__IsClean");
839     PhiZcut=TCut("abs(ZbMikko_pfJetDphiZ[0])>2.7 && ZbMikko__IsClean");
840 buchmann 1.7 recognized_scenario=true;
841     }
842    
843 buchmann 1.9 if(identifier=="Zb3010_SecEta3") {
844     LeadingB=TCut ("Zb3010_SecEta3_bTagProbCSVBP[0]>0.898");
845     EtaB=TCut("abs(Zb3010_SecEta3_pfJetGoodEta[0])<1.3");
846     PhiZcut=TCut("abs(Zb3010_SecEta3_pfJetDphiZ[0])>2.7");
847 buchmann 1.7 recognized_scenario=true;
848     }
849    
850 buchmann 1.9 if(identifier=="Zb3010_SecEta5") {
851     LeadingB=TCut ("Zb3010_SecEta5_bTagProbCSVBP[0]>0.898");
852     EtaB=TCut("abs(Zb3010_SecEta5_pfJetGoodEta[0])<1.3");
853     PhiZcut=TCut("abs(Zb3010_SecEta5_pfJetDphiZ[0])>2.7");
854 buchmann 1.7 recognized_scenario=true;
855     }
856    
857 buchmann 1.9 if(identifier=="Zb3010_p5Clean") {
858     LeadingB=TCut ("Zb3010_p5Clean_bTagProbCSVBP[0]>0.898 && Zb3010_p5Clean_IsClean");
859     EtaB=TCut("abs(Zb3010_p5Clean_pfJetGoodEta[0])<1.3 && Zb3010_p5Clean_IsClean");
860     PhiZcut=TCut("abs(Zb3010_p5Clean_pfJetDphiZ[0])>2.7 && Zb3010_p5Clean_IsClean");
861 buchmann 1.7 recognized_scenario=true;
862     }
863    
864 buchmann 1.9 if(identifier=="Zb3010CHS") {
865     LeadingB=TCut ("Zb3010CHS_bTagProbCSVBP[0]>0.898");
866     EtaB=TCut("abs(Zb3010CHS_pfJetGoodEta[0])<1.3");
867     PhiZcut=TCut("abs(Zb3010CHS_pfJetDphiZ[0])>2.7");
868 buchmann 1.7 recognized_scenario=true;
869     }
870    
871     assert(recognized_scenario);
872    
873     cout << "Cuts: " << endl;
874     cout << " " << (const char*) LeadingB << endl;
875     cout << " " << (const char*) EtaB << endl;
876     cout << " " << (const char*) PhiZcut << endl;
877    
878     cout << endl << endl;
879    
880 buchmann 1.9 // ScenarioComparison();
881     // ScenarioComparisonInclusive();
882    
883    
884     // DoPUStudy(identifier);
885    
886     // SpecialCutFlow(identifier);
887     // draw_kin_variable("pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",17000,30,200,"Z p_{T}","Official/"+identifier+"/Zpt__"+identifier,1);
888     // draw_kin_variable(identifier+"_pfJetGoodPt[0]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",40,0,200,"Leading Jet Pt (B)","Official/"+identifier+"/LeadingJetPt__"+identifier,1);
889     // draw_kin_variable(identifier+"_pfJetGoodPt[1]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",40,0,200,"Sub-Leading Jet Pt [GeV]","Official/"+identifier+"/Jet2Pt__"+identifier,1);
890     // draw_kin_variable(identifier+"_pfJetGoodPt[1]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,0,2,"p_{T}^{J2} / p_{T}^{Z}","Official/"+identifier+"/SubLeadingJetPt_Over_ZPt__"+identifier,1);
891     draw_kin_variable(identifier+"_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<50"),20,0,2,"#alpha","Official/"+identifier+"/alpha_pt_30_to_50__"+identifier,1);
892     draw_kin_variable(identifier+"_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>50&&pt<75"),20,0,2,"#alpha","Official/"+identifier+"/alpha_50_to_75__"+identifier,1);
893     draw_kin_variable(identifier+"_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>75&&pt<125"),20,0,2,"#alpha","Official/"+identifier+"/alpha_pt_75_to_125__"+identifier,1);
894     draw_kin_variable(identifier+"_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>125&&pt<1000"),20,0,2,"#alpha","Official/"+identifier+"/alpha_pt_125_to_1000__"+identifier,1);
895     draw_kin_variable(identifier+"_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut,20,0,2,"#alpha","Official/"+identifier+"/alpha_full__"+identifier,1);
896     // draw_kin_variable(identifier+"_pfBJetDphiZ[0]",ZplusBsel&&LeadingB&&"pt>10&&pt<1000&&pfBJetDphiZ[0]>0",30,0,3.2,"#delta#phi (Z,b lead)","Official/"+identifier+"/DeltaPhiZBlead__"+identifier,0);
897 buchmann 1.7
898 buchmann 1.9 // draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.3").c_str()),20,0,2,"Z+b MPF","Official/"+identifier+"/mpf_alpha_smaller_0p3___"+identifier,0,"#alpha<0.3, 10 GeV < p_{T}^{Z} < 1000 GeV");
899     // draw_kin_variable(identifier+"_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.3").c_str()),20,0,2,"p_{T} b-jet / p_{T} Z","Official/"+identifier+"/ptb_over_ptz___alpha_smaller_0p3______"+identifier,0,"#alpha<0.3, 10 GeV < p_{T}^{Z} < 1000 GeV");
900 buchmann 1.7
901    
902 buchmann 1.9 // */
903 buchmann 1.7
904     }
905    
906 buchmann 1.12 void do_basic_ZB_analysis(int do_inclusive) {
907    
908     //https://twiki.cern.ch/twiki/bin/viewauth/CMS/BtagPOG#2011_Data_and_MC
909     //https://twiki.cern.ch/twiki/pub/CMS/BtagPOG/SFb-ttbar_payload.txt
910 buchmann 1.9
911 buchmann 1.12 if(do_inclusive) {
912     LeadingB=TCut("Zb3010_bTagProbCSVBP[0]>=-10000");
913     cutWeight=TCut("(weight*(weight<1000)*(is_data+(!is_data)*((id1==id2&&id1==0)*0.95+(id1==id2&&id1==1)*0.88+(id1!=id2)*0.92)))");
914     write_info(__FUNCTION__,"Inclusive selection - no btag requirement, no btag efficiency/mistag correction (weight)");
915     } else {
916     LeadingB=TCut("Zb3010_bTagProbCSVBP[0]>0.898");
917     cutWeight=TCut("(weight*(weight<1000)*(is_data+(!is_data)*((id1==id2&&id1==0)*0.95+(id1==id2&&id1==1)*0.88+(id1!=id2)*0.92)*Zb3010_BTagWgt))");
918     write_info(__FUNCTION__,"Exclusive selection - switched on btag requirement and btag efficiency/mistag correction (weight)");
919     }
920    
921     bool doquick=true;
922 buchmann 1.9
923 buchmann 1.2 TCanvas *zbcanvas = new TCanvas("zbcanvas","zbcanvas");
924 buchmann 1.1
925 buchmann 1.9 // compare_selection("Zb3010_p5Clean");
926     // compare_selection("Zb3010CHS");
927     // compare_selection("Zb3010");
928     // compare_selection("Zb301030");
929     // compare_selection("Zb1530");
930     // compare_selection("Zb301010");
931     // compare_selection("Zb3010_SecEta3");
932     // compare_selection("Zb3010_SecEta5");
933     // compare_selection("Zb1510");
934     // compare_selection("Zb301010");
935     // compare_selection("ZbMikko");
936     // compare_selection("Zb3010");
937     // compare_selection("Zb40");
938 buchmann 1.8
939 buchmann 1.11 if(!doquick) {
940     print_all_b_yields();
941     draw_mpf_vars();
942     // DrawEvilCutFlow();
943    
944     draw_Zb_kin_vars();
945    
946     }
947 buchmann 1.7
948 buchmann 1.9 new_data_mc_agreement_2d();
949 buchmann 1.1
950 buchmann 1.12 write_error(__FUNCTION__,"Need to implement switch for inclusive result; inclusive result needs to have btag efficiency reweighting deactivated!");
951 buchmann 1.2 delete zbcanvas;
952 buchmann 1.1 }
953 buchmann 1.7
954    
955    
956    
957     //*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/* DELETE EVERYTHING BELOW THIS LINE /*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*//
958    
959    
960     const char* concatenate(string bla1, string bla2) {
961     stringstream bla;
962     bla << bla1 << bla2;
963     return bla.str().c_str();
964     }
965    
966     void SpecialCutFlow(string identifier) {
967     stringstream MegaCut;
968    
969    
970    
971     MegaCut << " 1*(" << (const char*) (ZplusBsel&&TCut(concatenate(identifier,"_pfJetGoodNumBtag>0"))) << ")";
972     MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB) << ")";
973     MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB) << ")";
974     MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut) << ")";
975 buchmann 1.9 MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000").c_str())) << ")";
976     MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.3").c_str())) << ")";
977     MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.2").c_str())) << ")";
978     MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.15").c_str())) << ")";
979     MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.1").c_str())) << ")";
980     MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.05").c_str())) << ")";
981 buchmann 1.7 MegaCut << " ";
982    
983     TCanvas *can = new TCanvas("can","can");
984     can->SetLogy(1);
985 buchmann 1.9 TCut basecut(ZplusBsel);
986 buchmann 1.7
987    
988    
989     TLegend *leg = allsamples.allbglegend();
990    
991    
992     TH1F *data = allsamples.Draw("data", MegaCut.str(),11,-0.5,10.5, "", "events", basecut,data,luminosity);
993     THStack mcm = allsamples.DrawStack("mc", MegaCut.str(),11,-0.5,10.5, "", "events", basecut,mc,luminosity);
994    
995     float runningsum=0;
996     for(int i=data->GetNbinsX();i>0;i--) {
997     runningsum+=data->GetBinContent(i);
998     data->SetBinContent(i,runningsum);
999     }
1000    
1001     float minimum=data->GetMaximum();
1002    
1003     TH1F* h;
1004     TIter nextOF(mcm.GetHists());
1005    
1006     float mcsum=0;
1007     float zb=0;
1008     while ( h = (TH1F*)nextOF() ) {
1009     float runningsum=0;
1010     minimum=data->GetMaximum();//done deliberately at each step so get the minimum of the last (!) contribution
1011     for(int i=h->GetNbinsX();i>0;i--) {
1012     runningsum+=h->GetBinContent(i);
1013     h->SetBinContent(i,runningsum);
1014     if(runningsum<minimum&&runningsum>0) minimum=runningsum;
1015     }
1016     if(Contains(h->GetName(),"Z_b_")) {
1017     zb=h->GetBinContent(h->GetNbinsX());
1018     }
1019     mcsum+=h->GetBinContent(h->GetNbinsX());
1020     }
1021    
1022    
1023    
1024     data->SetMinimum(0.05*minimum);
1025     can->cd(1)->SetBottomMargin(0.25);
1026     data->GetXaxis()->SetBinLabel(1,"Pre-selection");
1027     data->GetXaxis()->SetBinLabel(2,"Contains b jet");
1028     data->GetXaxis()->SetBinLabel(3,"Leading jet is b-jet");
1029     data->GetXaxis()->SetBinLabel(4,"|#eta_{b}|<1.3 ");
1030     data->GetXaxis()->SetBinLabel(5,"|#delta#phi(b,Z)|>2.7");
1031     data->GetXaxis()->SetBinLabel(6,"10<p_{T}^{Z}<1000 GeV");
1032     data->GetXaxis()->SetBinLabel(7,"#alpha < 0.3");
1033     data->GetXaxis()->SetBinLabel(8,"#alpha < 0.2");
1034     data->GetXaxis()->SetBinLabel(9,"#alpha < 0.15");
1035     data->GetXaxis()->SetBinLabel(10,"#alpha < 0.1");
1036     data->GetXaxis()->SetBinLabel(11,"#alpha < 0.05");
1037     data->GetXaxis()->LabelsOption("v");
1038    
1039     stringstream purityinfo;
1040     purityinfo << "Purity: " << std::setprecision(3) << 100*zb/mcsum << " %";
1041     stringstream neventsinfo;
1042     neventsinfo << "Nevents: " << data->GetBinContent(data->GetNbinsX());
1043     TH1F *crap = new TH1F("crap","",1,0,1);
1044     crap->SetLineColor(kWhite);
1045     leg->AddEntry(crap,purityinfo.str().c_str(),"l");
1046     leg->AddEntry(crap,neventsinfo.str().c_str(),"l");
1047    
1048     data->Draw("e1");
1049     mcm.Draw("same");
1050     data->Draw("same,e1,axis");
1051     data->Draw("same,e1");
1052     // data->Draw("same,TEXT");
1053    
1054     leg->Draw();
1055    
1056     DrawPrelim();
1057    
1058     CompleteSave(can,"CutFlow___"+identifier);
1059    
1060     delete crap;
1061     delete data;
1062    
1063    
1064     }