ra4_macros  bede988c286599a3a84b77a4d788ac0a971e89f9
plot_closure.cxx
Go to the documentation of this file.
1 #include <iostream>
2 #include <vector>
3 
4 #include "TString.h"
5 #include "TGraphErrors.h"
6 
7 #include "TChain.h"
8 #include "TH1D.h"
9 #include "TCanvas.h"
10 #include "TLegend.h"
11 #include "TLine.h"
12 #include "TColor.h"
13 #include "TMath.h"
14 #include "TF1.h"
15 #include "TLatex.h"
16 
17 #include "styles.hpp"
18 #include "utilities.hpp"
19 #include "utilities_macros.hpp"
20 
21 using namespace std;
22 
23 double GetWeight(sfeats Sample, TString luminosity);
24 void GetYield(double& yield, double& yield_err, TString baseline,TString external,TString mj_cut,TString mt_cut, sfeats Sample,TString luminosity);
25 void GetRatio(double& ratio, double& error, double num, double den, double num_err, double den_err, double weight);
26 void GetKappa(double& kappa, double& error, double n1, double n2, double n3, double n4, double n1_err, double n2_err, double n3_err, double n4_err, double weight);
27 void MakeGraphs(double rmj1[],double rmj2[],double rmj1_err[],double rmj2_err[],double rmt1[],double rmt2[],double rmt1_err[],double rmt2_err[],double kappa[],double kappa_err[], const vector<TString> binlabels,TString sampleName, TString baseline, TString mj_cut, TString mt_cut);
28 
29 int main(){
30  TString luminosity = "10";
31  // TString folder="/cms5r0/ald77/archive/2015_05_25/skim/";
32  // TString folder_noskim="/cms5r0/ald77/archive/2015_05_25/";
33  TString folder="/cms5r0/rohan/2015_05_25/skim_met100/";
34 
35  vector<TString> s_t1t;
36  s_t1t.push_back(folder+"*T1tttt*1500_*PU20*");
37  vector<TString> s_t1tc;
38  s_t1tc.push_back(folder+"*T1tttt*1200_*PU20*");
39  vector<TString> s_tt;
40  s_tt.push_back(folder+"*_TTJet*");
41  // vector<TString> s_tt_noskim;
42  // s_tt_noskim.push_back(folder_noskim+"*_TTJet*12.root");
43  vector<TString> s_wjets;
44  s_wjets.push_back(folder+"*WJetsToLNu_HT*");
45  vector<TString> s_singlet;
46  s_singlet.push_back(folder+"*_T*channel*");
47  vector<TString> s_ttv;
48  s_ttv.push_back(folder+"*TTW*");
49  s_ttv.push_back(folder+"*TTZ*");
50  vector<TString> s_other;
51  s_other.push_back(folder+"*QCD_HT*");
52  s_other.push_back(folder+"*_ZJet*");
53  s_other.push_back(folder+"*DY*");
54  s_other.push_back(folder+"*WH_HToBB*");
55 
56  //TString baseline = "ht>500&&met>100&&nbm==1&&(nmus+nels)==1&&njets>=7"; // MET 100, Nb=1
57  // TString baseline = "ht>500&&met>100&&nbm==1&&(nmus+nels)==1&&njets>=7&&ntks_chg==0"; // ITV, MET 100, Nb=1
58 
59  TString baseline = "ht>500&&met>100&&nbm>=2&&(nmus+nels)==1&&njets>=7"; // MET 100, Nb>=2
60  //TString baseline = "ht>500&&met>100&&nbm>=2&&(nmus+nels)==1&&njets>=7&&ntks_chg==0"; // ITV, MET 100, Nb>=2
61 
62  //TString baseline = "ht>500&&met>200&&nbm>=2&&(nmus+nels)==1&&njets>=7"; // MET 200, Nb>=2
63  //TString baseline = "ht>500&&met>200&&nbm>=2&&(nmus+nels)==1&&njets>=7&&ntks_chg==0"; // ITV, MET 200, Nb>=2
64 
65 
66  vector<sfeats> samples_MJ;
67  //samples_MJ.push_back(sfeats(s_tt, "t#bar{t}", 31, 1));
68  //samples_MJ.push_back(sfeats(s_tt, "2#it{l} t#bar{t}", 31,1,"ntruleps>=2"));
69  samples_MJ.push_back(sfeats(s_tt, "1#it{l} t#bar{t}", 31,1,"ntruleps<=1"));
70 
71  // METHOD 1
72  TString mj_cut = "mj>600";
73  TString mt_cut = "mt>140";
74  TString not_mt_cut = invertcut(mt_cut);
75  TString not_mj_cut = invertcut(mj_cut);
76 
77  const int nExt = 6; // Must match size of external vector
78  vector<TString> external_meth1;
79  external_meth1.push_back("njets>=7&&njets<=8&&met>100&&met<=200");
80  external_meth1.push_back("njets>=7&&njets<=8&&met>200&&met<=400");
81  external_meth1.push_back("njets>=7&&njets<=8&&met>400");
82  external_meth1.push_back("njets>=9&&met>100&&met<=200");
83  external_meth1.push_back("njets>=9&&met>200&&met<=400");
84  external_meth1.push_back("njets>=9&&met>400");
85 
86  double weight = GetWeight(samples_MJ.at(0), luminosity);
87 
88  double N1_meth1[nExt],N2_meth1[nExt],N3_meth1[nExt],N4_meth1[nExt];
89  double N1_meth1_err[nExt],N2_meth1_err[nExt],N3_meth1_err[nExt],N4_meth1_err[nExt];
90  for(int iR=0; iR<nExt; iR++){
91  GetYield(N1_meth1[iR], N1_meth1_err[iR], baseline, external_meth1.at(iR), not_mj_cut, not_mt_cut, samples_MJ.at(0), luminosity);
92  GetYield(N2_meth1[iR], N2_meth1_err[iR], baseline, external_meth1.at(iR), mj_cut, not_mt_cut, samples_MJ.at(0), luminosity);
93  GetYield(N3_meth1[iR], N3_meth1_err[iR], baseline, external_meth1.at(iR), not_mj_cut, mt_cut, samples_MJ.at(0), luminosity);
94  GetYield(N4_meth1[iR], N4_meth1_err[iR], baseline, external_meth1.at(iR), mj_cut, mt_cut, samples_MJ.at(0), luminosity);
95  }
96 
97  double rmj43[nExt], rmj21[nExt], rmj43_err[nExt], rmj21_err[nExt];
98  double rmt42[nExt], rmt31[nExt], rmt42_err[nExt], rmt31_err[nExt];
99  for(int iR=0; iR<nExt; iR++){
100  GetRatio(rmj21[iR], rmj21_err[iR], N2_meth1[iR], N1_meth1[iR], N2_meth1_err[iR], N1_meth1_err[iR], weight);
101  GetRatio(rmj43[iR], rmj43_err[iR], N4_meth1[iR], N3_meth1[iR], N4_meth1_err[iR], N3_meth1_err[iR], weight);
102  GetRatio(rmt31[iR], rmt31_err[iR], N3_meth1[iR], N1_meth1[iR], N3_meth1_err[iR], N1_meth1_err[iR], weight);
103  GetRatio(rmt42[iR], rmt42_err[iR], N4_meth1[iR], N2_meth1[iR], N4_meth1_err[iR], N2_meth1_err[iR], weight);
104  }
105 
106  double kappa[nExt], kappa_err[nExt];
107  for(int iR=0; iR<nExt; iR++){
108  GetKappa(kappa[iR], kappa_err[iR], N1_meth1[iR], N2_meth1[iR], N3_meth1[iR], N4_meth1[iR], N1_meth1_err[iR], N2_meth1_err[iR], N3_meth1_err[iR], N4_meth1_err[iR], weight);
109  }
110 
111  vector<TString> binlabels; //push back twice for each njets bin
112  binlabels.push_back("low MET");binlabels.push_back("med MET");binlabels.push_back("high MET");
113  binlabels.push_back("low MET");binlabels.push_back("med MET");binlabels.push_back("high MET");
114  TString binName = "method1";
115  if (baseline.Contains("met>100")) binName.Append("_met100");
116  if (baseline.Contains("nbm==1")) binName.Append("_nb1");
117  binlabels.push_back(binName); // Last bin is for method name.
118 
119  MakeGraphs(rmj43, rmj21, rmj43_err, rmj21_err, rmt42, rmt31, rmt42_err, rmt31_err, kappa, kappa_err, binlabels, samples_MJ.at(0).label, baseline, mj_cut, mt_cut);
120 
121  // METHOD 2
122  mj_cut = "mj>400";
123  mt_cut = "mt>140";
124  not_mt_cut = invertcut(mt_cut);
125  not_mj_cut = invertcut(mj_cut);
126 
127  const int nExt_2 = 6; // Must match size of external vector
128  vector<TString> externalA_meth2;
129  externalA_meth2.push_back("met>100&&met<=200");
130  externalA_meth2.push_back("met>200&&met<=400");
131  externalA_meth2.push_back("met>400");
132  externalA_meth2.push_back("met>100&&met<=200");
133  externalA_meth2.push_back("met>200&&met<=400");
134  externalA_meth2.push_back("met>400");
135  vector<TString> externalB_meth2;
136  externalB_meth2.push_back("njets>=7&&njets<=8&&met>100&&met<=200");
137  externalB_meth2.push_back("njets>=7&&njets<=8&&met>200&&met<=400");
138  externalB_meth2.push_back("njets>=7&&njets<=8&&met>400");
139  externalB_meth2.push_back("njets>=9&&met>100&&met<=200");
140  externalB_meth2.push_back("njets>=9&&met>200&&met<=400");
141  externalB_meth2.push_back("njets>=9&&met>400");
142 
143  double N1_meth2[nExt_2],N2_meth2[nExt_2],N3_meth2[nExt_2],N4_meth2[nExt_2];
144  double N1_meth2_err[nExt_2],N2_meth2_err[nExt_2],N3_meth2_err[nExt_2],N4_meth2_err[nExt_2];
145  for(int iR=0; iR<nExt_2; iR++){
146  GetYield(N1_meth2[iR], N1_meth2_err[iR], baseline, externalA_meth2.at(iR), not_mj_cut, not_mt_cut, samples_MJ.at(0), luminosity);
147  GetYield(N2_meth2[iR], N2_meth2_err[iR], baseline, externalB_meth2.at(iR), mj_cut, not_mt_cut, samples_MJ.at(0), luminosity);
148  GetYield(N3_meth2[iR], N3_meth2_err[iR], baseline, externalA_meth2.at(iR), not_mj_cut, mt_cut, samples_MJ.at(0), luminosity);
149  GetYield(N4_meth2[iR], N4_meth2_err[iR], baseline, externalB_meth2.at(iR), mj_cut, mt_cut, samples_MJ.at(0), luminosity);
150  }
151 
152  double rmj43_2[nExt_2], rmj21_2[nExt_2], rmj43_err_2[nExt_2], rmj21_err_2[nExt_2];
153  double rmt42_2[nExt_2], rmt31_2[nExt_2], rmt42_err_2[nExt_2], rmt31_err_2[nExt_2];
154  for(int iR=0; iR<nExt_2; iR++){
155  GetRatio(rmj21_2[iR], rmj21_err_2[iR], N2_meth2[iR], N1_meth2[iR], N2_meth2_err[iR], N1_meth2_err[iR], weight);
156  GetRatio(rmj43_2[iR], rmj43_err_2[iR], N4_meth2[iR], N3_meth2[iR], N4_meth2_err[iR], N3_meth2_err[iR], weight);
157  GetRatio(rmt31_2[iR], rmt31_err_2[iR], N3_meth2[iR], N1_meth2[iR], N3_meth2_err[iR], N1_meth2_err[iR], weight);
158  GetRatio(rmt42_2[iR], rmt42_err_2[iR], N4_meth2[iR], N2_meth2[iR], N4_meth2_err[iR], N2_meth2_err[iR], weight);
159  }
160 
161  double kappa_2[nExt_2], kappa_err_2[nExt_2];
162  for(int iR=0; iR<nExt_2; iR++){
163  GetKappa(kappa_2[iR], kappa_err_2[iR], N1_meth2[iR], N2_meth2[iR], N3_meth2[iR], N4_meth2[iR], N1_meth2_err[iR], N2_meth2_err[iR], N3_meth2_err[iR], N4_meth2_err[iR], weight);
164  }
165 
166  vector<TString> binlabels_2; //push back twice for each njets bin
167  binlabels_2.push_back("low MET");binlabels_2.push_back("med MET");binlabels_2.push_back("high MET");
168  binlabels_2.push_back("low MET");binlabels_2.push_back("med MET");binlabels_2.push_back("high MET");
169 
170  TString binName_2 = "method2";
171  if (baseline.Contains("met>100")) binName_2.Append("_met100");
172  if (baseline.Contains("nbm==1")) binName_2.Append("_nb1");
173  binlabels_2.push_back(binName_2); // Last bin is for method name.
174 
175  MakeGraphs(rmj43_2, rmj21_2, rmj43_err_2, rmj21_err_2, rmt42_2, rmt31_2, rmt42_err_2, rmt31_err_2, kappa_2, kappa_err_2, binlabels_2, samples_MJ.at(0).label, baseline, mj_cut, mt_cut);
176 
177  /* // METHOD 3
178  mj_cut = "mj>400";
179  mt_cut = "mt>140";
180  not_mt_cut = invertcut(mt_cut);
181  not_mj_cut = invertcut(mj_cut);
182 
183  const int nExt_3 = 6; // Must match size of largest external vector
184  vector<TString> externalA_meth3;
185  externalA_meth3.push_back("met>200&&met<=400&&nbm==2");
186  externalA_meth3.push_back("met>200&&met<=400&&nbm>=3");
187  externalA_meth3.push_back("met>400");
188  externalA_meth3.push_back("met>200&&met<=400&&nbm==2");
189  externalA_meth3.push_back("met>200&&met<=400&&nbm>=3");
190  externalA_meth3.push_back("met>400");
191  vector<TString> externalB_meth3;
192  externalB_meth3.push_back("njets>=7&&njets<=8&&met>200&&met<=400&&nbm==2");
193  externalB_meth3.push_back("njets>=7&&njets<=8&&met>200&&met<=400&&nbm>=3");
194  externalB_meth3.push_back("njets>=7&&njets<=8&&met>400");
195  externalB_meth3.push_back("njets>=9&&met>200&&met<=400&&nbm==2");
196  externalB_meth3.push_back("njets>=9&&met>200&&met<=400&&nbm>=3");
197  externalB_meth3.push_back("njets>=9&&met>400");
198 
199  double N1_meth3[nExt_3],N2_meth3[nExt_3],N3_meth3[nExt_3],N4_meth3[nExt_3];
200  double N1_meth3_err[nExt_3],N2_meth3_err[nExt_3],N3_meth3_err[nExt_3],N4_meth3_err[nExt_3];
201  for(int iR=0; iR<nExt_3; iR++){
202  GetYield(N1_meth3[iR], N1_meth3_err[iR], baseline, externalA_meth3.at(iR), not_mj_cut, not_mt_cut, samples_MJ.at(0), luminosity);
203  GetYield(N2_meth3[iR], N2_meth3_err[iR], baseline, externalB_meth3.at(iR), mj_cut, not_mt_cut, samples_MJ.at(0), luminosity);
204  GetYield(N3_meth3[iR], N3_meth3_err[iR], baseline, externalA_meth3.at(iR), not_mj_cut, mt_cut, samples_MJ.at(0), luminosity);
205  GetYield(N4_meth3[iR], N4_meth3_err[iR], baseline, externalB_meth3.at(iR), mj_cut, mt_cut, samples_MJ.at(0), luminosity);
206  }
207 
208  double rmj43_3[nExt_3], rmj21_3[nExt_3], rmj43_err_3[nExt_3], rmj21_err_3[nExt_3];
209  double rmt42_3[nExt_3], rmt31_3[nExt_3], rmt42_err_3[nExt_3], rmt31_err_3[nExt_3];
210  for(int iR=0; iR<nExt_3; iR++){
211  GetRatio(rmj21_3[iR], rmj21_err_3[iR], N2_meth3[iR], N1_meth3[iR], N2_meth3_err[iR], N1_meth3_err[iR], weight);
212  GetRatio(rmj43_3[iR], rmj43_err_3[iR], N4_meth3[iR], N3_meth3[iR], N4_meth3_err[iR], N3_meth3_err[iR], weight);
213  GetRatio(rmt31_3[iR], rmt31_err_3[iR], N3_meth3[iR], N1_meth3[iR], N3_meth3_err[iR], N1_meth3_err[iR], weight);
214  GetRatio(rmt42_3[iR], rmt42_err_3[iR], N4_meth3[iR], N2_meth3[iR], N4_meth3_err[iR], N2_meth3_err[iR], weight);
215  }
216 
217  double kappa_3[nExt_3], kappa_err_3[nExt_3];
218  for(int iR=0; iR<nExt_3; iR++){
219  GetKappa(kappa_3[iR], kappa_err_3[iR], N1_meth3[iR], N2_meth3[iR], N3_meth3[iR], N4_meth3[iR], N1_meth3_err[iR], N2_meth3_err[iR], N3_meth3_err[iR], N4_meth3_err[iR], weight);
220  }
221 
222  vector<TString> binlabels_3; //push back twice for each njets bin
223  binlabels_3.push_back("low MET, nb=2");binlabels_3.push_back("low MET, nb#geq3");binlabels_3.push_back("high MET");
224  binlabels_3.push_back("low MET, nb=2");binlabels_3.push_back("low MET, nb#geq3");binlabels_3.push_back("high MET");
225 
226  TString binName_3 = "method3";
227  if (baseline.Contains("met>100")) binName_3.Append("_met100");
228  if (baseline.Contains("nbm==1")) binName_3.Append("_nb1");
229  binlabels_3.push_back(binName_3); // Last bin is for method name.
230 
231  MakeGraphs(rmj43_3, rmj21_3, rmj43_err_3, rmj21_err_3, rmt42_3, rmt31_3, rmt42_err_3, rmt31_err_3, kappa_3, kappa_err_3, binlabels_3, samples_MJ.at(0).label, baseline, mj_cut, mt_cut);
232  */ //METHOD 3
233 
234  // METHOD 0
235  mj_cut = "mj>400";
236  mt_cut = "mt>140";
237  not_mt_cut = invertcut(mt_cut);
238  not_mj_cut = invertcut(mj_cut);
239 
240  const int nExt_0 = 6; // Must match size of external vector
241  vector<TString> externalA_meth0;
242  externalA_meth0.push_back("met>100");
243  externalA_meth0.push_back("met>100");
244  externalA_meth0.push_back("met>100");
245  externalA_meth0.push_back("met>100");
246  externalA_meth0.push_back("met>100");
247  externalA_meth0.push_back("met>100");
248  vector<TString> externalB_meth0;
249  externalB_meth0.push_back("njets>=7&&njets<=8&&met>100&&met<=200");
250  externalB_meth0.push_back("njets>=7&&njets<=8&&met>200&&met<=400");
251  externalB_meth0.push_back("njets>=7&&njets<=8&&met>400");
252  externalB_meth0.push_back("njets>=9&&met>100&&met<=200");
253  externalB_meth0.push_back("njets>=9&&met>200&&met<=400");
254  externalB_meth0.push_back("njets>=9&&met>400");
255 
256  double N1_meth0[nExt_0],N2_meth0[nExt_0],N3_meth0[nExt_0],N4_meth0[nExt_0];
257  double N1_meth0_err[nExt_0],N2_meth0_err[nExt_0],N3_meth0_err[nExt_0],N4_meth0_err[nExt_0];
258  for(int iR=0; iR<nExt_0; iR++){
259  GetYield(N1_meth0[iR], N1_meth0_err[iR], baseline, externalA_meth0.at(iR), not_mj_cut, not_mt_cut, samples_MJ.at(0), luminosity);
260  GetYield(N2_meth0[iR], N2_meth0_err[iR], baseline, externalB_meth0.at(iR), mj_cut, not_mt_cut, samples_MJ.at(0), luminosity);
261  GetYield(N3_meth0[iR], N3_meth0_err[iR], baseline, externalA_meth0.at(iR), not_mj_cut, mt_cut, samples_MJ.at(0), luminosity);
262  GetYield(N4_meth0[iR], N4_meth0_err[iR], baseline, externalB_meth0.at(iR), mj_cut, mt_cut, samples_MJ.at(0), luminosity);
263  }
264 
265  double rmj43_0[nExt_0], rmj21_0[nExt_0], rmj43_err_0[nExt_0], rmj21_err_0[nExt_0];
266  double rmt42_0[nExt_0], rmt31_0[nExt_0], rmt42_err_0[nExt_0], rmt31_err_0[nExt_0];
267  for(int iR=0; iR<nExt_0; iR++){
268  GetRatio(rmj21_0[iR], rmj21_err_0[iR], N2_meth0[iR], N1_meth0[iR], N2_meth0_err[iR], N1_meth0_err[iR], weight);
269  GetRatio(rmj43_0[iR], rmj43_err_0[iR], N4_meth0[iR], N3_meth0[iR], N4_meth0_err[iR], N3_meth0_err[iR], weight);
270  GetRatio(rmt31_0[iR], rmt31_err_0[iR], N3_meth0[iR], N1_meth0[iR], N3_meth0_err[iR], N1_meth0_err[iR], weight);
271  GetRatio(rmt42_0[iR], rmt42_err_0[iR], N4_meth0[iR], N2_meth0[iR], N4_meth0_err[iR], N2_meth0_err[iR], weight);
272  }
273 
274  double kappa_0[nExt_0], kappa_err_0[nExt_0];
275  for(int iR=0; iR<nExt_0; iR++){
276  GetKappa(kappa_0[iR], kappa_err_0[iR], N1_meth0[iR], N2_meth0[iR], N3_meth0[iR], N4_meth0[iR], N1_meth0_err[iR], N2_meth0_err[iR], N3_meth0_err[iR], N4_meth0_err[iR], weight);
277  }
278 
279  vector<TString> binlabels_0; //push back twice for each njets bin
280  binlabels_0.push_back("low MET");binlabels_0.push_back("med MET");binlabels_0.push_back("high MET");
281  binlabels_0.push_back("low MET");binlabels_0.push_back("med MET");binlabels_0.push_back("high MET");
282 
283  TString binName_0 = "method0";
284  if (baseline.Contains("met>100")) binName_0.Append("_met100");
285  if (baseline.Contains("nbm==1")) binName_0.Append("_nb1");
286  binlabels_0.push_back(binName_0); // Last bin is for method name.
287 
288  MakeGraphs(rmj43_0, rmj21_0, rmj43_err_0, rmj21_err_0, rmt42_0, rmt31_0, rmt42_err_0, rmt31_err_0, kappa_0, kappa_err_0, binlabels_0, samples_MJ.at(0).label, baseline, mj_cut, mt_cut);
289 }
290 
291 void MakeGraphs(double rmj1[],double rmj2[],double rmj1_err[],double rmj2_err[],double rmt1[],double rmt2[],double rmt1_err[],double rmt2_err[],double kappa[],double kappa_err[], const vector<TString> binlabels,TString sampleName, TString baseline, TString mj_cut, TString mt_cut){
292 
293  TString sampleNameText = sampleName;
294  if(sampleName.Contains("2#it{l} t#bar{t}")) { sampleNameText = sampleName; sampleNameText.ReplaceAll("2#it{l} t#bar{t}","ttbar_2l"); }
295  else if(sampleName.Contains("1#it{l} t#bar{t}")) { sampleNameText = sampleName; sampleNameText.ReplaceAll("1#it{l} t#bar{t}","ttbar_1l"); }
296  else if(sampleName.Contains("t#bar{t}")) { sampleNameText = sampleName; sampleNameText.ReplaceAll("t#bar{t}","ttbar"); }
297 
298 
299  vector<TGraphErrors*> graphsMJ;
300  vector<TGraphErrors*> graphsMT;
301  TGraphErrors* kappas;
302  const int nExt = binlabels.size()-1;
303 
304  vector<double> x, x1, x2, x_err;
305  for(int i=0; i<nExt; i++){
306  x.push_back(i);
307  x1.push_back(i-0.075);
308  x2.push_back(i+0.075);
309  x_err.push_back(0);
310  }
311 
312  TGraphErrors *mj1 = new TGraphErrors(nExt,&x1[0],rmj1,&x_err[0],rmj1_err); // &foo[0] convert vector "foo" to an array because TGraphErrors only accepts arrays
313  graphsMJ.push_back(mj1);
314 
315  TGraphErrors *mj2 = new TGraphErrors(nExt,&x2[0],rmj2,&x_err[0],rmj2_err); // &foo[0] converts vector "foo" to an array because TGraphErrors only accepts arrays
316  graphsMJ.push_back(mj2);
317 
318  float maxMJ = 0.01;
319  for(int imax = 0;imax<nExt;imax++){
320  if(rmj1[imax]>maxMJ) maxMJ=rmj1[imax];
321  if(rmj2[imax]>maxMJ) maxMJ=rmj2[imax];
322  }
323 
324  TGraphErrors *mt1 = new TGraphErrors(nExt,&x1[0],rmt1,&x_err[0],rmt1_err); // &foo[0] convert vector "foo" to an array because TGraphErrors only accepts arrays
325  graphsMT.push_back(mt1);
326 
327  TGraphErrors *mt2 = new TGraphErrors(nExt,&x2[0],rmt2,&x_err[0],rmt2_err); // &foo[0] converts vector "foo" to an array because TGraphErrors only accepts arrays
328  graphsMT.push_back(mt2);
329 
330  float maxMT = 0.01;
331  for(int imax = 0;imax<nExt;imax++){
332  if(rmt1[imax]>maxMT) maxMT=rmt1[imax];
333  if(rmt2[imax]>maxMT) maxMT=rmt2[imax];
334  }
335 
336  kappas = new TGraphErrors(nExt,&x[0],kappa,&x_err[0],kappa_err);
337 
338  styles style("RA4"); style.setDefaultStyle();
339  TCanvas can;
340 
341  //make hist to define pad with labels on x-axis
342  TH1F *h = new TH1F("for_axis_label",cuts2title(baseline),nExt,-0.5,nExt-0.5);
343  for (int ih=1;ih<=nExt;ih++){
344  h->GetXaxis()->SetBinLabel(ih,binlabels[ih-1]);
345  }
346 
347  h->SetMaximum(1.5*maxMJ);
348  h->SetMinimum(0.0);
349 
350  h->GetYaxis()->SetTitle("R_{MJ}");
351  h->GetXaxis()->SetLabelSize(0.045);
352  h->Draw();
353 
354  double legX = 0.65, legY = 0.89, legSingle = 0.14;
355  double legW = 0.22, legH = legSingle;
356  TLegend leg(legX, legY-legH, legX+legW, legY);
357  leg.SetTextSize(0.057); leg.SetFillColor(0); leg.SetFillStyle(0); leg.SetBorderSize(0);
358  leg.SetTextFont(132);
359 
360  graphsMJ.at(0)->SetMarkerStyle(20);
361  graphsMJ.at(0)->SetMarkerSize(1.2);
362  graphsMJ.at(0)->SetMarkerColor(31);
363  graphsMJ.at(0)->SetLineColor(31);
364  graphsMJ.at(0)->Draw("P");
365  leg.AddEntry(graphsMJ.at(0), sampleName+" "+cuts2title(mt_cut),"p");
366 
367  graphsMJ.at(1)->SetMarkerStyle(22);
368  graphsMJ.at(1)->SetMarkerSize(1.2);
369  graphsMJ.at(1)->SetMarkerColor(46);
370  graphsMJ.at(1)->SetLineColor(46);
371  graphsMJ.at(1)->Draw("P");
372  leg.AddEntry(graphsMJ.at(1), sampleName+" "+cuts2title(invertcut(mt_cut)),"p");
373 
374  leg.Draw("p");
375 
376  TLine line; line.SetLineColor(28); line.SetLineWidth(4); line.SetLineStyle(2);
377  line.DrawLine((nExt-1)/2., 0, (nExt-1)/2., 1.5*maxMJ);
378  TLatex *text78 = new TLatex(0.35,0.03,"n_{jets}= 7-8");
379  text78->SetNDC();
380  text78->SetTextSize(0.04);
381  text78->SetLineWidth(2);
382  text78->Draw();
383  TLatex *text9 = new TLatex(0.7,0.03,"n_{jets}#geq 9");
384  text9->SetNDC();
385  text9->SetTextSize(0.04);
386  text9->SetLineWidth(2);
387  text9->Draw();
388 
389  TString pname1 = "plots/closure/"+binlabels.at(nExt)+"_"+sampleNameText+"_rmj.pdf";
390  TString pname1root = "plots/closure/"+binlabels.at(nExt)+"_"+sampleNameText+"_rmj.root";
391  if(baseline.Contains("ntks_chg==0")) { pname1.ReplaceAll("_rmj","_rmj_ITV"); pname1root.ReplaceAll("_rmj","_rmj_ITV"); }
392  can.SaveAs(pname1);
393  can.SaveAs(pname1root);
394 
395 
396  graphsMJ.at(0)->Delete();
397  graphsMJ.at(1)->Delete();
398 
399  TCanvas can2;
400 
401  h->SetMaximum(1.5*maxMT);
402  h->SetMinimum(0.0);
403 
404  h->GetYaxis()->SetTitle("R_{mT}");
405  h->GetXaxis()->SetLabelSize(0.045);
406  h->Draw();
407 
408  TLegend leg2(legX, legY-legH, legX+legW, legY);
409  leg2.SetTextSize(0.057); leg2.SetFillColor(0); leg2.SetFillStyle(0); leg2.SetBorderSize(0);
410  leg2.SetTextFont(132);
411 
412  graphsMT.at(0)->SetMarkerStyle(20);
413  graphsMT.at(0)->SetMarkerSize(1.2);
414  graphsMT.at(0)->SetMarkerColor(31);
415  graphsMT.at(0)->SetLineColor(31);
416  graphsMT.at(0)->Draw("P");
417  leg2.AddEntry(graphsMT.at(0), sampleName+" "+cuts2title(mj_cut),"p");
418 
419  graphsMT.at(1)->SetMarkerStyle(22);
420  graphsMT.at(1)->SetMarkerSize(1.2);
421  graphsMT.at(1)->SetMarkerColor(46);
422  graphsMT.at(1)->SetLineColor(46);
423  graphsMT.at(1)->Draw("P");
424  leg2.AddEntry(graphsMT.at(1), sampleName+" "+cuts2title(invertcut(mj_cut)),"p");
425 
426  leg2.Draw();
427 
428  text78->Draw();
429  text9->Draw();
430  line.DrawLine((nExt-1)/2., 0, (nExt-1)/2., 1.5*maxMT);
431 
432  TString pname2 = "plots/closure/"+binlabels.at(nExt)+"_"+sampleNameText+"_rmt.pdf";
433  TString pname2root = "plots/closure/"+binlabels.at(nExt)+"_"+sampleNameText+"_rmt.root";
434  if(baseline.Contains("ntks_chg==0")) { pname2.ReplaceAll("_rmt","_rmt_ITV"); pname2root.ReplaceAll("_rmt","_rmt_ITV");}
435  can2.SaveAs(pname2);
436  can2.SaveAs(pname2root);
437 
438  graphsMT.at(0)->Delete();
439  graphsMT.at(1)->Delete();
440 
441  //now plot Kappa
442  TCanvas can3;
443 
444  //make hist to define pad with labels on x-axis
445  h->SetMaximum(2.0);
446  h->GetYaxis()->SetTitle("\\kappa");
447  h->Draw();
448 
449  kappas->SetMarkerStyle(20);
450  kappas->SetMarkerSize(1.2);
451  kappas->SetMarkerColor(kBlack);
452  kappas->SetLineColor(kBlack);
453  kappas->Draw("P");
454  line.DrawLine(h->GetBinLowEdge(1), 1, h->GetBinLowEdge(h->GetNbinsX()+1), 1);
455 
456  // line.DrawLine((nExt-1)/2., 0, (nExt-1)/2., 1.5*maxMT); //vertical line
457  text78->Draw();
458  text9->Draw();
459 
460  TString pname3 = "plots/closure/"+binlabels.at(nExt)+"_"+sampleNameText+"_kappa.pdf";
461  TString pname3root = "plots/closure/"+binlabels.at(nExt)+"_"+sampleNameText+"_kappa.root";
462  if(baseline.Contains("ntks_chg==0")) { pname3.ReplaceAll("_kappa","_kappa_ITV"); pname3root.ReplaceAll("_kappa","_kappa_ITV");}
463  can3.SaveAs(pname3);
464  can3.SaveAs(pname3root);
465 
466  kappas->Delete();
467  h->Delete();
468 }
469 
470 double GetWeight(sfeats Sample, TString luminosity){
471 
472  TChain * chain = new TChain("tree");
473  for(unsigned insam(0); insam < Sample.file.size(); insam++)
474  chain->Add(Sample.file[insam]);
475 
476  TH1D h_weight("weight", "",2000, 0, 2);
477  h_weight.Sumw2();
478 
479  chain->Project("weight", "weight*"+luminosity);
480 
481  double weight = h_weight.GetMean(); // This method could be improved. Only works if all the weights are the same.
482 
483  return weight;
484 }
485 
486 void GetYield(double& yield, double& yield_err, TString baseline,TString external,TString mj_cut,TString mt_cut, sfeats Sample,TString luminosity){
487 
488  TChain * chain = new TChain("tree");
489  for(unsigned insam(0); insam < Sample.file.size(); insam++)
490  chain->Add(Sample.file[insam]);
491 
492  TH1D h_yield("yields", "",100, 0, 10000);
493  h_yield.Sumw2();
494  TH1D h_weight("weight", "",2000, 0, 2);
495  h_weight.Sumw2();
496 
497  TString yieldCut = luminosity+"*weight*("+baseline+"&&"+external+"&&"+mj_cut+"&&"+mt_cut+"&&"+Sample.cut+")";
498 
499  chain->Project("yields", "met", yieldCut);
500  chain->Project("weight", "weight*"+luminosity);
501 
502  double weight = h_weight.GetMean(); // This method could be improved. Only works if all the weights are the same.
503 
504  yield = h_yield.IntegralAndError(0,101,yield_err);
505 
506  if(yield==0) yield_err = weight;
507 }
508 
509 void GetRatio(double& ratio, double& error, double num, double den, double num_err, double den_err, double weight){
510 
511  if(num==0 && den==0){
512  ratio = 0;
513  error = 0;
514  }
515  else if(num==0){ // Error on zero is ~1 event
516  num = 0;
517  num_err = weight;
518  }
519  else if(den==0){ //Pretend as if you saw 1 event.
520  den = weight*1.;
521  den_err = weight;
522  }
523  ratio = num/den;
524  error = sqrt(pow(num_err/den,2)+pow(num/pow(den,2)*den_err,2));
525 }
526 
527 void GetKappa(double& kappa, double& error, double n1, double n2, double n3, double n4, double n1_err, double n2_err, double n3_err, double n4_err, double weight){
528 
529  int nzeroes = 0;
530  if(n1==0) nzeroes++;
531  if(n2==0) nzeroes++;
532  if(n3==0) nzeroes++;
533  if(n4==0) nzeroes++;
534 
535  if(nzeroes>=2){ // If more than 1 zero, don't calculate since one R-factor is indeterminate
536  kappa = 0;
537  error = 0;
538  }
539  else if(n1==0){ // Error on zero is ~1 event
540  n1 = 0;
541  n1_err = weight;
542  }
543  else if(n4==0){ // Error on zero is ~1 event
544  n4 = 0;
545  n4_err = weight;
546  }
547  else if(n2==0){ // Pretend as if you saw 1 event
548  n2 = weight*1;
549  n2_err = weight;
550  }
551  else if(n3==0){ // Pretend as if you saw 1 event
552  n3 = weight*1;
553  n3_err = weight;
554  }
555  kappa = (n4*n1)/(n3*n2);
556  error = sqrt(pow((n4/(n2*n3))*n1_err,2)+pow((n1/(n2*n3))*n4_err,2)+pow(((n1*n4)/(pow(n2,2)*n3))*n2_err,2)+pow(((n1*n4)/(n2*pow(n3,2)))*n3_err,2));
557 }
GetKappa
void GetKappa(double &kappa, double &error, double n1, double n2, double n3, double n4, double n1_err, double n2_err, double n3_err, double n4_err, double weight)
Definition: plot_closure.cxx:527
styles::setDefaultStyle
void setDefaultStyle()
Definition: styles.cpp:36
rpv_bkg_syst.leg
tuple leg
Definition: rpv_bkg_syst.py:278
GetWeight
double GetWeight(sfeats Sample, TString luminosity)
Definition: plot_closure.cxx:470
cuts2title
TString cuts2title(TString title)
Definition: utilities_macros.cpp:640
std
STL namespace.
main
int main()
Definition: plot_closure.cxx:29
rpv::luminosity
TString luminosity
Definition: utilities_macros_rpv.hpp:3
sfeats::file
std::vector< TString > file
Definition: utilities_macros.hpp:81
MakeGraphs
void MakeGraphs(double rmj1[], double rmj2[], double rmj1_err[], double rmj2_err[], double rmt1[], double rmt2[], double rmt1_err[], double rmt2_err[], double kappa[], double kappa_err[], const vector< TString > binlabels, TString sampleName, TString baseline, TString mj_cut, TString mt_cut)
Definition: plot_closure.cxx:291
GetRatio
void GetRatio(double &ratio, double &error, double num, double den, double num_err, double den_err, double weight)
Definition: plot_closure.cxx:509
parse_card.kappa
tuple kappa
Definition: parse_card.py:264
sfeats::cut
TString cut
Definition: utilities_macros.hpp:82
GetYield
void GetYield(double &yield, double &yield_err, TString baseline, TString external, TString mj_cut, TString mt_cut, sfeats Sample, TString luminosity)
Definition: plot_closure.cxx:486
efficiency_map.num
tuple num
Definition: efficiency_map.py:180
utilities_macros.hpp
invertcut
TString invertcut(TString cut)
Definition: utilities_macros.cpp:712
Generated by   doxygen 1.8.11