ra4_draw  4bd0201e3d922d42bd545d4b045ed44db33454a4
table_2016-06-03_freeze_preapp_an.cxx
Go to the documentation of this file.
1 #include <fstream>
2 #include <iostream>
3 #include <vector>
4 #include <ctime>
5 #include <iomanip> // setw
6 
7 #include <unistd.h>
8 #include <stdlib.h>
9 #include <getopt.h>
10 
11 #include "TError.h" // Controls error level reporting
12 #include "RooStats/NumberCountingUtils.h"
13 
14 #include "utilities.hpp"
15 #include "baby.hpp"
16 #include "process.hpp"
17 #include "named_func.hpp"
18 #include "plot_maker.hpp"
19 #include "palette.hpp"
20 #include "table.hpp"
21 using namespace std;
22 
23 namespace{
24  TString method = "met200";
25  double lumi(0.815);
26  bool do_other(false);
27  bool full_lumi(false);
28  float syst = 0.001;
29  TString blind_s = "$\\spadesuit$";
30  bool really_unblind = false;
31  bool single_thread = false;
32 }
33 
34 template<typename T>
35 shared_ptr<Process> Proc(const string process_name, Process::Type type,
36  int color, const set<string> &files, const string &cut = "1"){
37  return make_shared<Process>(process_name, type, color,
38  unique_ptr<Baby>(new T(files)),
39  cut);
40 }
41 
42 void GetOptions(int argc, char *argv[]);
43 TString Zbi(double Nobs, double Nbkg, double Ebkg){
44  double Nsig = Nobs-Nbkg;
45  double zbi = RooStats::NumberCountingUtils::BinomialExpZ(Nsig, Nbkg, Ebkg/Nbkg);
46  if(Nbkg==0) zbi = RooStats::NumberCountingUtils::BinomialWithTauExpZ(Nsig, Nbkg, 1/Ebkg);
47  if(zbi<0) zbi=0;
48  TString zbi_s = "$"+RoundNumber(zbi,1)+"\\sigma$";
49  if(Nsig<=0) zbi_s = "-";
50  //cout<<"Zbi for Nobs "<<Nobs<<", Nbkg "<<Nbkg<<", Ebkg "<<Ebkg<<" is "<<zbi_s<<endl;
51  return zbi_s;
52 }
53 
54 int main(int argc, char *argv[]){
55  gErrorIgnoreLevel=6000; // Turns off ROOT errors due to missing branches
56  GetOptions(argc, argv);
57 
58  time_t begtime, endtime;
59  time(&begtime);
60 
62  string bfolder("");
63  string hostname = execute("echo $HOSTNAME");
64  if(Contains(hostname, "cms") || Contains(hostname, "compute-"))
65  bfolder = "/net/cms2"; // In laptops, you can't create a /net folder
66 
67  string foldermc(bfolder+"/cms2r0/babymaker/babies/2016_06_14/mc/merged_standard/");
68  string folderdata(bfolder+"/cms2r0/babymaker/babies/2016_06_21/data/skim_standard/");
69  if(method.Contains("met150")){
70  foldermc = bfolder+"/cms2r0/babymaker/babies/2016_06_14/mc/merged_1lht500met150nj5/";
71  folderdata = bfolder+"/cms2r0/babymaker/babies/2016_06_21/data/skim_1lmet150/";
72  //folderdata = bfolder+"/cms2r0/babymaker/babies/2016_06_14/data/merged_1lht500met150nj5/";
73  }
74  if(method.Contains("2015")){
75  foldermc = bfolder+"/cms2r0/babymaker/babies/2016_04_29/mc/merged_1lht500met200/";
76  folderdata = bfolder+"/cms2r0/babymaker/babies/2016_04_29/data/merged_1lht500met200/";
77  }
78 
79  Palette colors("txt/colors.txt", "default");
80 
81  auto tt = Proc<Baby_full>("t#bar{t}", Process::Type::background, colors("tt_1l"),
82  {foldermc+"*_TTJets*Lept*.root", foldermc+"*_TTJets_HT*.root",
83  foldermc+"*_WJetsToLNu*.root",foldermc+"*_ST_*.root",
84  foldermc+"*_TTW*.root",foldermc+"*_TTZ*.root",
85  foldermc+"*DYJetsToLL*.root",foldermc+"*QCD_HT*.root",
86  foldermc+"*_ZJet*.root",foldermc+"*_ttHJetTobb*.root",
87  foldermc+"*_TTGJets*.root",foldermc+"*_TTTT*.root",
88  foldermc+"*_WH_HToBB*.root",foldermc+"*_ZH_HToBB*.root",
89  foldermc+"*_WWTo*.root",foldermc+"*_WZ*.root",foldermc+"*_ZZ_*.root"},
90  "stitch");
91  // auto tt = Proc<Baby_full>("t#bar{t}", Process::Type::background, colors("tt_1l"),
92  // {foldermc+"*T1tttt_mGluino-1500_mLSP-100_*.root"},
93  // "stitch");
94  auto other = Proc<Baby_full>("Other", Process::Type::background, colors("other"),
95  {foldermc+"*_WJetsToLNu*",foldermc+"*_ST_*.root",
96  foldermc+"*_TTW*.root",foldermc+"*_TTZ*.root",
97  foldermc+"*DYJetsToLL*.root",foldermc+"*QCD_HT*.root",
98  foldermc+"*_ZJet*.root",foldermc+"*_ttHJetTobb*.root",
99  foldermc+"*_TTGJets*.root",foldermc+"*_TTTT*.root",
100  foldermc+"*_WH_HToBB*.root",foldermc+"*_ZH_HToBB*.root",
101  foldermc+"*_WWTo*.root",foldermc+"*_WZ*.root",foldermc+"*_ZZ_*.root"});
102 
103  string trigs = "(trig[4]||trig[8]||trig[13]||trig[33])";
104  if(method.Contains("2015")) trigs = "(trig[4]||trig[8]||trig[28]||trig[14])";
105  auto data = Proc<Baby_full>("Data", Process::Type::data, kBlack,
106  {folderdata+"*.root"},trigs);
107 
108  vector<shared_ptr<Process> > all_procs = {data, tt};
109  if (do_other) all_procs.push_back(other);
110 
111 
113  TString base_s = "mj14>250&&njets>=5&&stitch&&pass&&nonblind";
114 
115  vector<TString> njbcuts_stdnob = {"njets>=6&&njets<=8", "njets>=9", "njets>=6&&njets<=8", "njets>=9"};
116  vector<TString> njbcuts_2l = {"njets>=5&&njets<=7", "njets>=8", "njets>=5&&njets<=7", "njets>=8"};
117  vector<TString> njbcuts_5j = {"nbm==1&&njets==5", "nbm>=2&&njets==5", "nbm==1&&njets==5", "nbm>=2&&njets==5"};
118  vector<TString> njbcuts_m1lmet150nb12 = {"nbm==1&&njets>=6&&njets<=8", "nbm==1&&njets>=9",
119  "nbm>=2&&njets>=6&&njets<=8", "nbm>=2&&njets>=9"};
120  vector<TString> njbcuts_std = {"nbm==1&&njets>=6&&njets<=8", "nbm==1&&njets>=9",
121  "nbm==2&&njets>=6&&njets<=8", "nbm==2&&njets>=9",
122  "nbm>=3&&njets>=6&&njets<=8", "nbm>=3&&njets>=9",
123  "nbm==1&&njets>=6&&njets<=8", "nbm==1&&njets>=9",
124  "nbm==2&&njets>=6&&njets<=8", "nbm==2&&njets>=9",
125  "nbm>=3&&njets>=6&&njets<=8", "nbm>=3&&njets>=9"};
126  vector<TString> njbcuts_nb1 = {"nbm==1&&njets>=6&&njets<=8", "nbm==1&&njets>=9",
127  "nbm==2&&njets>=6&&njets<=8", "nbm==2&&njets>=9",
128  "nbm>=3&&njets>=6&&njets<=8", "nbm>=3&&njets>=9"};
129  vector<TString> njbcuts_met500 = {"nbm==1&&njets>=6&&njets<=8", "nbm==1&&njets>=9",
130  "nbm==2&&njets>=6&&njets<=8", "nbm==2&&njets>=9",
131  "nbm>=3&&njets>=6&&njets<=8", "nbm>=3&&njets>=9"};
132 
133  vector<TString> njbcuts_m1lnob = {"njets>=6&&njets<=8", "njets>=9"};
134  vector<TString> njbcuts_m2lnob = {"njets>=5&&njets<=7", "njets>=8"};
135  vector<TString> njbcuts_2lveto_lonj = {"njets>=5", "njets>=5"};
136  vector<TString> njbcuts_2lveto_hinj = {"njets>=8", "njets>=8"};
137  vector<TString> njbcuts_all = {"njets>=6"};
138  vector<TString> njbcuts_2l_lonj = {"njets>=5"};
139  vector<TString> njbcuts_2l_hinj = {"njets>=8"};
140 
141 
142  size_t ilowmet(2); // njbcuts index up to which metcuts[0] is applied
143  vector<TString> metcuts = {"met<=350", "met>350&&met<=500"};
144 
145  vector<TString> abcdcuts_std = {"mt<=140&&mj14<=400",
146  "mt<=140&&mj14>400",
147  "mt>140&&mj14<=400",
148  "mt>140&&mj14>400"};
149  vector<TString> abcdcuts_2lveto = {"mt<=140&&mj14<=400&&nleps==1&&nveto==0&&nbm>=1&&njets>=6",
150  "mt<=140&&mj14>400&&nleps==1&&nveto==0&&nbm>=1&&njets>=6",
151  "mj14<=400&&(nleps==2&&nbm<=2&&njets>=5 || nleps==1&&nveto==1&&nbm>=1&&nbm<=2&&mt>140&&njets>=6)",
152  "mj14> 400&&(nleps==2&&nbm<=2&&njets>=5 || nleps==1&&nveto==1&&nbm>=1&&nbm<=2&&mt>140&&njets>=6)"};
153 
154  vector<TString> abcdcuts_2lveto_el = {"mt<=140&&mj14<=400&&nels==1&&nveto==0&&nbm>=1&&njets>=6",
155  "mt<=140&&mj14>400&&nels==1&&nveto==0&&nbm>=1&&njets>=6",
156  "mj14<=400&&(nels==2&&nbm<=2&&njets>=5 || nels==1&&nveto==1&&nbm>=1&&nbm<=2&&mt>140&&njets>=6)",
157  "mj14> 400&&(nels==2&&nbm<=2&&njets>=5 || nels==1&&nveto==1&&nbm>=1&&nbm<=2&&mt>140&&njets>=6)"};
158 
159  vector<TString> abcdcuts_2lveto_mu = {"mt<=140&&mj14<=400&&nmus==1&&nveto==0&&nbm>=1&&njets>=6",
160  "mt<=140&&mj14>400&&nmus==1&&nveto==0&&nbm>=1&&njets>=6",
161  "mj14<=400&&(nmus==2&&nbm<=2&&njets>=5 || nmus==1&&nveto==1&&nbm>=1&&nbm<=2&&mt>140&&njets>=6)",
162  "mj14> 400&&(nmus==2&&nbm<=2&&njets>=5 || nmus==1&&nveto==1&&nbm>=1&&nbm<=2&&mt>140&&njets>=6)"};
163 
164  vector<TString> abcdcuts_2lveto_lonj = {"mt<=140&&mj14<=400&&nleps==1&&nveto==0&&nbm>=1&&njets>=6",
165  "mt<=140&&mj14>400&&nleps==1&&nveto==0&&nbm>=1&&njets>=6&&njets<=8",
166  "mj14<=400&&(nleps==2&&nbm<=2&&njets>=5 || nleps==1&&nveto==1&&nbm>=1&&nbm<=2&&mt>140&&njets>=6)",
167  "mj14>400&&(nleps==2&&nbm<=2&&njets>=5&&njets<=7 || nleps==1&&nveto==1&&nbm>=1&&nbm<=2&&mt>140&&njets>=6&&njets<=8)"};
168 
169 
170  vector<TString> abcdcuts_2lveto_hinj = {"mt<=140&&mj14<=400&&nleps==1&&nveto==0&&nbm>=1&&njets>=6",
171  "mt<=140&&mj14>400&&nleps==1&&nveto==0&&nbm>=1&&njets>=9",
172  "mj14<=400&&(nleps==2&&nbm<=2&&njets>=5 || nleps==1&&nveto==1&&nbm>=1&&nbm<=2&&mt>140&&njets>=6)",
173  "mj14>400&&(nleps==2&&nbm<=2&&njets>=8 || nleps==1&&nveto==1&&nbm>=1&&nbm<=2&&mt>140&&njets>=9)"};
174 
175  vector<TString> abcdcuts_veto = {"mt<=140&&mj14<=400&&njets>=6&&nbm>=1&&nleps==1&&nveto==0",
176  "mt<=140&&mj14>400&&nbm>=1&&nleps==1&&nveto==0",
177  "mt>140&&mj14<=400&&njets>=6&&nbm>=1&&nbm<=2&&nleps==1&&nveto==1",
178  "mt>140&&mj14>400&&njets>=6&&nbm>=1&&nbm<=2&&nleps==1&&nveto==1"};
179  vector<TString> abcdcuts_2l = {"mt<=140&&mj14<=400&&njets>=6&&nbm>=1&&nleps==1&&nveto==0",
180  "mt<=140&&mj14>400&&nbm>=1&&nleps==1&&nveto==0",
181  "mj14<=400&&nbm<=2&&nleps==2",
182  "mj14>400&&nbm<=2&&nleps==2"};
183 
184 
185  vector<TString> abcdcuts, njbcuts_himt = njbcuts_stdnob;
186  vector<TString> njbcuts = njbcuts_stdnob;
187  TString region_s = "R", method_s, base_all = "mj14>250&&pass&&nonblind&&stitch&&";
188  TString lumi_s = "815$ pb$^{-1}$";
189  bool unblind = false;
190 
191  if(full_lumi){
192  base_all = "mj14>250&&pass&&stitch&&";
193  if(method.Contains("2015")){
194  lumi = 2.3;
195  lumi_s = "2.3$ fb$^{-1}$";
196  } else {
197  lumi = 2.6;
198  lumi_s = "2.6$ fb$^{-1}$";
199  }
200  }
201 
202 
203  if(method=="m2l" || method=="m2l_2015") {
204  base_s = base_all+"njets>=5";
205  njbcuts_himt = njbcuts_2l;
206  abcdcuts = abcdcuts_2l;
207  region_s = "D";
208  method_s = "D3 and D4 have $2\\ell$";
209  } else if(method=="m2lveto" || method=="m2lveto_2015") {
210  base_s = base_all+"njets>=5";
211  njbcuts = njbcuts_2lveto_lonj;
212  njbcuts_himt = njbcuts_2lveto_lonj;
213  abcdcuts = abcdcuts_2lveto;
214  region_s = "D";
215  method_s = "D3 and D4 have $\\ell\\ell+\\ell t_{\\rm veto}$";
216  ilowmet = 1;
217  } else if(method=="m2lvetomet150" || method=="m2lvetomet150_2015") {
218  base_s = base_all+"njets>=5";
219  njbcuts = njbcuts_2lveto_lonj;
220  njbcuts_himt = njbcuts_2lveto_lonj;
221  abcdcuts = abcdcuts_2lveto;
222  region_s = "D";
223  method_s = "$150<E_{\\rm T}^{\\rm miss}<200$, D3 and D4 have $\\ell\\ell+\\ell t_{\\rm veto}$";
224  ilowmet = njbcuts.size();
225  } else if(method=="m2lveto_el" || method=="m2lveto_el_2015") {
226  base_s = base_all+"njets>=5&&nels>=1";
227  njbcuts = njbcuts_2lveto_lonj;
228  njbcuts_himt = njbcuts_2lveto_lonj;
229  abcdcuts = abcdcuts_2lveto_el;
230  region_s = "D";
231  method_s = "D3 and D4 have $ee+e t_{\\rm veto}$";
232  ilowmet = 1;
233  } else if(method=="m2lveto_mu" || method=="m2lveto_mu_2015") {
234  base_s = base_all+"njets>=5&&nmus>=1";
235  njbcuts = njbcuts_2lveto_lonj;
236  njbcuts_himt = njbcuts_2lveto_lonj;
237  abcdcuts = abcdcuts_2lveto_mu;
238  region_s = "D";
239  method_s = "D3 and D4 have $\\mu\\mu+\\mu t_{\\rm veto}$";
240  ilowmet = 1;
241  } else if(method=="m2lveto_lonj") {
242  base_s = base_all+"njets>=5";
243  njbcuts = njbcuts_2lveto_lonj;
244  njbcuts_himt = njbcuts_2lveto_lonj;
245  abcdcuts = abcdcuts_2lveto_lonj;
246  region_s = "D";
247  method_s = "D3 and D4 have $\\ell\\ell+\\ell t_{\\rm veto}$ - low $N_{\\rm jets}$";
248  ilowmet = 1;
249  } else if(method=="m2lveto_hinj") {
250  base_s = base_all+"njets>=5";
251  njbcuts = njbcuts_2lveto_hinj;
252  njbcuts_himt = njbcuts_2lveto_hinj;
253  abcdcuts = abcdcuts_2lveto_hinj;
254  region_s = "D";
255  method_s = "D3 and D4 have $\\ell\\ell+\\ell t_{\\rm veto}$ - high $N_{\\rm jets}$";
256  ilowmet = 1;
257  } else if(method=="m2lvetomet150_lonj") {
258  base_s = base_all+"njets>=5";
259  njbcuts = njbcuts_2l_lonj;
260  njbcuts_himt = njbcuts_2l_lonj;
261  abcdcuts = abcdcuts_2lveto_lonj;
262  region_s = "D";
263  method_s = "D3 and D4 have $\\ell\\ell+\\ell t_{\\rm veto}$ - low $N_{\\rm jets}$";
264  ilowmet = njbcuts.size();
265  } else if(method=="m2lvetomet150_hinj") {
266  base_s = base_all+"njets>=5";
267  njbcuts = njbcuts_2l_hinj;
268  njbcuts_himt = njbcuts_2l_hinj;
269  abcdcuts = abcdcuts_2lveto_hinj;
270  region_s = "D";
271  method_s = "D3 and D4 have $\\ell\\ell+\\ell t_{\\rm veto}$ - high $N_{\\rm jets}$";
272  ilowmet = njbcuts.size();
273  } else if(method=="mveto" || method=="mveto_2015") {
274  base_s = base_all+"njets>=6&&nbm>=1&&nleps==1";
275  njbcuts_himt = njbcuts_stdnob;
276  abcdcuts = abcdcuts_veto;
277  region_s = "D";
278  method_s = "D3 and D4 have $\\ell t_{\\rm veto}$";
279  } else if(method=="m5j") {
280  base_s = base_all+"njets==5&&nbm>=1&&nleps==1&&nveto==0";
281  njbcuts = njbcuts_5j;
282  njbcuts_himt = njbcuts_5j;
283  abcdcuts = abcdcuts_std;
284  method_s = "$N_{\\rm jets}=5$";
285  } else if(method=="m1lmet150") {
286  base_s = base_all+"njets>=6&&nbm>=1&&nleps==1&&nveto==0";
287  njbcuts = njbcuts_m1lmet150nb12;
288  njbcuts_himt = njbcuts_m1lmet150nb12;
289  abcdcuts = abcdcuts_std;
290  method_s = "$1\\ell, 150<E_{\\rm T}^{\\rm miss}<200$";
291  ilowmet = njbcuts.size();
292  } else if(method=="mvetomet150" || method=="mvetomet150_2015") {
293  base_s = base_all+"njets>=6&&nbm>=1&&nleps==1";
294  njbcuts = njbcuts_m1lnob;
295  njbcuts_himt = njbcuts_m1lnob;
296  abcdcuts = abcdcuts_veto;
297  region_s = "D";
298  method_s = "$150<E_{\\rm T}^{\\rm miss}<200$, D3 and D4 have $\\ell t_{\\rm veto}$";
299  } else if(method=="m2lmet150" || method=="m2lmet150_2015") {
300  base_s = base_all+"njets>=5";
301  njbcuts = njbcuts_m1lnob;
302  njbcuts_himt = njbcuts_m2lnob;
303  abcdcuts = abcdcuts_2l;
304  region_s = "D";
305  method_s = "$150<E_{\\rm T}^{\\rm miss}<200$, D3 and D4 have $\\ell\\ell+\\ell t_{\\rm veto}$";
306  } else if(method=="mall") {
307  base_s = base_all+"njets>=6&&nbm>=1&&met>200&&nleps==1&&nveto==0";
308  njbcuts = njbcuts_all;
309  njbcuts_himt = njbcuts_all;
310  abcdcuts = abcdcuts_std;
311  method_s = "All signal bins - $1\\ell$, $E_{\\rm T}^{\\rm miss}>200$";
312  metcuts[0] = "met>200";
313  ilowmet = njbcuts.size();
314  } else if(method=="met200") {
315  base_s = base_all+"njets>=6&&nbm>=1&&met>200&&nleps==1&&nveto==0";
316  njbcuts = njbcuts_std;
317  njbcuts_himt = njbcuts_std;
318  abcdcuts = abcdcuts_std;
319  method_s = "Signal bins - $1\\ell$, $200<E_{\\rm T}^{\\rm miss}\\leq500$";
320  ilowmet = 6;
321  } else if(method=="met500") {
322  base_s = base_all+"njets>=6&&nbm>=1&&met>500&&nleps==1&&nveto==0";
323  njbcuts = njbcuts_met500;
324  njbcuts_himt = njbcuts_met500;
325  abcdcuts = abcdcuts_std;
326  method_s = "Signal bins - $1\\ell$, $E_{\\rm T}^{\\rm miss}>500$";
327  metcuts[0] = "met>500";
328  ilowmet = njbcuts.size();
329  } else if(method=="met200nb1") {
330  base_s = base_all+"njets>=6&&nbm>=1&&met>200&&nleps==1&&nveto==0";
331  njbcuts = njbcuts_nb1;
332  njbcuts_himt = njbcuts_nb1;
333  abcdcuts = abcdcuts_std;
334  method_s = "Signal bins - $1\\ell$, separated N_{b}=1$, $200<E_{\\rm T}^{\\rm miss}\\leq350$";
335  metcuts[0] = "met>200&&met<=350&&nbm==1";
336  metcuts[1] = "met>200&&met<=350&&nbm>=2";
337  ilowmet = 2;
338  } else if(method=="met350nb1") {
339  base_s = base_all+"njets>=6&&nbm>=1&&met>350&&met<=500&&nleps==1&&nveto==0";
340  njbcuts = njbcuts_nb1;
341  njbcuts_himt = njbcuts_nb1;
342  abcdcuts = abcdcuts_std;
343  method_s = "Signal bins - $1\\ell$, separated N_{b}=1$, $350<E_{\\rm T}^{\\rm miss}\\leq500$";
344  metcuts[0] = "met>350&&met<=500&&nbm==1";
345  metcuts[1] = "met>350&&met<=500&&nbm>=2";
346  ilowmet = 2;
347  } else if(method=="met500nb1") {
348  base_s = base_all+"njets>=6&&nbm>=1&&met>500&&nleps==1&&nveto==0";
349  njbcuts = njbcuts_nb1;
350  njbcuts_himt = njbcuts_nb1;
351  abcdcuts = abcdcuts_std;
352  method_s = "Signal bins - $1\\ell$, separated N_{b}=1$, $E_{\\rm T}^{\\rm miss}>500$";
353  metcuts[0] = "met>500&&nbm==1";
354  metcuts[1] = "met>500&&nbm>=2";
355  ilowmet = 2;
356  }else if(method=="agg_himet"){
357  base_s = base_all+"met>500&&njets>=6&&nbm>=1";
358  njbcuts = vector<TString>{"nbm>=3&&njets>=6"};
359  njbcuts_himt = njbcuts;
360  abcdcuts = abcdcuts_std;
361  method_s = "Agg. Bin: $1\\ell$, MET500, $N_{j}\\geq6$, $N_{b}\\geq3$";
362  metcuts = vector<TString>{"met>500"};
363  ilowmet = 1;
364  }else if(method=="agg_mixed"){
365  base_s = base_all+"met>350&&njets>=6&&nbm>=1";
366  njbcuts = vector<TString>{"nbm>=2&&njets>=9"};
367  njbcuts_himt = njbcuts;
368  abcdcuts = abcdcuts_std;
369  method_s = "Agg. Bin: $1\\ell$, MET350, $N_{j}\\geq9$, $N_{b}\\geq2$";
370  metcuts = vector<TString>{"met>350"};
371  ilowmet = 1;
372  }else if(method=="agg_himult"){
373  base_s = base_all+"met>200&&njets>=6&&nbm>=1";
374  njbcuts = vector<TString>{"nbm>=3&&njets>=9"};
375  njbcuts_himt = njbcuts;
376  abcdcuts = abcdcuts_std;
377  method_s = "Agg. Bin: $1\\ell$, MET200, $N_{j}\\geq9$, $N_{b}\\geq3$";
378  metcuts = vector<TString>{"met>200"};
379  ilowmet = 1;
380  }else if(method=="agg_1b"){
381  base_s = base_all+"met>500&&njets>=6&&nbm>=1";
382  njbcuts = vector<TString>{"nbm>=1&&njets>=9"};
383  njbcuts_himt = njbcuts;
384  abcdcuts = abcdcuts_std;
385  method_s = "Agg. Bin: $1\\ell$, MET500, $N_{j}\\geq9$, $N_{b}\\geq1$";
386  metcuts = vector<TString>{"met>500"};
387  ilowmet = 1;
388  }else {
389  cout<<"Method "<<method<<" not available. Exiting"<<endl<<endl;
390  return 0;
391  }
392  if(method.Contains("2015")) method_s = "2015 data: "+method_s;
393 
394  if(really_unblind) unblind = true;
395 
397  vector<TableRow> bincuts;
398  for(size_t ind(0); ind<njbcuts.size(); ind++){
399  cout<<endl<<"New njb"<<endl;
400  for(size_t obs(0); obs < abcdcuts.size(); obs++){
401  TString totcut(abcdcuts[obs]+"&&"+metcuts[ind>=ilowmet]);
402  if(obs == 1) totcut += ("&&"+njbcuts[ind]);
403  if(obs == 3) totcut += ("&&"+njbcuts_himt[ind]);
404  cout<<base_s+"&&"+totcut<<endl;
405  bincuts.push_back(TableRow((base_s+"&&"+totcut).Data(),(base_s+"&&"+totcut).Data()));
406  } // Loop over observables going into kappa
407  } // Loop over signal bins
408 
409  PlotMaker pm;
410  pm.Push<Table>("preds", bincuts, all_procs);
411  if(single_thread) pm.multithreaded_ = false;
412  pm.MakePlots(lumi);
413 
414  Table * yield_table = static_cast<Table*>(pm.Figures().back().get());
415  vector<GammaParams> mcyield = yield_table->BackgroundYield(lumi);
416  vector<GammaParams> datayield = yield_table->DataYield(1);
417  vector<GammaParams> otheryield;
418  if (do_other) otheryield = yield_table->Yield(other, lumi);
419 
420 
421  vector<float> pow_pred;
422  pow_pred.push_back(-1); // mt<=140 mj14<=400 R1
423  pow_pred.push_back(1); // mt<=140 mj14>400 R2
424  pow_pred.push_back(1); // mt>140 mj14<=400 D3
425  vector<float> pow_tot;
426  pow_tot.push_back(-1); // mt<=140 mj14<=400 R1
427  pow_tot.push_back(1); // mt<=140 mj14>400 R2
428  pow_tot.push_back(1); // mt>140 mj14<=400 D3
429  pow_tot.push_back(1); // mt<=140 mj14<=400 R1
430  pow_tot.push_back(-1); // mt<=140 mj14>400 R2
431  pow_tot.push_back(-1); // mt>140 mj14<=400 D3
432  pow_tot.push_back(1); // mt>140 mj14>400 D4
433  vector<float> pow_k;
434  pow_k.push_back(1); // mt<=140 mj14<=400 R1
435  pow_k.push_back(-1); // mt<=140 mj14>400 R2
436  pow_k.push_back(-1); // mt>140 mj14<=400 D3
437  pow_k.push_back(1); // mt>140 mj14>400 D4
438 
439  float mSigma, pSigma, pred, pred_sys, mSigma_sys, pSigma_sys;
440  size_t nabcd(abcdcuts.size()), digits(2);
441  vector<vector<float> > preds, kappas, fractions;
442  for(size_t ind(0); ind<njbcuts.size(); ind++){
443  bool lowjets(ind%2==0);
444  vector<vector<float> > entries;
445  vector<vector<float> > weights;
446  for(size_t obs(0); obs < pow_pred.size(); obs++) {
447  size_t index(nabcd*ind+obs);
448  entries.push_back(vector<float>());
449  weights.push_back(vector<float>());
450  entries.back().push_back(datayield[index].Yield());
451  weights.back().push_back(1.);
452  } // Loop over observables for data
453  //pred = calcKappa(entries, weights, pow_pred, mSigma, pSigma);
454 
455  vector<vector<float> > kn, kw;
456  float k(1.), kup(1.), kdown(1.);
457  fractions.push_back(vector<float>());
458  for(size_t obs(0); obs < abcdcuts.size(); obs++){
459  size_t index(nabcd*ind+obs);
460  if(!do_other){
461  k *= pow(mcyield[index].Yield(), pow_tot[3+obs]);
462  entries.push_back(vector<float>());
463  weights.push_back(vector<float>());
464  entries.back().push_back(mcyield[index].NEffective());
465  weights.back().push_back(mcyield[index].Weight());
466 
467  kn.push_back(vector<float>());
468  kw.push_back(vector<float>());
469  kn.back().push_back(mcyield[index].NEffective());
470  kw.back().push_back(mcyield[index].Weight());
471  } else {
472  k *= pow(mcyield[index].Yield()+otheryield[index].Yield(), pow_tot[3+obs]);
473  float f(otheryield[index].Yield()/(mcyield[index].Yield()+otheryield[index].Yield()));
474  fractions.back().push_back(f*100);
475 
476  kup *= pow(mcyield[index].Yield()+exp(log(1+syst))*otheryield[index].Yield(), pow_tot[3+obs]);
477  kdown *= pow(mcyield[index].Yield()+exp(-log(1+syst))*otheryield[index].Yield(), pow_tot[3+obs]);
478  }
479  } // Loop over observables for MC
480  k = calcKappa(kn, kw, pow_k, kdown, kup);
481  cout<<"k = "<<k<<" +"<<kup<<" -"<<kdown<<endl;
482  if(!do_other){
483  // Calculating predictions without systematics
484  pred = calcKappa(entries, weights, pow_tot, mSigma, pSigma);
485  if(mSigma<0) mSigma = 0;
486 
487  // Calculating predictions with systematics
488  float totsys = (lowjets?0.51:1.07);
489  pred_sys = calcKappa(entries, weights, pow_tot, mSigma_sys, pSigma_sys, false, false, totsys);
490  if(mSigma_sys < 0) mSigma_sys = 0;
491  preds.push_back(vector<float>({pred, pSigma, mSigma, pred_sys, pSigma_sys, mSigma_sys, k, kup, kdown}));
492  }
493  if(do_other){
494  kup = (kup-k)/k*100;
495  kdown = (kdown-k)/k*100;
496  kappas.push_back(vector<float>({k, kup, kdown}));
497  cout<<"k = "<<RoundNumber(k,2)<<", up "<<setw(5)<<RoundNumber(kup,1)<<"%, down "<<setw(5)
498  <<RoundNumber(kdown,1)<<"%, other fractions ";
499  for(size_t oth(0); oth<fractions[ind].size(); oth++) cout <<setw(4)<<RoundNumber(fractions[ind][oth],1)<<"% ";
500  cout<<" => cuts "<<bincuts[4*ind+3].label_<<endl;
501  }
502 
503  } // Loop over signal bins
504 
506  TString outname = "txt/table_predictions_lumi0p815_"+method+".tex";
507  if(full_lumi) {
508  if(method.Contains("2015")) outname.ReplaceAll("lumi0p815", "lumi2p3");
509  else outname.ReplaceAll("lumi0p815", "lumi2p6");
510  }
511  if(unblind) outname.ReplaceAll("lumi", "unblind_lumi");
512 
513  if(do_other) outname.ReplaceAll("predictions", "other_sys");
514  ofstream out(outname);
515 
516  size_t Ncol = 7;
517  size_t index;
518  out << fixed << setprecision(digits);
519  out << "\\documentclass{article}\n";
520  out << "\\usepackage{amsmath,graphicx,rotating,geometry}\n";
521  out << "\\renewcommand{\\arraystretch}{1.3}\n";
522  out << "\\thispagestyle{empty}\n";
523  out << "\\begin{document}\n\n";
524  out << "\\begin{table}\n";
525  out << "\\centering\n";
526  //out << "\\resizebox{\\textwidth}{!}{\n";
527  if(!do_other){
528  out << "\\begin{tabular}[tbp!]{ l ";
529  for(size_t ind=0; ind<Ncol-1; ind++) out<<"c";
530  out<<"}\\hline\\hline\n";
531  out<<" \\multicolumn{"<<Ncol<<"}{c}{"<<method_s<<"} \\\\ \\hline"<<endl;
532  out<<"${\\cal L}="<<lumi_s<<" & $\\kappa$ & MC & Pred. & Obs. & MC/Obs & $Z_{\\rm bi}$ \\\\ \\hline\\hline\n";
533  if(method.Contains("met150")) out << " \\multicolumn{"<<Ncol<<"}{c}{$150<\\text{MET}\\leq 200$} \\\\ \\hline\n";
534  else if(method.Contains("met500")) out << " \\multicolumn{"<<Ncol<<"}{c}{$\\text{MET}> 500$} \\\\ \\hline\n";
535  else if(method.Contains("met200nb1")) out << " \\multicolumn{"<<Ncol<<"}{c}{$200<\\text{MET}\\leq 350, N_{b}=1$} \\\\ \\hline\n";
536  else if(method.Contains("met350nb1")) out << " \\multicolumn{"<<Ncol<<"}{c}{$350<\\text{MET}\\leq 500, N_{b}=1$} \\\\ \\hline\n";
537  else if(method.Contains("met200nb1")) out << " \\multicolumn{"<<Ncol<<"}{c}{$\\text{MET}> 500, N_{b}=1$} \\\\ \\hline\n";
538  else out << " \\multicolumn{"<<Ncol<<"}{c}{$200<\\text{MET}\\leq 350$} \\\\ \\hline\n";
539 
541  index = 0;
542  if(method.Contains("nb1")) out << "R1: $N_b=1,\\text{all }N_j$";
543  else out << "R1: all $N_b,N_j$";
544  out<<" & & "<<mcyield[index].Yield() <<" & & "
545  << setprecision(0) <<datayield[index].Yield()<<setprecision(digits)
546  <<" & "<<RoundNumber(mcyield[index].Yield(), digits, datayield[index].Yield())<<" & \\\\"<<endl;
548  for(size_t ind(0); ind<ilowmet; ind++){
549  index = nabcd*ind+1;
550  out<<"R2: "<<cuts2tex(njbcuts[ind])<<" & & "<<mcyield[index].Yield() <<" & & "
551  << setprecision(0) <<datayield[index].Yield()<<setprecision(digits)
552  <<" & "<<RoundNumber(mcyield[index].Yield(), digits, datayield[index].Yield())
553  <<" & \\\\"<<endl;
554  }
556  index = 2;
557  if(method.Contains("nb1")) out << region_s<<"3: $N_b=1,\\text{all }N_j$";
558  else out << region_s<<"3: all $N_b,N_j$";
559  out << " & & "<<mcyield[index].Yield() <<" & & "
560  << setprecision(0) << datayield[index].Yield() << setprecision(digits)
561  <<" & "<<RoundNumber(mcyield[index].Yield(), digits, datayield[index].Yield())<<" & \\\\"<<endl;
562  out << "\\hline"<<endl;
564  for(size_t ind(0); ind<ilowmet; ind++){
565  index = nabcd*ind+3;
566  out<<region_s<<"4: "<<cuts2tex(njbcuts_himt[ind])<<" & $"<<preds[ind][6]
567  << "^{+" << preds[ind][7] <<"}_{-" << preds[ind][8]
568  <<"}$ & "<<mcyield[index].Yield() <<" & $"<<preds[ind][0] << "^{+" << preds[ind][1]
569  <<"}_{-" << preds[ind][2] <<"}$ & ";
570  if(unblind) out << setprecision(0) <<datayield[index].Yield()<<setprecision(digits)
571  <<" & "<<RoundNumber(mcyield[index].Yield(), digits, datayield[index].Yield())
572  <<" & "<< Zbi(datayield[index].Yield(), preds[ind][0], preds[ind][1])<<" \\\\"<<endl;
573  else out << blind_s<<" & "<< blind_s<<" \\\\"<<endl;
574  }
575  if(ilowmet<njbcuts.size()){
576  out<<"\\hline\\hline\n ";
577  if(method.Contains("met200nb1"))out<<" \\multicolumn{"<<Ncol<<"}{c}{$200<\\text{MET}\\leq 350, N_{b}\\geq2$} \\\\ \\hline\n";
578  else if(method.Contains("met350nb1")) out << " \\multicolumn{"<<Ncol<<"}{c}{$350<\\text{MET}\\leq 500, N_{b}\\geq2$} \\\\ \\hline\n";
579  else if(method.Contains("met500nb1")) out << " \\multicolumn{"<<Ncol<<"}{c}{$\\text{MET}> 500, N_{b}\\geq2$} \\\\ \\hline\n";
580  else out << " \\multicolumn{"<<Ncol<<"}{c}{$350<\\text{MET}\\leq 500$} \\\\ \\hline\n";
582  index = nabcd*ilowmet;
583  if(method.Contains("nb1")) out << "R1: $N_b=1,\\text{all }N_j$";
584  else out << "R1: all $N_b,N_j$";
585  out << " & & "<<mcyield[index].Yield() <<" & & "
586  << setprecision(0) <<datayield[index].Yield()<<setprecision(digits)
587  <<" & "<<RoundNumber(mcyield[index].Yield(), digits, datayield[index].Yield())<<" & \\\\"<<endl;
589  for(size_t ind(ilowmet); ind<njbcuts.size(); ind++){
590  index = nabcd*ind+1;
591  out<<"R2: "<<cuts2tex(njbcuts[ind])<<" & & "<<mcyield[index].Yield() <<" & & "
592  << setprecision(0) <<datayield[index].Yield()<<setprecision(digits)
593  <<" & "<<RoundNumber(mcyield[index].Yield(), digits, datayield[index].Yield())<<" & \\\\"<<endl;
594  }
595 
597  index = nabcd*ilowmet+2;
598  if(method.Contains("nb1")) out << region_s<<"3: $N_b=1,\\text{all }N_j$";
599  else out << region_s<<"3: all $N_b,N_j$";
600  out << " & & "<<mcyield[index].Yield() <<" & & "
601  << setprecision(0) <<datayield[index].Yield()<<setprecision(digits)
602  <<" & "<<RoundNumber(mcyield[index].Yield(), digits, datayield[index].Yield())<<" & \\\\"<<endl;
603  out << "\\hline"<<endl;
605  for(size_t ind(ilowmet); ind<njbcuts.size(); ind++){
606  index = nabcd*ind+3;
607  out<<region_s<<"4: "<<cuts2tex(njbcuts_himt[ind])<<" & $"<<preds[ind][6]
608  << "^{+" << preds[ind][7] <<"}_{-" << preds[ind][8]
609  <<"}$ & "<<mcyield[index].Yield() <<" & $"<<preds[ind][0] << "^{+" << preds[ind][1]
610  <<"}_{-" << preds[ind][2] <<"}$ & ";
611  if(unblind) out<< setprecision(0) <<datayield[index].Yield() <<setprecision(digits)
612  <<" & "<<RoundNumber(mcyield[index].Yield(), digits, datayield[index].Yield())
613  <<" & "<< Zbi(datayield[index].Yield(), preds[ind][0], preds[ind][1])<<" \\\\"<<endl;
614  else out<< blind_s<<" & "<< blind_s<<" \\\\"<<endl;
615  }
616  }
617  } else {
618  out << "\n\\begin{tabular}[tbp!]{ l rrr |";
619  for(size_t obs(0); obs < abcdcuts.size(); obs++) out<<"r";
620  out << "}\\hline\\hline\n";
621  out << " & & non-$t\\bar{t}\\times"<<RoundNumber(1+syst,1)<<"$ & non-$t\\bar{t}\\times"
622  <<RoundNumber(exp(-log(1+syst)),1)<<"$ & \\multicolumn{4}{c}{Fraction of non-$t\\bar{t}$ bkg.} \\\\ \n";
623  out << "Bin & $\\kappa$ & $\\Delta\\kappa$ [\\%] & $\\Delta\\kappa$ [\\%]";
624  for(size_t obs(0); obs < abcdcuts.size(); obs++) out<<" & $f_{\\text{R"<<obs+1<<"}}$ [\\%]";
625  out << " \\\\ \\hline\\hline\n";
626  for(size_t ind(0); ind<njbcuts.size(); ind++){
627  out<<cuts2tex(njbcuts[ind])<<" & "<<RoundNumber(kappas[ind][0],2) <<" & "<<kappas[ind][1] <<" & "<<kappas[ind][2];
628  for(size_t obs(0); obs < abcdcuts.size(); obs++) out<<" & "<<fractions[ind][obs];
629  out<<" \\\\"<<endl;
630  if(ind==ilowmet-1) out<<"\\hline"<<endl;
631  }
632  } // if(!do_other)
633 
634  out<< "\\hline\\hline\n\\end{tabular}"<<endl;
635  //out << "}\n";
636  out << "\\end{table}\n\n";
637  out << "\\end{document}\n";
638  out.close();
639  TString pdfname(outname);
640  cout<<endl<<" pdflatex "<<outname<<endl;
641 
642  time(&endtime);
643  cout<<endl<<"Calculation took "<<difftime(endtime, begtime)<<" seconds"<<endl<<endl;
644 
645 }
646 
647 
648 void GetOptions(int argc, char *argv[]){
649  while(true){
650  static struct option long_options[] = {
651  {"method", required_argument, 0, 'm'},
652  {"unblind", no_argument, 0, 'u'},
653  {"full_lumi", no_argument, 0, 'f'},
654  {0, 0, 0, 0}
655  };
656 
657  char opt = -1;
658  int option_index;
659  opt = getopt_long(argc, argv, "m:uf", long_options, &option_index);
660  if(opt == -1) break;
661 
662  string optname;
663  switch(opt){
664  case 'm':
665  method = optarg;
666  break;
667  case 'u':
668  really_unblind = true;
669  break;
670  case 'f':
671  full_lumi = true;
672  break;
673  case 0:
674  break;
675  default:
676  printf("Bad option! getopt_long returned character code 0%o\n", opt);
677  break;
678  }
679  }
680 }
Proc
shared_ptr< Process > Proc(const string process_name, Process::Type type, int color, const set< string > &files, const string &cut="1")
Definition: table_2016-06-03_freeze_preapp_an.cxx:35
anonymous_namespace{table_2016-06-03_freeze_preapp_an.cxx}::method
TString method
Definition: table_2016-06-03_freeze_preapp_an.cxx:24
main
int main(int argc, char *argv[])
Definition: table_2016-06-03_freeze_preapp_an.cxx:54
Table::Yield
std::vector< GammaParams > Yield(const Process *process, double luminosity) const
Definition: table.cpp:158
anonymous_namespace{table_2016-06-03_freeze_preapp_an.cxx}::really_unblind
bool really_unblind
Definition: table_2016-06-03_freeze_preapp_an.cxx:30
anonymous_namespace{table_2016-06-03_freeze_preapp_an.cxx}::blind_s
TString blind_s
Definition: table_2016-06-03_freeze_preapp_an.cxx:29
anonymous_namespace{table_2016-06-03_freeze_preapp_an.cxx}::single_thread
bool single_thread
Definition: table_2016-06-03_freeze_preapp_an.cxx:31
GetOptions
void GetOptions(int argc, char *argv[])
Definition: table_2016-06-03_freeze_preapp_an.cxx:648
plot_maker.hpp
PlotMaker::Figures
const std::vector< std::unique_ptr< Figure > > & Figures() const
Definition: plot_maker.cpp:63
std
STL namespace.
Contains
bool Contains(const std::string &str, const std::string &pat)
Definition: utilities.cpp:44
Zbi
TString Zbi(double Nobs, double Nbkg, double Ebkg)
Definition: table_2016-06-03_freeze_preapp_an.cxx:43
PlotMaker::multithreaded_
bool multithreaded_
Definition: plot_maker.hpp:43
execute
std::string execute(const std::string &cmd)
Definition: utilities.cpp:65
named_func.hpp
calcKappa
double calcKappa(std::vector< std::vector< float > > &entries, std::vector< std::vector< float > > &weights, std::vector< float > &powers, float &mSigma, float &pSigma, bool do_data=false, bool verbose=false, double syst=-1., bool do_plot=false, int nrep=100000)
Definition: utilities.cpp:469
PlotMaker::Push
FigureType & Push(Args &&...args)
Definition: plot_maker.hpp:24
anonymous_namespace{table_2016-06-03_freeze_preapp_an.cxx}::lumi
double lumi(0.815)
Table::BackgroundYield
std::vector< GammaParams > BackgroundYield(double luminosity) const
Definition: table.cpp:175
anonymous_namespace{table_2016-06-03_freeze_preapp_an.cxx}::do_other
bool do_other(false)
PlotMaker
Organizes efficient production of plots with single loop over each process.
Definition: plot_maker.hpp:14
RoundNumber
TString RoundNumber(double num, int decimals, double denom=1.)
Definition: utilities.cpp:361
anonymous_namespace{plot_ratios.cxx}::other
Definition: plot_ratios.cxx:41
Table::DataYield
std::vector< GammaParams > DataYield() const
Definition: table.cpp:188
Table
Definition: table.hpp:15
anonymous_namespace{table_2016-06-03_freeze_preapp_an.cxx}::syst
float syst
Definition: table_2016-06-03_freeze_preapp_an.cxx:28
anonymous_namespace{table_preds.cxx}::unblind
bool unblind
Definition: table_preds.cxx:45
PlotMaker::MakePlots
void MakePlots(double luminosity, const std::string &subdir="")
Prints all added plots with given luminosity.
Definition: plot_maker.cpp:54
anonymous_namespace{table_2016-06-03_freeze_preapp_an.cxx}::full_lumi
bool full_lumi(false)
Palette
Loads colors from a text configuration file.
Definition: palette.hpp:8
Generated by   doxygen 1.8.11