ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/ZbTools.C
(Generate patch)

Comparing UserCode/cbrown/Development/Plotting/Modules/ZbTools.C (file contents):
Revision 1.14 by buchmann, Mon Nov 26 12:06:26 2012 UTC vs.
Revision 1.24 by buchmann, Fri Jan 25 14:01:03 2013 UTC

# Line 1 | Line 1
1   #include <iostream>
2   #include <vector>
3   #include <sys/stat.h>
4 + #include <math.h>
5  
6   #include <TCut.h>
7   #include <TROOT.h>
# Line 21 | Line 22 | using namespace std;
22   using namespace PlottingSetup;
23  
24   namespace ZbData {
24 vector<float> data_over_mc;
25
25   TCut ZplusBsel("pt1>20&&pt2>20&&mll>60&&mll<120&&id1==id2");
26 < TCut LeadingB("Zb3010_bTagProbCSVBP[0]>0.898");
27 < TCut EtaB("abs(Zb3010_pfJetGoodEta[0])<1.3");
28 < TCut PhiZcut("abs(Zb3010_pfJetDphiZ[0])>2.7");
26 > TCut LeadingB("ZbCHS3010_bTagProbCSVBP[0]>0.898");
27 > TCut EtaB("abs(ZbCHS3010_pfJetGoodEta[0])<1.3");
28 > TCut PhiZcut("abs(ZbCHS3010_pfJetDphiZ[0])>2.7");
29    
30   }
31  
32 + TCut STDWEIGHT;
33 +
34 + class ZbPointResult {
35 + public:
36 +  ZbPointResult(string var, float ptlow, float pthigh, float alphalow, float alphahigh, Value Data, Value MC);
37 +  float PtLow;
38 +  float PtHigh;
39 +  float AlphaLow;
40 +  float AlphaHigh;
41 +  string variable;
42 +  Value Data;
43 +  Value MC;
44 +  Value Ratio;
45 + };
46 +
47 + ZbPointResult::ZbPointResult(string var, float ptlow, float pthigh, float alphalow, float alphahigh, Value data, Value mc) :
48 +  variable(var), PtLow(ptlow), PtHigh(pthigh), AlphaLow(alphalow), AlphaHigh(alphahigh), Data(data), MC(mc) {
49 +    Ratio = data/mc;
50 + }
51 +  
52 + std::ostream &operator<<(std::ostream &ostr, ZbPointResult &b)
53 + {
54 +  ostr << b.PtLow << ";" << b.PtHigh << ";" << b.AlphaLow << ";" << b.AlphaHigh << ";" << b.Data.getValue() << ";" << b.Data.getError() << ";" << b.MC.getValue() << ";" << b.MC.getError() << ";" << b.Ratio.getValue() << ";" << b.Ratio.getError();
55 +  return ostr;
56 + }
57 +
58 +
59 + class ZbResultContainer {
60 + public:
61 +  string name;  
62 +  vector<ZbPointResult> MPF_Rabs_location;
63 +  Value MPF_Result_FaceValue;
64 +  Value Rabs_Result_FaceValue;
65 +  Value MPF_Result_Extrapolated;
66 +  Value Rabs_Result_Extrapolated;
67 +  
68 +  ZbResultContainer(string name);
69 +  ZbResultContainer(const ZbResultContainer &old);
70 + };
71 +
72 +
73 + ZbResultContainer::ZbResultContainer(string _name) {
74 +  name=_name;
75 +  cout<< "Result Container for \"" << name << "\" has been initialized" << endl;
76 +  
77 + }
78 +
79 + std::ostream &operator<<(std::ostream &ostr, ZbResultContainer &b)
80 + {
81 +  ostr << "Results for " << b.name << endl;
82 +  ostr << " MPF and Rabs locations: " << endl;
83 +  for(int i=0;i<b.MPF_Rabs_location.size();i++) ostr << "   " << b.MPF_Rabs_location[i] << endl;
84 +  ostr << " Results: " << endl;
85 +  ostr << "   Face Value: " << endl;
86 +  ostr << "     MPF " << b.MPF_Result_FaceValue << endl;
87 +  ostr << "     Rabs " << b.Rabs_Result_FaceValue << endl;
88 +  ostr << "   Extrapolated: " << endl;
89 +  ostr << "     MPF " << b.MPF_Result_Extrapolated << endl;
90 +  ostr << "             Rabs " << b.Rabs_Result_Extrapolated << endl;
91 +  
92 +  return ostr;
93 + }
94 +
95 +
96 +
97   using namespace ZbData;
98  
99  
# Line 39 | Line 103 | void SpecialCutFlow(string identifier);
103   //*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*            AND THIS LINE             /*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*//
104  
105  
106 < void print_yield(TCut cut) {
106 > class CutYields{
107 > public:
108 >  float data;
109 >  float Zb;
110 >  float Zc;
111 >  float Zl;
112 >  float tt;
113 >  float SiTo;
114 >  float Dibosons;
115 >  float WJets;
116 >  float Other;
117 >  
118 >  CutYields();
119 > };
120 >
121 > CutYields::CutYields() {
122 >  this->data=0;
123 >  this->Zb=0;
124 >  this->Zc=0;
125 >  this->Zl=0;
126 >  this->tt=0;
127 >  this->SiTo=0;
128 >  this->Dibosons=0;
129 >  this->WJets=0;
130 >  this->Other=0;
131 > }
132 >
133 > std::ostream &operator<<(std::ostream &ostr, CutYields y)
134 > {
135 >  return ostr << y.data << ";" << y.Zb << ";" << y.Zc << ";" << y.Zl << ";" << y.tt << ";" << y.SiTo << ";" << y.Dibosons << ";" << y.WJets << ";" << endl;
136 > }
137 >
138 >
139 > CutYields print_yield(TCut cut) {
140 >  
141 >  CutYields Yield;
142    TH1F *data  = allsamples.Draw("data",    "mll",1,0,10000000, "m_{ll} [GeV]", "events", cut,data,luminosity);
143    dout << "  Yields for " << (const char*) cut << endl;
144 <  dout << "     data: " << data->Integral() << endl;
144 > //  dout << "     data: " << data->Integral() << endl;
145 >  
146 >  Yield.data=data->Integral();
147    THStack mcm = allsamples.DrawStack("mc", "mll",1,0,10000000, "m_{ll} [GeV]", "events", cut,mc,luminosity);
148    int nhists = mcm.GetHists()->GetEntries();
48  dout << "Going to loop over " << nhists << " histograms " << endl;
149    TFile *test = new TFile("test.root","RECREATE");
150    mcm.Write();
151    test->Close();
152    for (int i = 0; i < nhists; i++) {
153      TH1* hist = static_cast<TH1*>(mcm.GetHists()->At(i));
154 <    cout << "     " << hist->GetName() << ": " << hist->Integral() << endl;
154 >    string name=hist->GetName();
155 >    if(Contains(name,"t_bar_t")) Yield.tt+=hist->Integral();
156 >    if(Contains(name,"W_Jets")) Yield.WJets+=hist->Integral();
157 >    if(Contains(name,"Single_top")) Yield.SiTo+=hist->Integral();
158 >    if(Contains(name,"Dibosons")) Yield.Dibosons+=hist->Integral();
159 >    if(Contains(name,"Z_l_")) Yield.Zl+=hist->Integral();
160 >    if(Contains(name,"Z_c_")) Yield.Zc+=hist->Integral();
161 >    if(Contains(name,"Z_b_")) Yield.Zb+=hist->Integral();
162    }
163 <    
163 >  
164 >  cout << Yield << endl;
165    cout << endl << endl;
166    
167    delete data;
168 +  CleanLegends();
169 +  DeleteStack(mcm);
170 +
171 +  return Yield;
172   }
173  
174   void draw_kin_variable(string variable, TCut cut, int nbins, float min, float max, string xlabel, string saveas, bool dologscale, string speciallabel="") {
63  
175    TCanvas *ckin = new TCanvas("ckin","kin variable canvas");
176 <  ckin->SetLogy(dologscale);
176 >  if(dologscale) ckin->SetLogy(1);
177    TH1F *datahisto = allsamples.Draw("datahisto",variable,nbins,min,max, xlabel, "events",cut,data,luminosity);
178    datahisto->SetMarkerSize(DataMarkerSize);
179    THStack mcstack = allsamples.DrawStack("mcstack",variable,nbins,min,max, xlabel, "events",cut,mc,luminosity);
180    if(dologscale) datahisto->SetMaximum(3*datahisto->GetMaximum());
181    else datahisto->SetMaximum(1.5*datahisto->GetMaximum());
182 <  TText *speciall;
183 <  if(speciallabel!="") {
184 <    speciall = new TText(0.2,0.8,speciallabel.c_str());
185 <    speciall->Draw();
182 >  
183 >  float SumMC=0.0;
184 >  TIter nextH(mcstack.GetHists());
185 >  TH1F* h;
186 >  while ( h = (TH1F*)nextH() ) {
187 >    SumMC += h->Integral();
188 >  }
189 >  float scale = datahisto->Integral()/SumMC;
190 >
191 >  // Re-scale stack to data integral
192 >  nextH = TIter(mcstack.GetHists());
193 >  
194 >  while ( h = (TH1F*)nextH() ) {
195 >    h->Scale(scale);
196    }
197 +  
198 +  TText *speciall;
199    datahisto->Draw("e1");
200    ckin->Update();
201 <  mcstack.Draw("same");
201 >  mcstack.Draw("histo,same");
202    datahisto->Draw("same,e1");
203    TLegend *kinleg = allsamples.allbglegend();
204 +
205    kinleg->Draw();
206  
207    TPad *kinpad = new TPad("kinpad","kinpad",0,0,1,1);
208    kinpad->cd();
209 <  kinpad->SetLogy(dologscale);
209 >  if(dologscale) kinpad->SetLogy(1);
210    datahisto->Draw("e1");
211 <  mcstack.Draw("same");
211 >  mcstack.Draw("histo,same");
212    datahisto->Draw("same,e1");
213    datahisto->Draw("same,axis");
214    kinleg->Draw();
215    DrawPrelim();
216    saveas="kin/"+saveas;
217 <  save_with_ratio(datahisto,mcstack,kinpad->cd(),saveas);
94 <
95 <  
217 >  kinpad->cd();
218    
219 +  if(speciallabel!="") {
220 +    speciall = write_title(speciallabel);
221 +    speciall->SetX(0.18);
222 +    speciall->SetTextAlign(11);
223 +    speciall->SetTextSize(0.03);
224 +    speciall->SetY(0.85);
225 +    speciall->Draw();
226 +  }
227    
228 <  ZbData::data_over_mc.push_back(datahisto->Integral()/CollapseStack(mcstack)->Integral());
228 >  Save_With_Ratio(datahisto,mcstack,kinpad->cd(),saveas);
229    
230    datahisto->Delete();
231 +  delete kinleg;
232    delete ckin;
233 +  CleanLegends();
234 +  DeleteStack(mcstack);
235   }
236 <  
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 <  
236 >
237   void print_all_b_yields() {
238 +  vector<CutYields> yields;
239    cout << "Basic selection with a b jet" << endl;
240 <  print_yield(ZplusBsel&&"Zb3010_pfJetGoodNumBtag>0");
240 >  yields.push_back(print_yield(ZplusBsel&&"ZbCHS3010_pfJetGoodNumBtag>0"));
241    cout << "Leading jet is a b-jet " << endl;
242 <  print_yield(ZplusBsel&&LeadingB);
242 >  yields.push_back(print_yield(ZplusBsel&&LeadingB));
243    cout << "|#eta_{b}|<1.3 " << endl;
244 <  print_yield(ZplusBsel&&LeadingB&&EtaB);
244 >  yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB));
245    cout << "|#delta#phi(b<Z)|>2.7" << endl;
246 <  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut);
247 <  cout << "10<ptZ<1000 GeV" << endl;
248 <  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000");
246 >  yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut));
247 >  cout << "30<ptZ<1000 GeV" << endl;
248 >  yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"));
249    cout << "#alpha < 0.35" << endl;
250 <  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.35");
250 >  yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"(ZbCHS3010_alpha)<0.35"));
251    cout << "#alpha < 0.3" << endl;
252 <  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.3");
253 <  cout << "#alpha < 0.2" << endl;
254 <  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.2");
252 >  yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"(ZbCHS3010_alpha)<0.3"));
253 > /*  cout << "#alpha < 0.2" << endl;
254 >  yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"(ZbCHS3010_alpha)<0.2"));
255    cout << "#alpha < 0.15" << endl;
256 <  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.15");
256 >  yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"(ZbCHS3010_alpha)<0.15"));
257    cout << "#alpha < 0.1" << endl;
258 <  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.1");
259 < //  cout << "#alpha < 0.05" << endl;
260 < //  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"(Zb3010_alpha)<0.05");
258 >  yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"(ZbCHS3010_alpha)<0.1"));*/
259 >
260 >  for(int iy=0;iy<yields.size();iy++) {
261 >    cout << yields[iy] << endl;
262 >  }
263   }  
264  
265   void DrawEvilCutFlow() {
266    stringstream MegaCut;
267    
268 <  MegaCut << "   1*(" << (const char*) (ZplusBsel&&TCut("Zb3010_pfJetGoodNumBtag>0")) << ")";
268 >  MegaCut << "   1*(" << (const char*) (ZplusBsel&&TCut("ZbCHS3010_pfJetGoodNumBtag>0")) << ")";
269    MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB) << ")";
270    MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB) << ")";
271    MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut) << ")";
272 <  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>10&&pt<1000")) << ")";
273 <  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>10&&pt<1000&&Zb3010_alpha<0.3")) << ")";
274 <  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>10&&pt<1000&&Zb3010_alpha<0.2")) << ")";
275 <  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>10&&pt<1000&&Zb3010_alpha<0.15")) << ")";
276 <  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>10&&pt<1000&&Zb3010_alpha<0.1")) << ")";
277 <  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>10&&pt<1000&&Zb3010_alpha<0.05")) << ")";
272 >  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000")) << ")";
273 >  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.3")) << ")";
274 >  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.2")) << ")";
275 >  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.15")) << ")";
276 >  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.1")) << ")";
277 >  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.05")) << ")";
278    MegaCut << " ";
279    
280    TCanvas *can = new TCanvas("can","can");
# Line 159 | Line 287 | void DrawEvilCutFlow() {
287    
288    
289    TH1F *data  = allsamples.Draw("data",    MegaCut.str(),11,-0.5,10.5, "", "events", basecut,data,luminosity);
290 <   THStack mcm = allsamples.DrawStack("mc", MegaCut.str(),11,-0.5,10.5, "", "events", basecut,mc,luminosity);
290 >  THStack mcm = allsamples.DrawStack("mc", MegaCut.str(),11,-0.5,10.5, "", "events", basecut,mc,luminosity);
291  
292    float runningsum=0;
293    for(int i=data->GetNbinsX();i>0;i--) {
# Line 190 | Line 318 | void DrawEvilCutFlow() {
318    data->GetXaxis()->SetBinLabel(3,"Leading jet is b-jet");
319    data->GetXaxis()->SetBinLabel(4,"|#eta_{b}|<1.3 ");
320    data->GetXaxis()->SetBinLabel(5,"|#delta#phi(b,Z)|>2.7");
321 <  data->GetXaxis()->SetBinLabel(6,"10<p_{T}^{Z}<1000 GeV");
321 >  data->GetXaxis()->SetBinLabel(6,"30<p_{T}^{Z}<1000 GeV");
322    data->GetXaxis()->SetBinLabel(7,"#alpha < 0.3");
323    data->GetXaxis()->SetBinLabel(8,"#alpha < 0.2");
324    data->GetXaxis()->SetBinLabel(9,"#alpha < 0.15");
# Line 199 | Line 327 | void DrawEvilCutFlow() {
327    data->GetXaxis()->LabelsOption("v");
328    
329    data->Draw("e1");
330 <  mcm.Draw("same");
330 >  mcm.Draw("histo,same");
331    data->Draw("same,e1,axis");
332    data->Draw("same,e1");
333   //  data->Draw("same,TEXT");
# Line 209 | Line 337 | void DrawEvilCutFlow() {
337    DrawPrelim();
338    
339    CompleteSave(can,"CutFlow");
340 +  CleanLegends();
341 +  DeleteStack(mcm);
342 +
343   }
344  
345   void draw_Zb_kin_vars() {
346 <   draw_kin_variable("pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",25,0,200,"Z p_{T}","Official/Zpt",1);
347 <   draw_kin_variable("Zb3010_pfJetGoodPt[1]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",25,0,200,"Sub-Leading Jet Pt","Official/Jet2Pt",1);
348 <   draw_kin_variable("Zb3010_pfJetGoodPt[0]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",25,0,200,"Leading Jet Pt (B)","Official/LeadingJetPt",1);
349 <   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);
350 <   draw_kin_variable("Zb3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<50",   40,0,2,"#alpha","Official/alpha_pt_30_to_50",1);
351 <   draw_kin_variable("Zb3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>50&&pt<75",  40,0,2,"#alpha","Official/alpha_50_to_75",1);
352 <   draw_kin_variable("Zb3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>75&&pt<125", 40,0,2,"#alpha","Official/alpha_pt_75_to_125",1);
353 <   draw_kin_variable("Zb3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>125&&pt<1000",40,0,2,"#alpha","Official/alpha_pt_125_to_1000",1);
354 <    draw_kin_variable("Zb3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000",40,0,2,"#alpha","Official/alpha_pt_30_to_1000",1);
355 <  
356 <  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);
357 <  
227 < /*
228 <  draw_kin_variable("pt1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,0,200,"p_{T}^{l1}","Lep/pt1",1);
229 <  draw_kin_variable("pt1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id1==0",20,0,200,"p_{T}^{e1}","Lep/pt1_e",1);
230 <  draw_kin_variable("pt1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id1==1",20,0,200,"p_{T}^{#mu1}","Lep/pt1_m",1);
231 <  draw_kin_variable("pt2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,0,200,"p_{T}^{l2}","Lep/pt2",1);
232 <  draw_kin_variable("pt2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id2==0",20,0,200,"p_{T}^{e2}","Lep/pt2_e",1);
233 <  draw_kin_variable("pt2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000&&id2==1",20,0,200,"p_{T}^{#mu2}","Lep/pt2_m",1);
234 <  draw_kin_variable("eta1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,-3.1,3.1,"#eta^{l1}","Lep/eta1",1);
235 <  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);
236 <  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);
237 <  draw_kin_variable("eta2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,-3.1,3.1,"#eta^{l1}","Lep/eta2",1);
238 <  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);
239 <  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);
240 <  draw_kin_variable("numVtx",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",30,0,30,"N(Vtx)","Lep/nVtx",1);
241 <  */
346 >  
347 >   draw_kin_variable("pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000"),25,0,200,"Z p_{T}","Official/Zpt",1);
348 >   draw_kin_variable("ZbCHS3010_pfJetGoodPt[1]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000"),25,0,200,"Sub-Leading Jet Pt","Official/Jet2Pt",1);
349 >   draw_kin_variable("ZbCHS3010_pfJetGoodPt[0]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000"),25,0,200,"Leading Jet Pt (B)","Official/LeadingJetPt",1);
350 >   draw_kin_variable("ZbCHS3010_pfJetGoodPt[1]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000"),20,0,2,"p_{T}^{J2} / p_{T}^{Z}","Official/SubLeadingJetPt_Over_ZPt",1);
351 >   draw_kin_variable("ZbCHS3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<50"),   40,0,2,"#alpha","Official/alpha_pt_30_to_50",1);
352 >   draw_kin_variable("ZbCHS3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>50&&pt<75"),  40,0,2,"#alpha","Official/alpha_50_to_75",1);
353 >   draw_kin_variable("ZbCHS3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>75&&pt<125"), 40,0,2,"#alpha","Official/alpha_pt_75_to_125",1);
354 >   draw_kin_variable("ZbCHS3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>125&&pt<1000"),40,0,2,"#alpha","Official/alpha_pt_125_to_1000",1);
355 >   draw_kin_variable("ZbCHS3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000"),40,0,2,"#alpha","Official/alpha_pt_30_to_1000",1);
356 >  
357 >   draw_kin_variable("ZbCHS3010_pfBJetDphiZ[0]",ZplusBsel&&LeadingB&&TCut("pt>30&&pt<1000&&pfBJetDphiZ[0]>0"),30,0,3.2,"#delta#phi (Z,b lead)","DeltaPhiZBlead",0);
358   }
359  
360   void draw_mpf_vars() {
361 <  ZbData::data_over_mc.clear();
362 <  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");
363 <  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");
364 <  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");
365 <  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");
366 <  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");
367 <  
368 <  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");
369 <  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");
370 <  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");
371 <  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");
256 <  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");
361 >  draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.3"), 20,0.0,2.0,"Z+b MPF","mpf_alpha_smaller_0p3",0,"#alpha<0.3, 30 GeV < p_{T}^{Z} < 1000 GeV");
362 >  draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.2"), 20,0.0,2.0,"Z+b MPF","mpf_alpha_smaller_0p2",0,"#alpha<0.2, 30 GeV < p_{T}^{Z} < 1000 GeV");
363 >  draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.15"),20,0.0,2.0,"Z+b MPF","mpf_alpha_smaller_0p15",0,"#alpha<0.15, 30 GeV < p_{T}^{Z} < 1000 GeV");
364 >  draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.1"), 20,0.0,2.0,"Z+b MPF","mpf_alpha_smaller_0p1",0,"#alpha<0.1, 30 GeV < p_{T}^{Z} < 1000 GeV");
365 >  draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.05"),20,0.0,2.0,"Z+b MPF","mpf_alpha_smaller_0p05",0,"#alpha<0.05, 30 GeV < p_{T}^{Z} < 1000 GeV");
366 >  
367 >  draw_kin_variable("ZbCHS3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.3"), 20,0.0,2.0,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p3",    0,"#alpha<0.3, 30 GeV < p_{T}^{Z} < 1000 GeV");
368 >  draw_kin_variable("ZbCHS3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.2"), 20,0.0,2.0,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p2",    0,"#alpha<0.2, 30 GeV < p_{T}^{Z} < 1000 GeV");
369 >  draw_kin_variable("ZbCHS3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.15"),20,0.0,2.0,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p15",   0,"#alpha<0.15, 30 GeV < p_{T}^{Z} < 1000 GeV");
370 >  draw_kin_variable("ZbCHS3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.1"), 20,0.0,2.0,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p1",    0,"#alpha<0.1, 30 GeV < p_{T}^{Z} < 1000 GeV");
371 >  draw_kin_variable("ZbCHS3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.05"),20,0.0,2.0,"p_{T} b-jet / p_{T} Z","ptb_over_ptz__mpf_alpha_smaller_0p05",0,"#alpha<0.05, 30 GeV < p_{T}^{Z} < 1000 GeV");
372   }  
373  
374 < /*double GetMedian(TH1F *histo) {
260 <    int numBins = histo->GetXaxis()->GetNbins();
261 <    Double_t x[numBins];
262 <    Double_t y[numBins];
263 <    for (int i = 1; i <= numBins; i++) {
264 <        x[i] = histo->GetBinCenter(i);
265 <        y[i] = histo->GetBinContent(i);
266 <    }
267 <    return TMath::Median(numBins, &x, &y);
268 < }*/
269 <
270 < Value get_Zb_data_over_mc(string variable, string variable2, TCut cut, string saveas) {
374 > Value get_Zb_data_over_mc(string ContainerName, ZbResultContainer &res, string variable, string variable2, TCut cut, string saveas, float pt1, float pt2, float a1, float a2) {
375   //  write_warning(__FUNCTION__,"Debugging this function, therefore always returning 3.1415+/-1.010101"); return Value(3.1415,1.010101);
376    TCanvas *cn = new TCanvas("cn","cn");
377    string varname="MPF";
274  if(Contains(variable,"pfJetGood")) varname="R_{abs}";
275  TH1F *hdata  = allsamples.Draw("hdata",variable,1200,0,30, varname, "events",  cut,data,luminosity);
276  TH1F *hmc    = allsamples.Draw("hmc" , variable,1200,0,30, varname, "events", cut, mc,  luminosity);
277  hmc->SetLineColor(kBlue);
378    
379 +  bool DoFit=false;
380 +  bool UseFit=false;
381    
382    
383 +  if(Contains(variable,"pfJetGood")) varname="R_{abs}";
384 + //  TH1F *hdata  = allsamples.Draw("hdata",variable,1200,0,30, varname, "events",  cut,data,luminosity);
385 + //  TH1F *hmc    = allsamples.Draw("hmc" , variable,1200,0,30, varname, "events", cut, mc,  luminosity);
386 +  TH1F *hdata  = allsamples.Draw("hdata",variable,1200,0.5,1.5, varname, "events",  cut,data,luminosity); // making it more precise
387 +  TH1F *hmc    = allsamples.Draw("hmc" , variable,1200,0.5,1.5, varname, "events", cut, mc,  luminosity); // making it more precise
388 +  hmc->SetLineColor(kBlue);
389    
390    float a=hdata->GetMean();
391    float b=hmc->GetMean();
# Line 287 | Line 395 | Value get_Zb_data_over_mc(string variabl
395    float error = TMath::Sqrt(  (1/(b*b)) * (da*da) + ((a*a)/(b*b))*db*db);
396    
397    TF1 *gaus_data = new TF1("gaus_data","gaus",0.7,1.3);
290  gaus_data->SetParameter(0,hdata->GetMaximum());
291  gaus_data->SetParameter(2,hdata->GetMean());
292  gaus_data->SetParameter(1,hdata->GetRMS());
293  hdata->Fit(gaus_data);
398    TF1 *gaus_mc = new TF1("gaus_mc","gaus",0.7,1.3);
399 <  gaus_mc->SetParameter(0,hmc->GetMaximum());
400 <  gaus_mc->SetParameter(1,hmc->GetMean());
401 <  gaus_mc->SetParameter(2,hmc->GetRMS());
402 <  hmc->Fit(gaus_mc);
403 <  
399 >  if(DoFit) {
400 >    gaus_data->SetParameter(0,hdata->GetMaximum());
401 >    gaus_data->SetParameter(2,hdata->GetMean());
402 >    gaus_data->SetParameter(1,hdata->GetRMS());
403 >    hdata->Fit(gaus_data);
404 >    
405 >    gaus_mc->SetParameter(0,hmc->GetMaximum());
406 >    gaus_mc->SetParameter(1,hmc->GetMean());
407 >    gaus_mc->SetParameter(2,hmc->GetRMS());
408 >    hmc->Fit(gaus_mc);
409 >  }
410    
411    cn->cd();
412 <  hdata->GetYaxis()->SetTitle("events / 0.1");
412 >  hdata->Rebin(0.05*1200/(1.5-0.5));
413 >  hmc->Rebin(0.05*1200/(1.5-0.5));
414 >  hdata->GetYaxis()->SetTitle("events / 0.05");
415    hdata->SetMarkerColor(kRed);
416 < /*  hdata->Rebin(4);
417 <  hmc->Rebin(4);
418 <  
419 <  gaus_mc->SetParameter(0,gaus_mc->GetParameter(0)*4);//rebinning
308 <  gaus_data->SetParameter(0,gaus_data->GetParameter(0)*4);//rebinning
309 < */
310 <  gaus_mc->SetLineColor(kBlue);
311 <  gaus_data->SetLineColor(kRed);
416 >  if(DoFit) {
417 >    gaus_mc->SetLineColor(kBlue);
418 >    gaus_data->SetLineColor(kRed);
419 >  }
420    
421    hdata->GetXaxis()->SetRangeUser(0,2);
422 +  hdata->GetYaxis()->SetTitleOffset(1.3);
423    hdata->Draw("e1");
424 +  
425 +  hmc->Scale(hdata->Integral()/hmc->Integral());
426 +  if(hmc->GetMaximum()>hdata->GetMaximum()) hdata->SetMaximum((hmc->GetMaximum()+sqrt((double)hmc->GetMaximum()))*1.2);
427 +  
428    hmc->Draw("histo,same");
429    hdata->Draw("e1,same");
430 <  gaus_data->Draw("same");
431 <  gaus_mc->Draw("same");
430 >  if(DoFit) {
431 >    gaus_data->Draw("same");
432 >    gaus_mc->Draw("same");
433 >  }
434 >  
435    stringstream summary;
436 <  summary << "#splitline{" << varname << "}{#splitline{data: " << std::setprecision(4) << a << " #pm " << std::setprecision(4) << da << "}{";
436 >  summary << "#splitline{#bf{" << varname << "}}{#splitline{" << pt1 << " < p_{T}^{Z} < " << pt2 << ", " << a1 << " #leq #alpha < " << a2 << "}{#splitline{data: " << std::setprecision(4) << a << " #pm " << std::setprecision(4) << da << "}{";
437    summary << "#splitline{mc:   " << std::setprecision(4) << b << " #pm " << std::setprecision(4) << db << "}{#splitline{ratio: " << factor;
438 <  summary << " #pm " << error << "}{#splitline{}{#splitline{data fit: " << std::setprecision(4) << gaus_data->GetParameter(1) << "+/-";
439 <  summary << std::setprecision(4) << gaus_data->GetParError(1) << "}{#splitline{";
440 <  summary << "mc fit:" << std::setprecision(4) << gaus_mc->GetParameter(1) << "+/-" << std::setprecision(4) << gaus_mc->GetParError(1) << "}{ratio: ";
441 <  summary << gaus_data->GetParameter(1)/gaus_mc->GetParameter(1) << "+/-" << gaus_data->GetParameter(1)/gaus_mc->GetParameter(1) * sqrt((gaus_data->GetParError(1) * gaus_data->GetParError(1))/(gaus_data->GetParameter(1)*gaus_data->GetParameter(1)) + (gaus_mc->GetParError(1) * gaus_mc->GetParError(1))/(gaus_mc->GetParameter(1)*gaus_mc->GetParameter(1)) ) << "}}}}}}}";
438 >  summary << " #pm " << error;
439 >  if(DoFit) {
440 >    summary << "}{#splitline{}{#splitline{data fit: " << std::setprecision(4) << gaus_data->GetParameter(1) << "+/-";
441 >    summary << std::setprecision(4) << gaus_data->GetParError(1) << "}{#splitline{";
442 >    summary << "mc fit:" << std::setprecision(4) << gaus_mc->GetParameter(1) << "#pm" << std::setprecision(4) << gaus_mc->GetParError(1) << "}{ratio: ";
443 >    summary << gaus_data->GetParameter(1)/gaus_mc->GetParameter(1) << "#pm" << gaus_data->GetParameter(1)/gaus_mc->GetParameter(1) * sqrt((gaus_data->GetParError(1) * gaus_data->GetParError(1))/(gaus_data->GetParameter(1)*gaus_data->GetParameter(1)) + (gaus_mc->GetParError(1) * gaus_mc->GetParError(1))/(gaus_mc->GetParameter(1)*gaus_mc->GetParameter(1)) ) << "}}}}";
444 >    summary << "}}}}";
445 >  } else {
446 >    summary << "}{}}}}}";
447 >  }
448 >  
449 >  if(UseFit) {
450 >    write_error(__FUNCTION__,"Requested using fit result for final result - not implemented yet!");
451 >    assert(0);
452 >  }
453 >    
454 >  
455    TText *infobox = write_title(summary.str());
456    infobox->SetX(0.75);
457    infobox->SetTextSize(0.03);
# Line 331 | Line 460 | Value get_Zb_data_over_mc(string variabl
460      
461      
462    DrawPrelim();
463 <  CompleteSave(cn,"ResponsePlots/"+saveas);
463 >  CompleteSave(cn,ContainerName+"/ResponsePlots/"+saveas);
464    
465 <  cout << "Data;" << a << ";"<<da<<";MC"<<b<<";"<<db<<endl;
465 >  res.MPF_Rabs_location.push_back(ZbPointResult(variable,pt1,pt2,a1,a2,Value(a,da),Value(b,db)));
466 >  cout << "Data;" << a << ";"<<da<<";MC;"<<b<<";"<<db<<endl;
467    delete cn;
468    delete hdata;
469    delete hmc;
470  
341 //    cout << a << "+/-" << da << "        " << b << "=/-" << db << endl;
342
471    return Value(factor,error);
472   }
473    
474 + int nFakes=1;
475 +
476 + float PurityInclusive=0.05;
477 + float PurityLoose=0.20;
478 + float PurityMedium=0.54;
479 + float PurityTight=0.82;
480 +
481 +
482 + void DeterminePurity(TCut basecut, int nptcuts, float ptcuts[], string AlphaVariable, string ContainerName, float FaceValueToBeTakenAt) {
483 +  
484 +  bool ComputePurity=false;
485 +  
486 +  if(ContainerName=="inclusive") ComputePurity=true;  
487 +  if(ContainerName=="loose") ComputePurity=true;
488 +  if(ContainerName=="medium") ComputePurity=true;
489 +  if(ContainerName=="reference") ComputePurity=true;
490 +  
491 +  if(!ComputePurity) return;
492 +  
493 +  stringstream specialcut;
494 +  specialcut << "((pt>30&& pt< 1000) && (" << AlphaVariable << "<" << FaceValueToBeTakenAt << "))";
495 +  
496 +  THStack mcm    = allsamples.DrawStack("mcm" , "mpf",1200,0,30, "MPF", "events", TCut(basecut && specialcut.str().c_str()), mc,  luminosity);
497 +  int nhists = mcm.GetHists()->GetEntries();
498 +  
499 +  float total=0;
500 +  float Zb=0;
501 +  
502 +  for (int i = 0; i < nhists; i++) {
503 +    TH1* hist = static_cast<TH1*>(mcm.GetHists()->At(i));
504 +    string name=hist->GetName();
505 +    if(Contains(name,"t_bar_t")) total+=hist->Integral();
506 +    if(Contains(name,"W_Jets")) total+=hist->Integral();
507 +    if(Contains(name,"Single_top")) total+=hist->Integral();
508 +    if(Contains(name,"Dibosons")) total+=hist->Integral();
509 +    if(Contains(name,"Z_l_")) total+=hist->Integral();
510 +    if(Contains(name,"Z_c_")) total+=hist->Integral();
511 +    if(Contains(name,"Z_b_")) {
512 +      total+=hist->Integral();
513 +      Zb+=hist->Integral();
514 +    }
515 +  }
516    
517 +  float purity=Zb/total;
518 +  dout << "Determined a purity of " << 100*purity << " % for the container " << ContainerName << endl;
519 +  
520 +  if(ContainerName=="inclusive") PurityInclusive=purity;
521 +  if(ContainerName=="loose") PurityLoose=purity;
522 +  if(ContainerName=="medium") PurityMedium=purity;
523 +  if(ContainerName=="reference") PurityTight=purity;
524 +  
525 +  DeleteStack(mcm);
526 + }
527 +
528 + string NiceAlphaPrint(float alpha) {
529 +  //get in : 0.3, return 0p3 ...
530 +  stringstream res;
531 +  res << int(alpha) << "p" << int(100*alpha);
532 +  return res.str();
533 + }
534  
535 < void new_data_mc_agreement_2d() {
535 > ZbResultContainer new_data_mc_agreement_2d(bool gofast, string AlphaVariable, string ContainerName, TCut BTagCut, TCut BTagWgt, float FaceValueToBeTakenAt) {
536 >  
537 >  
538 >  if(gofast) write_info(__FUNCTION__,"Fast mode activated for "+ContainerName);
539 >  ZbResultContainer results(ContainerName);
540 >  
541 >  cutWeight=STDWEIGHT*BTagWgt;
542 >  LeadingB=BTagCut;
543 >  
544    
545    gStyle->SetOptFit(0);
546    TCut basecut(ZplusBsel&&LeadingB&&EtaB&&PhiZcut);
547    const int nalphacuts=6;
548    float alphacuts[nalphacuts] = {0.1,0.15,0.2,0.25,0.3,0.35};
549    
550 +  
551    const int nptcuts=5;
552    float ptcuts[nptcuts]={30,50,75,125,1000};
553    
554 +  DeterminePurity(basecut, nptcuts,ptcuts, AlphaVariable, ContainerName, FaceValueToBeTakenAt);
555 +  
556    float MPF_Results[nalphacuts][nptcuts];
557    float MPF_Errors[nalphacuts][nptcuts];
558    float RABS_Results[nalphacuts][nptcuts];
559    float RABS_Errors[nalphacuts][nptcuts];
560    
561    
562 <  for(int ia=0;ia<nalphacuts;ia++) {
562 >  for(int ia=0;ia<nalphacuts && !gofast;ia++) {
563      for(int ipt=0;ipt<nptcuts-1;ipt++) {
564        stringstream specialcut;
565 <      if(ia>0) specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (Zb3010_alpha<" << alphacuts[ia] << " && Zb3010_alpha>" << alphacuts[ia-1] << "))";
566 <      else     specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (Zb3010_alpha<" << alphacuts[ia] << "))";
565 >      float alow,ahigh;
566 >      if(ia>0) {
567 >        specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (" << AlphaVariable << "<" << alphacuts[ia] << " && " << AlphaVariable << ">" << alphacuts[ia-1] << "))";
568 >        alow=alphacuts[ia-1];
569 >        ahigh=alphacuts[ia];
570 >      } else {
571 >        specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (" << AlphaVariable << "<" << alphacuts[ia] << "))";
572 >        alow=0;
573 >        ahigh=alphacuts[ia];
574 >      }
575        cout << specialcut.str() << endl;
576 <      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]));
577 <      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]));
576 >      Value MPF_data_over_mc = get_Zb_data_over_mc(ContainerName,results,"mpf","",TCut(basecut && specialcut.str().c_str()),"MPF___pt_"+any2string(ptcuts[ipt])+"_"+any2string(ptcuts[ipt+1])+"__alpha_"+any2string(NiceAlphaPrint(alphacuts[ia])), ptcuts[ipt], ptcuts[ipt+1],alow,ahigh);
577 >      Value RABS_data_over_mc = get_Zb_data_over_mc(ContainerName,results,"ZbCHS3010_pfJetGoodPt[0]/pt","",TCut(basecut && specialcut.str().c_str()),"Rabs___pt_"+any2string(ptcuts[ipt])+"_"+any2string(ptcuts[ipt+1])+"__alpha_"+any2string(NiceAlphaPrint(alphacuts[ia])), ptcuts[ipt], ptcuts[ipt+1],alow,ahigh);
578        
579        MPF_Results[ia][ipt]=MPF_data_over_mc.getValue();
580        MPF_Errors[ia][ipt]=MPF_data_over_mc.getError();
# Line 402 | Line 608 | void new_data_mc_agreement_2d() {
608      int rabs_counter=0;
609      int mpf_counter=0;
610      
611 <    for(int ia=0;ia<nalphacuts;ia++) {
611 >    for(int ia=0;ia<nalphacuts && !gofast;ia++) {
612        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;
613        if(MPF_Results[ia][ipt]==MPF_Results[ia][ipt] && MPF_Errors[ia][ipt]==MPF_Errors[ia][ipt]) { // remove nan's
614         mpf_gr->SetPoint(mpf_counter,alphacuts[ia],MPF_Results[ia][ipt]);
# Line 425 | Line 631 | void new_data_mc_agreement_2d() {
631      }
632      
633      stringstream specialcut;
634 <    specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (Zb3010_alpha<0.3))";
635 <    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])+"__INCLUSIVEalpha_0.3");
636 <    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])+"__INCLUSIVEalpha_0.3");
634 >    specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (" << AlphaVariable << "<" << FaceValueToBeTakenAt << "))";
635 >    Value MPF_data_over_mc = get_Zb_data_over_mc(ContainerName,results,"mpf","",TCut(basecut && specialcut.str().c_str()),"MPF___pt_"+any2string(ptcuts[ipt])+"_"+any2string(ptcuts[ipt+1])+"__INCLUSIVEalpha_"+any2string(FaceValueToBeTakenAt),ptcuts[ipt],ptcuts[ipt+1],0,FaceValueToBeTakenAt);
636 >    Value RABS_data_over_mc = get_Zb_data_over_mc(ContainerName,results,"ZbCHS3010_pfJetGoodPt[0]/pt","",TCut(basecut && specialcut.str().c_str()),"Rabs___pt_"+any2string(ptcuts[ipt])+"_"+any2string(ptcuts[ipt+1])+"__INCLUSIVEalpha_"+any2string(FaceValueToBeTakenAt),ptcuts[ipt],ptcuts[ipt+1],0,FaceValueToBeTakenAt);
637      cout << "About to add a point in pt (" << (ptcuts[ipt]+ptcuts[ipt+1])/2 << " ) to face value : " << MPF_data_over_mc << " (MPF) , " <<  RABS_data_over_mc << "(RABS)" << endl;
638      MPF_FaceValueAtPoint3->SetPoint(ipt,(ptcuts[ipt]+ptcuts[ipt+1])/2,MPF_data_over_mc.getValue());
639      MPF_FaceValueAtPoint3->SetPointError(ipt,0,MPF_data_over_mc.getError());
640      RABS_FaceValueAtPoint3->SetPoint(ipt,(ptcuts[ipt]+ptcuts[ipt+1])/2,RABS_data_over_mc.getValue());
641      RABS_FaceValueAtPoint3->SetPointError(ipt,0,RABS_data_over_mc.getError());
642      
437    can->cd();
438    mpf_gr->SetMarkerStyle(21);
439    mpf_gr->SetMarkerColor(kBlue);
440    mpf_gr->Draw("AP");
441    mpf_gr->Fit("pol1","");
442    mpf_gr->GetXaxis()->SetTitle("#alpha");
443    mpf_gr->GetXaxis()->CenterTitle();
444    mpf_gr->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
445    mpf_gr->GetYaxis()->CenterTitle();
446    TF1 *mpf_pol=(TF1*)mpf_gr->GetFunction("pol1");
447    mpf_pol->SetLineWidth(0);
448    TF1 *mpf_pol_complete = new TF1("mpf_pol_complete","[0]+[1]*x");
449    mpf_pol_complete->SetParameters(mpf_pol->GetParameters());
450    mpf_pol_complete->SetLineWidth(1);
451    mpf_pol_complete->SetLineColor(kBlue);
452    float min,max;
453    FindMinMax(mpf_gr,min,max);
454    TH1F *histo = new TH1F("histo","histo",1,0,0.4);
455    histo->GetXaxis()->SetTitle("#alpha");
456    histo->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
457    histo->GetXaxis()->CenterTitle();
458    histo->GetYaxis()->CenterTitle();
459    histo->SetStats(0);
460    histo->GetXaxis()->SetRangeUser(0,0.4);
461    histo->GetYaxis()->SetRangeUser(min,max);
462    mpf_gr->GetYaxis()->SetRangeUser(min,max);
463    
464    TPolyLine *mpf_fit_uncert = GetFitUncertaintyShape(mpf_gr, "pol1", min, max,0.0,0.4);
465    mpf_pol->SetLineColor(TColor::GetColor("#0101DF"));
466    mpf_fit_uncert->SetFillColor(TColor::GetColor("#5353FC"));
467    histo->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
468    histo->Draw();
469    mpf_fit_uncert->Draw("F");
470    mpf_fit_uncert->Draw("");
471    mpf_gr->Draw("P");
472    mpf_pol_complete->Draw("same");
473    
474    float mpf_a=mpf_pol->GetParameter(0);
475    float mpf_b=mpf_pol->GetParameter(1);
476    MPF_FinalGraph->SetPoint(ipt,0.5*(ptcuts[ipt]+ptcuts[ipt+1]),mpf_a);
477    MPF_FinalGraph->SetPointError(ipt,0,mpf_pol->GetParError(0));
478    
479    MPF_ExtraPolatedResults[ipt]=mpf_a;
480    
481    stringstream mpf_info;
482    mpf_info <<  "#splitline{#splitline{#splitline{" << ptcuts[ipt] << " GeV < p_{T}^{Z} < " << ptcuts[ipt+1] << " GeV }{ (MPF) }}{Extrapolated value: " << std::setprecision(4) << mpf_a << "}}{#chi^{2}/NDF : " << mpf_pol->GetChisquare() << " / " << mpf_pol->GetNDF() << "}";
483    
484    TText *mpf_mark = write_title(mpf_info.str());
485    mpf_mark->SetX(0.75);
486    mpf_mark->SetTextSize(0.03);
487    mpf_mark->SetY(0.75);
488    mpf_mark->Draw();
643      
490    string filenamebkp=filename.str();
491    filename << "__MPF";
492    DrawPrelim();
493    CompleteSave(can,filename.str());
494    cout << "MPF : " << mpf_a << " + " << mpf_b << " * alpha " << endl;
644      
645 <    filename.str("");
645 >    if(!gofast) {
646 >      can->cd();
647 >      mpf_gr->SetMarkerStyle(21);
648 >      mpf_gr->SetMarkerColor(kBlue);
649 >      mpf_gr->Draw("AP");
650 >      mpf_gr->Fit("pol1","");
651 >      mpf_gr->GetXaxis()->SetTitle("#alpha");
652 >      mpf_gr->GetXaxis()->CenterTitle();
653 >      mpf_gr->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
654 >      mpf_gr->GetYaxis()->CenterTitle();
655 >      TF1 *mpf_pol=(TF1*)mpf_gr->GetFunction("pol1");
656 >      mpf_pol->SetLineWidth(0);
657 >      TF1 *mpf_pol_complete = new TF1("mpf_pol_complete","[0]+[1]*x");
658 >      mpf_pol_complete->SetParameters(mpf_pol->GetParameters());
659 >      mpf_pol_complete->SetLineWidth(1);
660 >      mpf_pol_complete->SetLineColor(kBlue);
661 >      
662 >      float min,max;
663 >      FindMinMax(mpf_gr,min,max);
664 >      TH1F *histo = new TH1F("histo","histo",1,0,0.4);
665 >      histo->GetXaxis()->SetTitle("#alpha");
666 >      histo->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
667 >      histo->GetXaxis()->CenterTitle();
668 >      histo->GetYaxis()->CenterTitle();
669 >      histo->GetYaxis()->SetTitleOffset(1.3);
670 >      histo->SetStats(0);
671 >      histo->GetXaxis()->SetRangeUser(0,0.4);
672 >      histo->GetYaxis()->SetRangeUser(min,max);
673 >      mpf_gr->GetYaxis()->SetRangeUser(min,max);
674 >      
675 >      TPolyLine *mpf_fit_uncert = GetFitUncertaintyShape(mpf_gr, "pol1", min, max,0.0,0.4);
676 >      mpf_pol->SetLineColor(TColor::GetColor("#0101DF"));
677 >      mpf_fit_uncert->SetFillColor(TColor::GetColor("#5353FC"));
678 >      histo->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
679 >      histo->Draw();
680 >      mpf_fit_uncert->Draw("F");
681 >      mpf_fit_uncert->Draw("");
682 >      mpf_gr->Draw("P");
683 >      mpf_pol_complete->Draw("same");
684 >      
685 >      float mpf_a=mpf_pol->GetParameter(0);
686 >      float mpf_b=mpf_pol->GetParameter(1);
687 >      MPF_FinalGraph->SetPoint(ipt,0.5*(ptcuts[ipt]+ptcuts[ipt+1]),mpf_a);
688 >      MPF_FinalGraph->SetPointError(ipt,0,mpf_pol->GetParError(0));
689 >      
690 >      MPF_ExtraPolatedResults[ipt]=mpf_a;
691 >      
692 >      stringstream mpf_info;
693 >      mpf_info <<  "#splitline{#splitline{#splitline{" << ptcuts[ipt] << " GeV < p_{T}^{Z} < " << ptcuts[ipt+1] << " GeV }{ (MPF) }}{Extrapolated value: " << std::setprecision(4) << mpf_a << "}}{#chi^{2}/NDF : " << mpf_pol->GetChisquare() << " / " << mpf_pol->GetNDF() << "}";
694 >      
695 >      TText *mpf_mark = write_title(mpf_info.str());
696 >      mpf_mark->SetX(0.75);
697 >      mpf_mark->SetTextSize(0.03);
698 >      mpf_mark->SetY(0.75);
699 >      mpf_mark->Draw();
700 >      
701 >      string filenamebkp=filename.str();
702 >      filename << "__MPF";
703 >      DrawPrelim();
704 >      CompleteSave(can,ContainerName+"/"+filename.str());
705 >      cout << "MPF : " << mpf_a << " + " << mpf_b << " * alpha " << endl;
706 >      
707 >      filename.str("");
708 >      
709 >      rabs_gr->SetMarkerStyle(21);
710 >      rabs_gr->SetMarkerColor(kRed);
711 >      rabs_gr->Draw("AP");
712 >      rabs_gr->Fit("pol1","");
713 >      rabs_gr->GetXaxis()->SetTitle("#alpha");
714 >      rabs_gr->GetXaxis()->CenterTitle();
715 >      rabs_gr->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
716 >      rabs_gr->GetYaxis()->CenterTitle();
717 >      TF1 *rabs_pol=(TF1*)rabs_gr->GetFunction("pol1");
718 >      rabs_pol->SetLineWidth(0);
719 >      TF1 *rabs_pol_complete = new TF1("rabs_pol_complete","[0]+[1]*x");
720 >      rabs_pol_complete->SetParameters(rabs_pol->GetParameters());
721 >      rabs_pol_complete->SetLineColor(kRed);
722 >      rabs_pol_complete->SetLineWidth(1);
723 >      FindMinMax(rabs_gr,min,max);
724 >      rabs_gr->GetYaxis()->SetRangeUser(min,max);
725 >      histo->GetYaxis()->SetRangeUser(min,max);
726 >      TPolyLine *rabs_fit_uncert = GetFitUncertaintyShape(rabs_gr, "pol1", min, max,0.0,0.4);
727 >      rabs_pol->SetLineColor(TColor::GetColor("#FA5858"));
728 >      rabs_pol->SetFillColor(TColor::GetColor("#FA5858"));
729 >      rabs_fit_uncert->SetLineColor(TColor::GetColor("#FA5858"));
730 >      rabs_fit_uncert->SetFillColor(TColor::GetColor("#FA5858"));
731 >      histo->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
732 >      histo->Draw();
733 >      rabs_fit_uncert->Draw("F");
734 >      rabs_fit_uncert->Draw("");
735 >      rabs_gr->Draw("P");
736 >      rabs_pol_complete->Draw("same");
737 >      
738 >      float rabs_a=rabs_pol->GetParameter(0);
739 >      float rabs_b=rabs_pol->GetParameter(1);
740 >      RABS_FinalGraph->SetPoint(ipt,0.5*(ptcuts[ipt]+ptcuts[ipt+1]),rabs_a);
741 >      RABS_FinalGraph->SetPointError(ipt,0,rabs_pol->GetParError(0));
742 >      
743 >      cout << "!+!+!+!+!+!+!+!+!+ Just added a point to the final plot : " << rabs_a << " +/- " << rabs_pol->GetParError(0) << endl;
744 >      
745 >      stringstream rabs_info;
746 >      rabs_info << "#splitline{#splitline{#splitline{" << ptcuts[ipt] << " GeV < p_{T}^{Z} < " << ptcuts[ipt+1] << " GeV }{ (R_{abs}) }}{Extrapolated value: " << std::setprecision(4) << rabs_a << "}}{#chi^{2}/NDF : " << rabs_pol->GetChisquare() << " / " << rabs_pol->GetNDF() << "}";
747 >      TText *rabs_mark = write_title(rabs_info.str());
748 >      rabs_mark->SetX(0.75);
749 >      rabs_mark->SetTextSize(0.03);
750 >      rabs_mark->SetY(0.75);
751 >      rabs_mark->Draw();
752 >      
753 >      RABS_ExtraPolatedResults[ipt]=rabs_a;
754 >      filename << filenamebkp << "__RABS";
755 >      DrawPrelim();
756 >      CompleteSave(can,ContainerName+"/"+filename.str());
757 >      cout << "RABS : " << rabs_a << " + " << rabs_b << " * alpha " << endl;
758 >    }
759 >  }
760 >  
761 >  float MPFResult;
762 >  float MPFResultError;
763 >  
764 >  float RabsResult;
765 >  float RabsResultError;
766 >
767 >  if(!gofast) {
768 >    MPF_FinalGraph->SetMarkerStyle(21);
769 >    MPF_FinalGraph->SetMarkerColor(kBlue);
770 >    MPF_FinalGraph->Fit("pol0",""); // quiet, use minos
771 >    
772 >    TF1 *mpf_pol0=(TF1*)MPF_FinalGraph->GetFunction("pol0");
773 >    TF1 *mpf_pol0_complete = new TF1("mpf_pol0_complete","[0]+[1]*x");
774 >    mpf_pol0_complete->SetLineWidth(1);
775 >    mpf_pol0_complete->SetLineColor(kBlue);
776 >    mpf_pol0->SetLineColor(kBlue);
777 >    MPF_FinalGraph->Draw("AP");
778 >    MPF_FinalGraph->GetXaxis()->SetTitle("p_{T}^{Z}");
779 >    MPF_FinalGraph->GetXaxis()->CenterTitle();
780 >    MPF_FinalGraph->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
781 >    MPF_FinalGraph->GetYaxis()->CenterTitle();
782 >    MPF_FinalGraph->Draw("AP");
783 >    mpf_pol0_complete->Draw("same");
784 >    
785 >    stringstream mpf_result;
786 >    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)}";
787 >    TText *rmpf_final_mark = write_title(mpf_result.str());
788 >    rmpf_final_mark->SetX(0.75);
789 >    rmpf_final_mark->SetY(0.75);
790 >    rmpf_final_mark->SetTextSize(0.03);
791 >    rmpf_final_mark->Draw();
792 >    DrawPrelim();
793 >    MPFResult=1/mpf_pol0->GetParameter(0);
794 >    MPFResultError=mpf_pol0->GetParError(0)/(mpf_pol0->GetParameter(0)*mpf_pol0->GetParameter(0));
795      
796 <    rabs_gr->SetMarkerStyle(21);
797 <    rabs_gr->SetMarkerColor(kRed);
798 <    rabs_gr->Draw("AP");
799 <    rabs_gr->Fit("pol1","");
800 <    rabs_gr->GetXaxis()->SetTitle("#alpha");
801 <    rabs_gr->GetXaxis()->CenterTitle();
802 <    rabs_gr->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
803 <    rabs_gr->GetYaxis()->CenterTitle();
804 <    TF1 *rabs_pol=(TF1*)rabs_gr->GetFunction("pol1");
805 <    rabs_pol->SetLineWidth(0);
806 <    TF1 *rabs_pol_complete = new TF1("rabs_pol_complete","[0]+[1]*x");
807 <    rabs_pol_complete->SetParameters(rabs_pol->GetParameters());
808 <    rabs_pol_complete->SetLineColor(kRed);
809 <    rabs_pol_complete->SetLineWidth(1);
810 <    FindMinMax(rabs_gr,min,max);
811 <    rabs_gr->GetYaxis()->SetRangeUser(min,max);
812 <    histo->GetYaxis()->SetRangeUser(min,max);
813 <    TPolyLine *rabs_fit_uncert = GetFitUncertaintyShape(rabs_gr, "pol1", min, max,0.0,0.4);
814 <    rabs_pol->SetLineColor(TColor::GetColor("#FA5858"));
815 <    rabs_pol->SetFillColor(TColor::GetColor("#FA5858"));
816 <    rabs_fit_uncert->SetLineColor(TColor::GetColor("#FA5858"));
817 <    rabs_fit_uncert->SetFillColor(TColor::GetColor("#FA5858"));
818 <    histo->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
521 <    histo->Draw();
522 <    rabs_fit_uncert->Draw("F");
523 <    rabs_fit_uncert->Draw("");
524 <    rabs_gr->Draw("P");
525 <    rabs_pol_complete->Draw("same");
526 <
527 <    
528 <    float rabs_a=rabs_pol->GetParameter(0);
529 <    float rabs_b=rabs_pol->GetParameter(1);
530 <    RABS_FinalGraph->SetPoint(ipt,0.5*(ptcuts[ipt]+ptcuts[ipt+1]),rabs_a);
531 <    RABS_FinalGraph->SetPointError(ipt,0,rabs_pol->GetParError(0));
532 <    
533 <    cout << "!+!+!+!+!+!+!+!+!+ Just added a point to the final plot : " << rabs_a << " +/- " << rabs_pol->GetParError(0) << endl;
534 <    
535 <    stringstream rabs_info;
536 <    rabs_info << "#splitline{#splitline{#splitline{" << ptcuts[ipt] << " GeV < p_{T}^{Z} < " << ptcuts[ipt+1] << " GeV }{ (R_{abs}) }}{Extrapolated value: " << std::setprecision(4) << rabs_a << "}}{#chi^{2}/NDF : " << rabs_pol->GetChisquare() << " / " << rabs_pol->GetNDF() << "}";
537 <    TText *rabs_mark = write_title(rabs_info.str());
538 <    rabs_mark->SetX(0.75);
539 <    rabs_mark->SetTextSize(0.03);
540 <    rabs_mark->SetY(0.75);
541 <    rabs_mark->Draw();
796 >    stringstream filename2;
797 >    filename2 << "Extrapolation/FINAL_RESULT_MPF";
798 >    CompleteSave(can,ContainerName+"/"+filename2.str());
799 >    
800 >    RABS_FinalGraph->SetMarkerStyle(21);
801 >    RABS_FinalGraph->SetMarkerStyle(21);
802 >    RABS_FinalGraph->SetMarkerColor(kRed);
803 >    RABS_FinalGraph->Fit("pol0","QE"); // quiet, use minos
804 >    RABS_FinalGraph->Draw("AP");
805 >    RABS_FinalGraph->GetXaxis()->SetTitle("p_{T}^{Z}");
806 >    RABS_FinalGraph->GetXaxis()->CenterTitle();
807 >    RABS_FinalGraph->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
808 >    RABS_FinalGraph->GetYaxis()->CenterTitle();
809 >    RABS_FinalGraph->Draw("AP");
810 >    TF1 *rabs_pol0=(TF1*)RABS_FinalGraph->GetFunction("pol0");
811 >    stringstream rabs_result;
812 >    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})}";;
813 >    TText *rabs_final_mark = write_title(rabs_result.str());
814 >    rabs_final_mark->SetX(0.75);
815 >    rabs_final_mark->SetY(0.75);
816 >    rabs_final_mark->SetTextSize(0.03);
817 >    rabs_final_mark->Draw();
818 >    DrawPrelim();
819      
820 +    RabsResult=1/rabs_pol0->GetParameter(0);
821 +    RabsResultError=rabs_pol0->GetParError(0)/(rabs_pol0->GetParameter(0)*rabs_pol0->GetParameter(0));
822      
823 <    RABS_ExtraPolatedResults[ipt]=rabs_a;
824 <    filename << filenamebkp << "__RABS";
825 <    DrawPrelim();
547 <    CompleteSave(can,filename.str());
548 <    cout << "RABS : " << rabs_a << " + " << rabs_b << " * alpha " << endl;
823 >    filename2.str("");
824 >    filename2 << "Extrapolation/FINAL_RESULT_RABS";
825 >    CompleteSave(can,ContainerName+"/"+filename2.str());
826    }
827 <  
551 <  MPF_FinalGraph->SetMarkerStyle(21);
552 <  MPF_FinalGraph->SetMarkerColor(kBlue);
553 <  MPF_FinalGraph->Fit("pol0",""); // quiet, use minos
554 <  
555 <  TF1 *mpf_pol0=(TF1*)MPF_FinalGraph->GetFunction("pol0");
556 <  TF1 *mpf_pol0_complete = new TF1("mpf_pol0_complete","[0]+[1]*x");
557 <  mpf_pol0_complete->SetLineWidth(1);
558 <  mpf_pol0_complete->SetLineColor(kBlue);
559 <  mpf_pol0->SetLineColor(kBlue);
560 <  MPF_FinalGraph->Draw("AP");
561 <  MPF_FinalGraph->GetXaxis()->SetTitle("p_{T}^{Z}");
562 <  MPF_FinalGraph->GetXaxis()->CenterTitle();
563 <  MPF_FinalGraph->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
564 <  MPF_FinalGraph->GetYaxis()->CenterTitle();
565 <  MPF_FinalGraph->Draw("AP");
566 <  mpf_pol0_complete->Draw("same");
567 <  
568 <  stringstream mpf_result;
569 <  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)}";
570 <  TText *rmpf_final_mark = write_title(mpf_result.str());
571 <  rmpf_final_mark->SetX(0.75);
572 <  rmpf_final_mark->SetY(0.75);
573 <  rmpf_final_mark->SetTextSize(0.03);
574 <  rmpf_final_mark->Draw();
575 <  DrawPrelim();
576 <  float MPFResult=1/mpf_pol0->GetParameter(0);
577 <  float MPFResultError=mpf_pol0->GetParError(0)/(mpf_pol0->GetParameter(0)*mpf_pol0->GetParameter(0));
578 <  
579 <  stringstream filename2;
580 <  filename2 << "Extrapolation/FINAL_RESULT_MPF";
581 <  CompleteSave(can,filename2.str());
582 <  
583 <  
584 <  RABS_FinalGraph->SetMarkerStyle(21);
585 <  RABS_FinalGraph->SetMarkerStyle(21);
586 <  RABS_FinalGraph->SetMarkerColor(kRed);
587 <  RABS_FinalGraph->Fit("pol0","QE"); // quiet, use minos
588 <  RABS_FinalGraph->Draw("AP");
589 <  RABS_FinalGraph->GetXaxis()->SetTitle("p_{T}^{Z}");
590 <  RABS_FinalGraph->GetXaxis()->CenterTitle();
591 <  RABS_FinalGraph->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
592 <  RABS_FinalGraph->GetYaxis()->CenterTitle();
593 <  RABS_FinalGraph->Draw("AP");
594 <  TF1 *rabs_pol0=(TF1*)RABS_FinalGraph->GetFunction("pol0");
595 <  stringstream rabs_result;
596 <  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})}";;
597 <  TText *rabs_final_mark = write_title(rabs_result.str());
598 <  rabs_final_mark->SetX(0.75);
599 <  rabs_final_mark->SetY(0.75);
600 <  rabs_final_mark->SetTextSize(0.03);
601 <  rabs_final_mark->Draw();
602 <  DrawPrelim();
603 <  
604 <  float RabsResult=1/rabs_pol0->GetParameter(0);
605 <  float RabsResultError=rabs_pol0->GetParError(0)/(rabs_pol0->GetParameter(0)*rabs_pol0->GetParameter(0));
606 <  
607 <  filename2.str("");
608 <  filename2 << "Extrapolation/FINAL_RESULT_RABS";
609 <  CompleteSave(can,filename2.str());
610 <  
827 >
828    can->cd();
829    MPF_FaceValueAtPoint3->SetMarkerStyle(21);
830    MPF_FaceValueAtPoint3->SetMarkerColor(kBlue);
# Line 615 | Line 832 | void new_data_mc_agreement_2d() {
832    MPF_FaceValueAtPoint3->Draw("AP");
833    MPF_FaceValueAtPoint3->GetXaxis()->SetTitle("p_{T}^{Z}");
834    MPF_FaceValueAtPoint3->GetXaxis()->CenterTitle();
835 <  MPF_FaceValueAtPoint3->GetYaxis()->SetTitle("<MPF>__{data}/<MPF>_{mc} for #alpha<0.3");
835 >  MPF_FaceValueAtPoint3->GetYaxis()->SetTitle(((string)"<MPF>__{data}/<MPF>_{mc} for #alpha<"+any2string(FaceValueToBeTakenAt)).c_str());
836    MPF_FaceValueAtPoint3->GetYaxis()->CenterTitle();
837    MPF_FaceValueAtPoint3->Draw("AP");
838    TF1 *faceval_pol0=(TF1*)MPF_FaceValueAtPoint3->GetFunction("pol0");
# Line 625 | Line 842 | void new_data_mc_agreement_2d() {
842    TF1 *faceval_pol0_complete = new TF1("faceval_pol0_complete","[0]+[1]*x");
843    faceval_pol0_complete->SetParameters(faceval_pol0->GetParameters());
844    stringstream faceval_result;
845 <  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)}";;
845 >  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<" << FaceValueToBeTakenAt << ")}";;
846    TText *faceval_final_mark = write_title(faceval_result.str());
847    faceval_final_mark->SetX(0.75);
848    faceval_final_mark->SetY(0.75);
# Line 636 | Line 853 | void new_data_mc_agreement_2d() {
853    float FaceValResult=1/faceval_pol0->GetParameter(0);
854    float FaceValResultError=faceval_pol0->GetParError(0)/(faceval_pol0->GetParameter(0)*faceval_pol0->GetParameter(0));
855    
856 +  stringstream filename2;
857    filename2.str("");
858    filename2 << "Extrapolation/FINAL_RESULT_MPF_FaceValp3";
859 <  CompleteSave(can,filename2.str());
859 >  CompleteSave(can,ContainerName+"/"+filename2.str());
860    
861    can->cd();
862    RABS_FaceValueAtPoint3->SetMarkerStyle(21);
# Line 647 | Line 865 | void new_data_mc_agreement_2d() {
865    RABS_FaceValueAtPoint3->Draw("AP");
866    RABS_FaceValueAtPoint3->GetXaxis()->SetTitle("p_{T}^{Z}");
867    RABS_FaceValueAtPoint3->GetXaxis()->CenterTitle();
868 <  RABS_FaceValueAtPoint3->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc} for #alpha<0.3");
868 >  RABS_FaceValueAtPoint3->GetYaxis()->SetTitle(((string)"<R_{abs}>_{data}/<R_{abs}>_{mc} for #alpha<"+any2string(FaceValueToBeTakenAt)).c_str());
869    RABS_FaceValueAtPoint3->GetYaxis()->CenterTitle();
870    RABS_FaceValueAtPoint3->Draw("AP");
871    TF1 *faceval_pol0RABS=(TF1*)RABS_FaceValueAtPoint3->GetFunction("pol0");
# Line 657 | Line 875 | void new_data_mc_agreement_2d() {
875    faceval_pol0RABS->SetLineColor(kRed);
876    faceval_pol0RABS_complete->Draw("same");
877    stringstream RABSfaceval_result;
878 <  RABSfaceval_result << "#splitline{C_{abs}= " << std::setprecision(5) << 1/faceval_pol0RABS->GetParameter(0) << " #pm " << std::setprecision(2) << faceval_pol0RABS->GetParError(0)/(faceval_pol0RABS->GetParameter(0)*faceval_pol0RABS->GetParameter(0)) << "}{ (R_{abs} at #alpha<0.3)}";;
878 >  RABSfaceval_result << "#splitline{C_{abs}= " << std::setprecision(5) << 1/faceval_pol0RABS->GetParameter(0) << " #pm " << std::setprecision(2) << faceval_pol0RABS->GetParError(0)/(faceval_pol0RABS->GetParameter(0)*faceval_pol0RABS->GetParameter(0)) << "}{ (R_{abs} at #alpha<" << FaceValueToBeTakenAt << ")}";;
879    TText *RABSfaceval_final_mark = write_title(RABSfaceval_result.str());
880    RABSfaceval_final_mark->SetX(0.75);
881    RABSfaceval_final_mark->SetY(0.75);
# Line 670 | Line 888 | void new_data_mc_agreement_2d() {
888    
889    filename2.str("");
890    filename2 << "Extrapolation/FINAL_RESULT_RABS_FaceValp3";
891 <  CompleteSave(can,filename2.str());
891 >  CompleteSave(can,ContainerName+"/"+filename2.str());
892    
893    cout << "FINAL RESULTS: " << endl;
894 <  cout << "MPF : " << MPFResult  << " +/- " << MPFResultError  << endl;
895 <  cout << "Rabs: " << RabsResult << " +/- " << RabsResultError << endl;
896 <  cout << "Face value at 0.3 :" << FaceValResult << " +/- " << FaceValResultError << "  (MPF) " << endl;
897 <  cout << "Face value at 0.3 :" << RABSFaceValResult << " +/- " << RABSFaceValResultError << "  (RABS) " << endl;
894 >  if(!gofast) cout << "MPF : " << MPFResult  << " +/- " << MPFResultError  << endl;
895 >  if(!gofast) cout << "Rabs: " << RabsResult << " +/- " << RabsResultError << endl;
896 >  cout << "Face value at " << FaceValueToBeTakenAt << " :" << FaceValResult << " +/- " << FaceValResultError << "  (MPF) " << endl;
897 >  cout << "Face value at " << FaceValueToBeTakenAt << " :" << RABSFaceValResult << " +/- " << RABSFaceValResultError << "  (RABS) " << endl;
898    
899    delete can;
900    
901 <  
901 >  results.MPF_Result_FaceValue=Value(FaceValResult,FaceValResultError);
902 >  results.Rabs_Result_FaceValue=Value(RABSFaceValResult,RABSFaceValResultError);
903 >  if(!gofast) results.MPF_Result_Extrapolated=Value(MPFResult,MPFResultError);
904 >  if(!gofast) results.Rabs_Result_Extrapolated=Value(RabsResult,RabsResultError);
905 >  cout << results << endl;
906 >  return results;
907   }
908  
909  
# Line 703 | Line 926 | void DoPUStudy(string identifier) {
926        cout << numVtxcuts[i-1] << ";";
927      }
928      
929 +    sSpecialCut << " && " << identifier << "_alpha<0.3";
930      TCut SpecialCut(sSpecialCut.str().c_str());
931      
932      
933 <    malpha[i] = allsamples.Draw("malpha"+any2string(i),"mpf",1000,0,10,"MPF","events",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&SpecialCut,mc,luminosity);
934 <    dalpha[i] = allsamples.Draw("dalpha"+any2string(i),"mpf",1000,0,10,"MPF","events",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&SpecialCut,data,luminosity);
933 >    malpha[i] = allsamples.Draw("malpha"+any2string(i),"mpf",1000,0,2,"MPF","events",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&SpecialCut,mc,luminosity);
934 >    dalpha[i] = allsamples.Draw("dalpha"+any2string(i),"mpf",1000,0,2,"MPF","events",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&SpecialCut,data,luminosity);
935      
936      
937      float a=dalpha[i]->GetMean();
# Line 718 | Line 942 | void DoPUStudy(string identifier) {
942      float error = TMath::Sqrt(  (1/(b*b)) * (da*da) + ((a*a)/(b*b))*db*db);
943    
944      cout << malpha[i]->GetMean() << ";" << malpha[i]->GetMeanError() << ";" << dalpha[i]->GetMean() << ";" << dalpha[i]->GetMeanError() << ";" << malpha[i]->Integral()<<";"<<dalpha[i]->Integral() << ";"<<factor << ";" << error <<  endl;
721    malpha[i]->SetLineColor(i+1);
722    if(i==0) malpha[i]->Draw("histo");
723    else malpha[i]->Draw("histo,same");
724  }
725  
726  CompleteSave(can,"AlphaOverview_"+identifier);
727  delete can;
728  
729 /*  for(int i=1;i<8;i++) {
945      delete malpha[i];
946      delete dalpha[i];
947    }
948 <  */
948 >  
949 >  delete can;
950    
951   }
952  
953   void ScenarioComparison() {
954   //very dumb way to do this.
955 +
956    TGraphAsymmErrors *gr[5];
957    TF1 *fit[5];
958 <  bool dofit=false;
958 >  bool dofit=true;
959    
960    TCanvas *can = new TCanvas("can","can",500,500);
961    
# Line 746 | Line 963 | void ScenarioComparison() {
963    leg->SetBorderSize(0);
964    leg->SetFillColor(kWhite);
965    
966 <  for(int i=0;i<5;i++) {
966 >  for(int i=0;i<6;i++) {
967      gr[i] = new TGraphAsymmErrors();
968 <    double x[4]   = {5,12.5,17.5,20};
969 <    double dxl[4] = {5,2.5,2.5,2.5};
968 >    double x[4]   = {5,12.5,17.5,25};
969 >    double dxl[4] = {5,2.5,2.5,5};
970      double dxh[4] = {5,2.5,2.5,10};
971      
972      string name="";
973      string color="";
974      
975      if(i==0) {
976 <      double y[4]   = {1.008,1.013,1.003,1.007};
977 <      double dyl[4] = {0.011,0.009,0.013,0.020};
978 <      double dyh[4] = {0.011,0.009,0.013,0.020};
976 >      double y[4]   = {1.01935,1.01558,0.996438,1.01901};
977 >      double dyl[4] = {0.0132973,0.0131661,0.017629,0.0319767};
978 >      double dyh[4] = {0.0132973,0.0131661,0.017629,0.0319767};
979        name="30/30";
980        gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
981        color="#a01d99";
982      }
983      if(i==1) {
984 <      double y[4]   = {1.008+0.001,1.013+0.001,1.003+0.001,1.007+0.001};
985 <      double dyl[4] = {0.011,0.009,0.013,0.020};
986 <      double dyh[4] = {0.011,0.009,0.013,0.020};
984 >      double y[4]   = {1.00577,1.00706,0.991029,0.954098};
985 >      double dyl[4] = {0.0135575,0.0136548,0.0170778,0.0263919};
986 >      double dyh[4] = {0.0135575,0.0136548,0.0170778,0.0263919};
987        gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
988        name="30/15";
989        color="#2464bc";
990      }
991      if(i==2) {
992 <      double y[4]   = {1.014,1.026,1.014,1.059};
993 <      double dyl[4] = {0.011,0.010,0.014,0.023};
994 <      double dyh[4] = {0.011,0.010,0.014,0.023};
992 >      double y[4]   = {1.023,1.00648,1.02217,1.03181};
993 >      double dyl[4] = {0.0157372,0.0187778,0.0261512,0.0422295};
994 >      double dyh[4] = {0.0157372,0.0187778,0.0261512,0.0422295};
995        gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
996        name="15/15";
997        color="#f49230";
998      }
999      if(i==3) {
1000 <      double y[4]   = {1.015,1.020,1.008,1.050};
1001 <      double dyl[4] = {0.011,0.010,0.014,0.022};
1002 <      double dyh[4] = {0.011,0.010,0.014,0.022};
1000 >      double y[4]   = {1.02029,1.01338,0.999769,1.00357};
1001 >      double dyl[4] = {0.0129509,0.0120127,0.0160537,0.0268972};
1002 >      double dyh[4] = {0.0129509,0.0120127,0.0160537,0.0268972};
1003        gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
1004 <      name="20/15";
1004 >      name="CHS 30/30";
1005        color="#c72f3c";
1006      }
1007      if(i==4) {
1008 <      double y[4]   = {1.015+0.001,1.020+0.001,1.008+0.001,1.050+0.001};
1009 <      double dyl[4] = {0.011,0.010,0.014,0.022};
1010 <      double dyh[4] = {0.011,0.010,0.014,0.022};
1008 >      double y[4]   = {1.00274,1.01056,0.993214,0.997963};
1009 >      double dyl[4] = {0.0127179,0.0140344,0.0162814,0.0278914};
1010 >      double dyh[4] = {0.0127179,0.0140344,0.0162814,0.0278914};
1011        gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
1012 <      name="20/20";
1012 >      name="CHS 30/15";
1013        color="#7cb433";
1014      }
1015 <    
1016 <    
1015 >    if(i==5) {
1016 >      double y[4]   = {1.04103,1.03726,1.04649,1.05383};
1017 >      double dyl[4] = {0.0156465,0.0168194,0.0231982,0.0376114};
1018 >      double dyh[4] = {0.0156465,0.0168194,0.0231982,0.0376114};
1019 >      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
1020 >      name="CHS 15/15";
1021 >      color="#DF013A";
1022 >    }
1023 >    
1024 >    
1025 >    gStyle->SetOptFit(0);
1026      gr[i]->SetTitle();
1027      gr[i]->GetXaxis()->SetTitle("N(Vtx)");
1028      gr[i]->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
# Line 832 | Line 1058 | void ScenarioComparisonInclusive() {
1058    
1059    TCanvas *can = new TCanvas("can","can",500,500);
1060    
1061 <  TLegend *leg = new TLegend(0.2,0.2,0.5,0.45);
1061 >  TLegend *leg = new TLegend(0.2,0.2,0.65,0.45);
1062    leg->SetBorderSize(0);
1063    leg->SetFillColor(kWhite);
1064    
1065 <  for(int i=0;i<5;i++) {
1065 >  for(int i=0;i<6;i++) {
1066      gr[i] = new TGraphAsymmErrors();
1067 <    double x[4]   = {5,12.5,17.5,20};
1068 <    double dxl[4] = {5,2.5,2.5,2.5};
1069 <    double dxh[4] = {5,2.5,2.5,10};
1067 >    double x[4]   = {5,12.5,17.5,25};
1068 >    double dxl[4] = {5,2.5,2.5,5};
1069 >    double dxh[4] = {5,2.5,2.5,5};
1070      
1071      string name="";
1072      string color="";
1073      
1074      if(i==0) {
1075 <      double y[4]   = {1.01103,1.01496,1.01486,1.02785};
1076 <      double dyl[4] = {0.00235044,0.00225264,0.00313574,0.00541047};
1077 <      double dyh[4] = {0.00235044,0.00225264,0.00313574,0.00541047};
1075 >      double y[4]   = {1.01093,1.01379,1.0147,1.02766};
1076 >      double dyl[4] = {0.00243596,0.00242534,0.00341941,0.00600663};
1077 >      double dyh[4] = {0.00243596,0.00242534,0.00341941,0.00600663};
1078        name="30/30";
1079        gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
1080        color="#a01d99";
1081      }
1082      if(i==1) {
1083 <      double y[4]   = {1.01103+0.001,1.01496+0.001,1.01486+0.001,1.02785+0.001};
1084 <      double dyl[4] = {0.00235044,0.00225264,0.00313574,0.00541047};
1085 <      double dyh[4] = {0.00235044,0.00225264,0.00313574,0.00541047};
1083 >      double y[4]   = {1.0068,1.00779,0.999692,0.994484};
1084 >      double dyl[4] = {0.00249233,0.00257169,0.00371399,0.00631138};
1085 >      double dyh[4] = {0.00249233,0.00257169,0.00371399,0.00631138};
1086        gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
1087        name="30/15";
1088        color="#2464bc";
1089      }
1090      if(i==2) {
1091 <      double y[4]   = {1.01575,1.02879,1.0325,1.05456};
1092 <      double dyl[4] = {0.00288281,0.00278536,0.00367451,0.00590421};
1093 <      double dyh[4] = {0.00288281,0.00278536,0.00367451,0.00590421};
1091 >      double y[4]   = {1.01208,1.01991,1.02185,1.03479};
1092 >      double dyl[4] = {0.00366745,0.00400001,0.00584054,0.0104661};
1093 >      double dyh[4] = {0.00366745,0.00400001,0.00584054,0.0104661};
1094        gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
1095        name="15/15";
1096        color="#f49230";
1097      }
1098      if(i==3) {
1099 <      double y[4]   = {1.01221,1.01922,1.02396,1.04677};
1100 <      double dyl[4] = {0.00263821,0.00260853,0.00353673,0.00584728};
1101 <      double dyh[4] = {0.00263821,0.00260853,0.00353673,0.00584728};
1099 >      double y[4]   = {1.01177,1.01568,1.01616,1.03211};
1100 >      double dyl[4] = {0.00217317,0.00203014,0.00276838,0.00465643};
1101 >      double dyh[4] = {0.00217317,0.00203014,0.00276838,0.00465643};
1102        gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
1103 <      name="20/15";
1103 >      name="CHS 30/30";
1104        color="#c72f3c";
1105      }
1106      if(i==4) {
1107 <      double y[4]   = {1.01221+0.001,1.01922+0.001,1.02396+0.001,1.04677+0.001};
1108 <      double dyl[4] = {0.00263821,0.00260853,0.00353673,0.00584728};
1109 <      double dyh[4] = {0.00263821,0.00260853,0.00353673,0.00584728};
1107 >      double y[4]   = {1.00751,1.00964,1.00804,1.01073};
1108 >      double dyl[4] = {0.00224851,0.00215563,0.0029592,0.00499805};
1109 >      double dyh[4] = {0.00224851,0.00215563,0.0029592,0.00499805};
1110        gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
1111 <      name="20/20";
1111 >      name="CHS 30/15";
1112        color="#7cb433";
1113      }
1114 +    if(i==5) {
1115 +      double y[4]   = {1.01557,1.0318,1.03535,1.06331};
1116 +      double dyl[4] = {0.00300783,0.00302486,0.00421064,0.00739416};
1117 +      double dyh[4] = {0.00300783,0.00302486,0.00421064,0.00739416};
1118 +      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
1119 +      name="CHS 15/15";
1120 +      color="#DF013A";
1121 +    }
1122      
1123 <    
1123 >    gStyle->SetOptFit(0);
1124      gr[i]->SetTitle();
1125      gr[i]->GetXaxis()->SetTitle("N(Vtx)");
1126      gr[i]->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
# Line 902 | Line 1136 | void ScenarioComparisonInclusive() {
1136        fit[i] = (TF1*)gr[i]->GetFunction("pol1");
1137        fit[i]->SetLineColor(TColor::GetColor(color.c_str()));
1138      }
1139 <    leg->AddEntry(gr[i],name.c_str(),"l");
1139 >    
1140 >    stringstream nameplus;
1141 >    nameplus << name << " (slope: " << std::setprecision(3) << fit[i]->GetParameter(1) << ")";
1142 >    
1143 >    leg->AddEntry(gr[i],nameplus.str().c_str(),"l");
1144      if(i==0) gr[i]->Draw("AP*");
1145      else gr[i]->Draw("P");
1146    }
# Line 918 | Line 1156 | void compare_selection(string identifier
1156    bool recognized_scenario=false;
1157    
1158    cout << "Running with identifier " << identifier << endl;
1159 +  if(identifier=="Zb30") {
1160 +      LeadingB=TCut ("Zb30_bTagProbCSVBP[0]>0.898");
1161 +      EtaB=TCut("abs(Zb30_pfJetGoodEta[0])<1.3");
1162 +      PhiZcut=TCut("abs(Zb30_pfJetDphiZ[0])>2.7");
1163 +      recognized_scenario=true;
1164 +  }
1165    if(identifier=="Zb3010") {
1166        LeadingB=TCut ("Zb3010_bTagProbCSVBP[0]>0.898");
1167        EtaB=TCut("abs(Zb3010_pfJetGoodEta[0])<1.3");
# Line 925 | Line 1169 | void compare_selection(string identifier
1169        recognized_scenario=true;
1170    }
1171    
1172 <  if(identifier=="Zb3010") {
1173 <      LeadingB=TCut ("Zb3010_bTagProbCSVBP[0]>0.898");
1174 <      EtaB=TCut("abs(Zb3010_pfJetGoodEta[0])<1.3");
1175 <      PhiZcut=TCut("abs(Zb3010_pfJetDphiZ[0])>2.7");
1172 >  if(identifier=="ZbCHS3010") {
1173 >      LeadingB=TCut ("ZbCHS1510_bTagProbCSVBP[0]>0.898");
1174 >      EtaB=TCut("abs(ZbCHS1510_pfJetGoodEta[0])<1.3");
1175 >      PhiZcut=TCut("abs(ZbCHS1510_pfJetDphiZ[0])>2.7");
1176        recognized_scenario=true;
1177    }
1178    
# Line 939 | Line 1183 | void compare_selection(string identifier
1183        recognized_scenario=true;
1184    }
1185    
1186 <  if(identifier=="Zb301030") {
1187 <      LeadingB=TCut ("Zb301030_bTagProbCSVBP[0]>0.898");
1188 <      EtaB=TCut("abs(Zb301030_pfJetGoodEta[0])<1.3");
1189 <      PhiZcut=TCut("abs(Zb301030_pfJetDphiZ[0])>2.7");
1186 >  if(identifier=="ZbCHS301030") {
1187 >      LeadingB=TCut ("ZbCHS301030_bTagProbCSVBP[0]>0.898");
1188 >      EtaB=TCut("abs(ZbCHS301030_pfJetGoodEta[0])<1.3");
1189 >      PhiZcut=TCut("abs(ZbCHS301030_pfJetDphiZ[0])>2.7");
1190        recognized_scenario=true;
1191    }
1192    
# Line 953 | Line 1197 | void compare_selection(string identifier
1197        recognized_scenario=true;
1198    }
1199    
1200 <  if(identifier=="Zb301010") {
1201 <      LeadingB=TCut ("Zb301010_bTagProbCSVBP[0]>0.898");
1202 <      EtaB=TCut("abs(Zb301010_pfJetGoodEta[0])<1.3");
1203 <      PhiZcut=TCut("abs(Zb301010_pfJetDphiZ[0])>2.7");
1200 >  if(identifier=="ZbCHS301010") {
1201 >      LeadingB=TCut ("ZbCHS301010_bTagProbCSVBP[0]>0.898");
1202 >      EtaB=TCut("abs(ZbCHS301010_pfJetGoodEta[0])<1.3");
1203 >      PhiZcut=TCut("abs(ZbCHS301010_pfJetDphiZ[0])>2.7");
1204        recognized_scenario=true;
1205    }
1206    
# Line 967 | Line 1211 | void compare_selection(string identifier
1211        recognized_scenario=true;
1212    }
1213    
1214 <  if(identifier=="Zb301010") {
1215 <      LeadingB=TCut ("Zb301010_bTagProbCSVBP[0]>0.898");
1216 <      EtaB=TCut("abs(Zb301010_pfJetGoodEta[0])<1.3");
1217 <      PhiZcut=TCut("abs(Zb301010_pfJetDphiZ[0])>2.7");
1214 >  if(identifier=="ZbCHS301010") {
1215 >      LeadingB=TCut ("ZbCHS301010_bTagProbCSVBP[0]>0.898");
1216 >      EtaB=TCut("abs(ZbCHS301010_pfJetGoodEta[0])<1.3");
1217 >      PhiZcut=TCut("abs(ZbCHS301010_pfJetDphiZ[0])>2.7");
1218        recognized_scenario=true;
1219    }
1220    
# Line 981 | Line 1225 | void compare_selection(string identifier
1225        recognized_scenario=true;
1226    }
1227    
1228 <  if(identifier=="Zb3010_SecEta3") {
1229 <      LeadingB=TCut ("Zb3010_SecEta3_bTagProbCSVBP[0]>0.898");
1230 <      EtaB=TCut("abs(Zb3010_SecEta3_pfJetGoodEta[0])<1.3");
1231 <      PhiZcut=TCut("abs(Zb3010_SecEta3_pfJetDphiZ[0])>2.7");
1228 >  if(identifier=="ZbCHS3010_SecEta3") {
1229 >      LeadingB=TCut ("ZbCHS3010_SecEta3_bTagProbCSVBP[0]>0.898");
1230 >      EtaB=TCut("abs(ZbCHS3010_SecEta3_pfJetGoodEta[0])<1.3");
1231 >      PhiZcut=TCut("abs(ZbCHS3010_SecEta3_pfJetDphiZ[0])>2.7");
1232        recognized_scenario=true;
1233    }
1234    
1235 <  if(identifier=="Zb3010_SecEta5") {
1236 <      LeadingB=TCut ("Zb3010_SecEta5_bTagProbCSVBP[0]>0.898");
1237 <      EtaB=TCut("abs(Zb3010_SecEta5_pfJetGoodEta[0])<1.3");
1238 <      PhiZcut=TCut("abs(Zb3010_SecEta5_pfJetDphiZ[0])>2.7");
1235 >  if(identifier=="ZbCHS3010_SecEta5") {
1236 >      LeadingB=TCut ("ZbCHS3010_SecEta5_bTagProbCSVBP[0]>0.898");
1237 >      EtaB=TCut("abs(ZbCHS3010_SecEta5_pfJetGoodEta[0])<1.3");
1238 >      PhiZcut=TCut("abs(ZbCHS3010_SecEta5_pfJetDphiZ[0])>2.7");
1239        recognized_scenario=true;
1240    }
1241    
1242 <  if(identifier=="Zb3010_p5Clean") {
1243 <      LeadingB=TCut ("Zb3010_p5Clean_bTagProbCSVBP[0]>0.898 && Zb3010_p5Clean_IsClean");
1244 <      EtaB=TCut("abs(Zb3010_p5Clean_pfJetGoodEta[0])<1.3 && Zb3010_p5Clean_IsClean");
1245 <      PhiZcut=TCut("abs(Zb3010_p5Clean_pfJetDphiZ[0])>2.7 && Zb3010_p5Clean_IsClean");
1242 >  if(identifier=="ZbCHS3010_p5Clean") {
1243 >      LeadingB=TCut ("ZbCHS3010_p5Clean_bTagProbCSVBP[0]>0.898 && ZbCHS3010_p5Clean_IsClean");
1244 >      EtaB=TCut("abs(ZbCHS3010_p5Clean_pfJetGoodEta[0])<1.3 && ZbCHS3010_p5Clean_IsClean");
1245 >      PhiZcut=TCut("abs(ZbCHS3010_p5Clean_pfJetDphiZ[0])>2.7 && ZbCHS3010_p5Clean_IsClean");
1246        recognized_scenario=true;
1247    }
1248    
1249 <  if(identifier=="Zb3010CHS") {
1250 <      LeadingB=TCut ("Zb3010CHS_bTagProbCSVBP[0]>0.898");
1251 <      EtaB=TCut("abs(Zb3010CHS_pfJetGoodEta[0])<1.3");
1252 <      PhiZcut=TCut("abs(Zb3010CHS_pfJetDphiZ[0])>2.7");
1249 >  if(identifier=="ZbCHS3010") {
1250 >      LeadingB=TCut ("ZbCHS3010_bTagProbCSVBP[0]>0.898");
1251 >      EtaB=TCut("abs(ZbCHS3010_pfJetGoodEta[0])<1.3");
1252 >      PhiZcut=TCut("abs(ZbCHS3010_pfJetDphiZ[0])>2.7");
1253 >      recognized_scenario=true;
1254 >  }
1255 >  if(identifier=="ZbCHS3010") {
1256 >      LeadingB=TCut ("ZbCHS3010_bTagProbCSVBP[0]>0.898");
1257 >      EtaB=TCut("abs(ZbCHS3010_pfJetGoodEta[0])<1.3");
1258 >      PhiZcut=TCut("abs(ZbCHS3010_pfJetDphiZ[0])>2.7");
1259 >      recognized_scenario=true;
1260 >  }
1261 >  if(identifier=="ZbCHS30") {
1262 >      LeadingB=TCut ("ZbCHS30_bTagProbCSVBP[0]>0.898");
1263 >      EtaB=TCut("abs(ZbCHS30_pfJetGoodEta[0])<1.3");
1264 >      PhiZcut=TCut("abs(ZbCHS30_pfJetDphiZ[0])>2.7");
1265        recognized_scenario=true;
1266    }
1267    
# Line 1017 | Line 1273 | void compare_selection(string identifier
1273    cout << "   " << (const char*) PhiZcut << endl;
1274    
1275    cout << endl << endl;
1276 +
1277 + //  ScenarioComparison();
1278 + //  ScenarioComparisonInclusive();
1279 +
1280 + //    DoPUStudy(identifier);
1281 + //    SpecialCutFlow(identifier);
1282    
1283 < //    ScenarioComparison();
1284 < //    ScenarioComparisonInclusive();
1285 <  
1286 <  
1025 < //     DoPUStudy(identifier);
1026 <  
1027 < //   SpecialCutFlow(identifier);
1028 < //    draw_kin_variable("pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",17000,30,200,"Z p_{T}","Official/"+identifier+"/Zpt__"+identifier,1);
1029 < //   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);
1030 < //   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);
1031 < //   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);
1283 > //    draw_kin_variable("pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000",17000,30,200,"Z p_{T}","Official/"+identifier+"/Zpt__"+identifier,1);
1284 > //    draw_kin_variable(identifier+"_pfJetGoodPt[0]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000",40,0,200,"Leading Jet Pt (B)","Official/"+identifier+"/LeadingJetPt__"+identifier,1);
1285 > //    draw_kin_variable(identifier+"_pfJetGoodPt[1]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000",40,0,200,"Sub-Leading Jet Pt [GeV]","Official/"+identifier+"/Jet2Pt__"+identifier,1);
1286 > //    draw_kin_variable(identifier+"_pfJetGoodPt[1]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000",20,0,2,"p_{T}^{J2} / p_{T}^{Z}","Official/"+identifier+"/SubLeadingJetPt_Over_ZPt__"+identifier,1);
1287    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);
1288    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);
1289    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);
1290    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);
1291    draw_kin_variable(identifier+"_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut,20,0,2,"#alpha","Official/"+identifier+"/alpha_full__"+identifier,1);
1292 < //   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);
1292 >   draw_kin_variable(identifier+"_pfBJetDphiZ[0]",ZplusBsel&&LeadingB&&"pt>30&&pt<1000&&pfBJetDphiZ[0]>0",30,0,3.2,"#delta#phi (Z,b lead)","Official/"+identifier+"/DeltaPhiZBlead__"+identifier,0);
1293    
1294 < //   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");
1295 < //   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");
1296 <
1042 <  
1043 < //   */
1294 >   draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>30&&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, 30 GeV < p_{T}^{Z} < 1000 GeV");
1295 >   draw_kin_variable(identifier+"_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>30&&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, 30 GeV < p_{T}^{Z} < 1000 GeV");
1296 >   draw_kin_variable(identifier+"_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>30&&pt<50&&"+identifier+"_alpha<0.3").c_str()),20,0,2,"p_{T} b-jet / p_{T} Z","Official/"+identifier+"/ptb_over_ptz___alpha_smaller_0p3___pt3050___"+identifier,0,"#alpha<0.3, 30 GeV < p_{T}^{Z} < 1000 GeV}{30<p_{T}^{Z}<50 GeV}");
1297    
1298   }
1299  
1300   void GetNumberEventsInsideOutsideAlphaWindow(TCut specialcut, string title) {
1301    
1302 < //  TCut cut =  ZplusBsel&&LeadingB&&PhiZcut&&specialcut&&TCut("Zb3010_alpha<0.3");
1302 > //  TCut cut =  ZplusBsel&&LeadingB&&PhiZcut&&specialcut&&TCut("ZbCHS3010_alpha<0.3");
1303    TCut cut = specialcut;
1304    TCut in  =  EtaB;
1305 <  TCut out =  TCut("abs(Zb3010_pfJetGoodEta[0])>1.3");
1305 >  TCut out =  TCut("abs(ZbCHS3010_pfJetGoodEta[0])>1.3");
1306  
1307    TH1F *data_IN  = allsamples.Draw("data_IN",    "mll",1,0,150, "nothing [GeV]", "events", cut&&in, data,luminosity);
1308    TH1F *data_OUT = allsamples.Draw("data_OUT",   "mll",1,0,150, "nothing [GeV]", "events", cut&&out,data,luminosity);
# Line 1068 | Line 1321 | void GetNumberEventsInsideOutsideAlphaWi
1321    delete mc_OUT;
1322   }
1323  
1324 < void GetNumberEventsInsideOutsideAlphaWindow() {
1324 > void TEST_GetNumberEventsInsideOutsideAlphaWindow() {
1325    TCanvas *ca = new TCanvas("ca","ca");
1326   //   GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==0"), "ee");
1327   //   GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==1"), "mm");
# Line 1095 | Line 1348 | void GetNumberEventsInsideOutsideAlphaWi
1348    GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2"), "ee/mm");
1349    
1350    cout << "alpha cut" << endl;
1351 <  GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==0")&&TCut("Zb3010_alpha<0.3"), "ee");
1352 <  GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==1")&&TCut("Zb3010_alpha<0.3"), "mm");
1353 <  GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2")&&TCut("Zb3010_alpha<0.3"), "ee/mm");
1351 >  GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==0")&&TCut("ZbCHS3010_alpha<0.3"), "ee");
1352 >  GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==1")&&TCut("ZbCHS3010_alpha<0.3"), "mm");
1353 >  GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2")&&TCut("ZbCHS3010_alpha<0.3"), "ee/mm");
1354    
1355    delete ca;
1356    
# Line 1105 | Line 1358 | void GetNumberEventsInsideOutsideAlphaWi
1358   }
1359    
1360  
1361 <
1362 <    
1363 < void do_basic_ZB_analysis(int do_inclusive) {
1361 > void MetSpectrumDepletionIllustration(int isdata) {
1362 >  
1363 >  TCanvas *can = new TCanvas("can","can");
1364 >  TCut basecut(ZplusBsel&&LeadingB&&EtaB&&PhiZcut);
1365 >  
1366 >  cutWeight=STDWEIGHT*TCut("ZbCHS3010_BTagWgtTUp*(1.0)*(ZbCHS3010_bTagProbCSVBP[0]>0.898)");
1367 >  TH1F *standard  = allsamples.Draw("standard","met[1]",80,0,400, "MET [GeV]", "events", basecut,isdata,luminosity);
1368 >  
1369 >  cutWeight=STDWEIGHT*TCut("ZbCHS3010_BTagWgtTUp*(1.0*softMuon)*(ZbCHS3010_bTagProbCSVBP[0]>0.898)");
1370 >  TH1F *Neutrinos  = allsamples.Draw("Neutrinos","met[1]",80,0,400, "MET [GeV]", "events", basecut,isdata,luminosity);
1371 >  
1372 >  cutWeight=STDWEIGHT*TCut("ZbCHS3010_BTagWgtTUp*(1.0*(1.0-softMuon))*(ZbCHS3010_bTagProbCSVBP[0]>0.898)");
1373 >  TH1F *NoNeutrinos  = allsamples.Draw("NoNeutrinos","met[1]",80,0,400, "MET [GeV]", "events", basecut,isdata,luminosity);
1374 >  
1375 >  standard->Scale(1.0/standard->Integral());
1376 >  Neutrinos->Scale(1.0/Neutrinos->Integral());
1377 >  NoNeutrinos->Scale(1.0/NoNeutrinos->Integral());
1378 >  
1379 >  Neutrinos->SetLineColor(kGreen);
1380 >  NoNeutrinos->SetLineColor(kRed);
1381 >  
1382 >  can->SetLogy(1);
1383 >  TLegend *leg = make_legend();
1384 >  leg->AddEntry(standard,"Standard sample","ep");
1385 >  leg->AddEntry(Neutrinos,"#nu enriched sample","l");
1386 >  leg->AddEntry(NoNeutrinos,"#nu depleted sample","l");
1387 >  leg->SetX1(0.9*leg->GetX1());
1388 >  standard->Draw("e1");
1389 >  Neutrinos->Draw("histo,same");
1390 >  NoNeutrinos->Draw("histo,same");
1391 >  if(isdata==data) leg->SetHeader("Data:");
1392 >  else leg->SetHeader("MC:");
1393 >  leg->Draw();
1394 >  DrawPrelim();
1395 >  
1396 >  if(isdata==data) CompleteSave(can,"NeutrinoMetSpectrum__Data");
1397 >  else CompleteSave(can,"NeutrinoMetSpectrum__MC");
1398 >  delete standard;
1399 >  delete Neutrinos;
1400 >  delete NoNeutrinos;
1401 >  delete can;
1402 >  CleanLegends();
1403 > }
1404 >  
1405 >  
1406 >  
1407    
1112  //https://twiki.cern.ch/twiki/bin/viewauth/CMS/BtagPOG#2011_Data_and_MC
1113  //https://twiki.cern.ch/twiki/pub/CMS/BtagPOG/SFb-ttbar_payload.txt
1408    
1409 +  
1410 +    
1411 + void do_basic_ZB_analysis(bool DoCompleteAnalysis) {
1412 +  STDWEIGHT=TCut(cutWeight);
1413 +  cutWeight=TCut(STDWEIGHT*TCut("ZbCHS3010_BTagWgtT"));
1414    cout << "The lepton requirement is " << (const char*) leptoncut << endl;
1415 <  if(do_inclusive) {
1416 <    LeadingB=TCut("Zb3010_bTagProbCSVBP[0]>=-10000");
1417 <    cutWeight=TCut("(weight*(weight<1000)*(is_data+(!is_data)))");
1418 <    write_info(__FUNCTION__,"Inclusive selection - no btag requirement, no btag efficiency/mistag correction (weight)");
1419 <  } else {
1420 <    LeadingB=TCut("Zb3010_bTagProbCSVBP[0]>0.898");
1421 <    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))");
1422 <    write_info(__FUNCTION__,"Exclusive selection - switched on btag requirement and btag efficiency/mistag correction (weight)");
1423 <  }
1424 <  
1425 <  bool doquick=true;
1426 <  
1427 <  TCanvas *zbcanvas = new TCanvas("zbcanvas","zbcanvas");
1428 <  
1429 < //     compare_selection("Zb3010_p5Clean");
1430 < //     compare_selection("Zb3010CHS");
1431 < //     compare_selection("Zb3010");
1432 < //     compare_selection("Zb301030");
1433 < //     compare_selection("Zb1530");
1434 < //     compare_selection("Zb301010");
1435 < //     compare_selection("Zb3010_SecEta3");
1436 < //     compare_selection("Zb3010_SecEta5");
1437 < //     compare_selection("Zb1510");
1438 < //     compare_selection("Zb301010");
1439 < //     compare_selection("ZbMikko");
1141 < //     compare_selection("Zb3010");
1142 < //     compare_selection("Zb40");
1415 >  bool doquick=false;
1416 >  
1417 >   TCanvas *zbcanvas = new TCanvas("zbcanvas","zbcanvas");
1418 >  
1419 >  
1420 >   MetSpectrumDepletionIllustration(data);
1421 >   MetSpectrumDepletionIllustration(mc);
1422 >   assert(0);
1423 >  
1424 > //   compare_selection("ZbCHS3010_p5Clean");
1425 > //   compare_selection("ZbCHS3010");
1426 > //   compare_selection("Zb30");
1427 > //   compare_selection("Zb3010");
1428 > //   compare_selection("Zb1510");
1429 > //   compare_selection("ZbCHS30");
1430 > //   compare_selection("ZbCHS3010");
1431 > //   compare_selection("ZbCHS3010");
1432 > //   compare_selection("ZbCHS301010");
1433 > //   compare_selection("ZbCHS3010_SecEta3");
1434 > //   compare_selection("ZbCHS3010_SecEta5");
1435 > //   compare_selection("Zb1510");
1436 > //   compare_selection("ZbCHS301010");
1437 > //   compare_selection("ZbMikko");
1438 > //   compare_selection("ZbCHS3010");
1439 > //   compare_selection("Zb40");
1440  
1441    if(!doquick) {
1442 <    print_all_b_yields();
1443 <    draw_mpf_vars();
1444 <    //  DrawEvilCutFlow();
1442 >    draw_kin_variable("mll","",28,60.0,200.0,"m_{ll}","mll",1,"");//nice and quick test to see if the normalization is ok (i.e. all data is processed etc.)
1443 >    //print_all_b_yields();
1444 > //     draw_mpf_vars();
1445 >    //DrawEvilCutFlow();
1446          
1447 <    draw_Zb_kin_vars();
1447 > //     draw_Zb_kin_vars();
1448          
1449     }
1450    
1451   //  GetNumberEventsInsideOutsideAlphaWindow();
1452 <  new_data_mc_agreement_2d();
1452 >  
1453 >  if(DoCompleteAnalysis) {
1454 >    /* need to run the following scenarios:
1455 >     * - normal
1456 >     * - inclusive
1457 >     * - medium WP
1458 >     * - alpha up
1459 >     * - alpha down
1460 >     */
1461 >    
1462 >    bool fast=true;bool slow=false;//don't change these
1463 >    
1464 >    ZbResultContainer inclusive = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","inclusive",TCut("mll>0"),TCut("1.0"),0.3);
1465 >    ZbResultContainer Reference = new_data_mc_agreement_2d(slow,"ZbCHS3010_alpha","reference",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtT"),0.3);
1466 >    
1467 >    ZbResultContainer NeutrinoQ = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","NeutrinoQ",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtTUp*(1.0*softMuon)"),0.3);
1468 >    ZbResultContainer ANeutrinoQ = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","ANeutrinoQ",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtTUp*(1.0*(1.0-softMuon))"),0.3);
1469 >    
1470 >    ZbResultContainer effmistagUp = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","effmistagUp",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtTUp"),0.3);
1471 >    ZbResultContainer effmistagDown = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","effmistagDown",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtTDown"),0.3);
1472 >    ZbResultContainer medium = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","medium",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.679"),TCut("ZbCHS3010_BTagWgtM"),0.3);
1473 >    ZbResultContainer loose = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","loose",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.244"),TCut("ZbCHS3010_BTagWgtL"),0.3);
1474 >    ZbResultContainer AlphaUp = new_data_mc_agreement_2d(fast,"ZbCHS3010_alphaUp","AlphaUp",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtT"),0.3);
1475 >    ZbResultContainer AlphaDown = new_data_mc_agreement_2d(fast,"ZbCHS3010_alphaDown","AlphaDown",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtT"),0.3);
1476 >    
1477 >    ZbResultContainer L5inclusive = new_data_mc_agreement_2d(fast,"ZbCHS3010_alphaL5","L5inclusive",TCut("mll>0"),TCut("1.0"),0.3);
1478 >    ZbResultContainer L5Reference = new_data_mc_agreement_2d(fast,"ZbCHS3010_alphaL5","L5reference",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtT"),0.3);
1479 >    
1480 >    ZbResultContainer AlphaThresholdVariationDown10 =    new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","AlphaThresholdVariationDown10",   TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),   TCut("ZbCHS3010_BTagWgtT"),0.3-0.1*0.3);
1481 >    ZbResultContainer AlphaThresholdVariationUp10 =      new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","AlphaThresholdVariationUp10",     TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),   TCut("ZbCHS3010_BTagWgtT"),0.3+0.1*0.3);
1482 >    ZbResultContainer AlphaThresholdVariationDown30 =    new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","AlphaThresholdVariationDown30",   TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),   TCut("ZbCHS3010_BTagWgtT"),0.3-0.3*0.3);
1483 >    ZbResultContainer AlphaThresholdVariationUp30 =      new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","AlphaThresholdVariationUp30",     TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),   TCut("ZbCHS3010_BTagWgtT"),0.3+0.3*0.3);
1484 >    ZbResultContainer AlphaThresholdVariationDown10inc = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","AlphaThresholdVariationDown10inc",TCut("mll>0"),                              TCut("1.0"),0.3-0.1*0.3);
1485 >    ZbResultContainer AlphaThresholdVariationUp10inc =   new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","AlphaThresholdVariationUp10inc",  TCut("mll>0"),                              TCut("1.0"),0.3+0.1*0.3);
1486 >    ZbResultContainer AlphaThresholdVariationDown30inc = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","AlphaThresholdVariationDown30inc",TCut("mll>0"),                              TCut("1.0"),0.3-0.3*0.3);
1487 >    ZbResultContainer AlphaThresholdVariationUp30inc =   new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","AlphaThresholdVariationUp30inc",  TCut("mll>0"),                              TCut("1.0"),0.3+0.3*0.3);
1488 >
1489 >    //and now let's compute some systematics!
1490 >    //1 alpha variation
1491 >    float SysMPF_Alpha_down = Reference.MPF_Result_FaceValue.getValue()-AlphaDown.MPF_Result_FaceValue.getValue();
1492 >    float SysMPF_Alpha_up = Reference.MPF_Result_FaceValue.getValue()-AlphaUp.MPF_Result_FaceValue.getValue();
1493 >    
1494 >    //2 btagging efficiency/mistag correction
1495 >    float SysMPF_EFF_down = Reference.MPF_Result_FaceValue.getValue()-effmistagDown.MPF_Result_FaceValue.getValue();
1496 >    float SysMPF_EFF_up = Reference.MPF_Result_FaceValue.getValue()-effmistagUp.MPF_Result_FaceValue.getValue();
1497 >    
1498 >    //3 Presence of soft muons (->neutrino question)
1499 >    float SysMPF_Neutrino_up = Reference.MPF_Result_FaceValue.getValue()-NeutrinoQ.MPF_Result_FaceValue.getValue();
1500 >    float SysMPF_Neutrino_down = Reference.MPF_Result_FaceValue.getValue()-ANeutrinoQ.MPF_Result_FaceValue.getValue();
1501 >    
1502 >    //4 purity
1503 >    zbcanvas->cd();
1504 >    gStyle->SetOptFit(0);
1505 >    write_info(__FUNCTION__,"Don't know the exact purity at the moment, still needs to be determined!");
1506 >    TH1F *purityhisto = new TH1F("purityhisto","purityhisto",1,0,1);
1507 >    purityhisto->GetXaxis()->SetTitle("Purity");
1508 >    purityhisto->GetYaxis()->SetTitle("C_{abs}");
1509 >    purityhisto->GetXaxis()->CenterTitle();
1510 >    purityhisto->GetYaxis()->CenterTitle();
1511 >    
1512 >    TGraphErrors *gSysMPF_Purity = new TGraphErrors();
1513 >    
1514 >    gSysMPF_Purity->SetPoint(0,PurityInclusive,inclusive.MPF_Result_FaceValue.getValue());
1515 >    gSysMPF_Purity->SetPointError(0,0.0,inclusive.MPF_Result_FaceValue.getError());
1516 >    
1517 >    gSysMPF_Purity->SetPoint(1,PurityLoose,loose.MPF_Result_FaceValue.getValue());//x value?
1518 >    gSysMPF_Purity->SetPointError(1,0.0,loose.MPF_Result_FaceValue.getError());
1519 >    gSysMPF_Purity->SetPoint(2,PurityMedium,medium.MPF_Result_FaceValue.getValue());//x value?
1520 >    gSysMPF_Purity->SetPointError(2,0.0,medium.MPF_Result_FaceValue.getError());
1521 >    gSysMPF_Purity->SetPoint(3,PurityTight,Reference.MPF_Result_FaceValue.getValue());//x value?
1522 >    gSysMPF_Purity->SetPointError(3,0.0,Reference.MPF_Result_FaceValue.getError());
1523 >    
1524 >    
1525 >    float min,max;
1526 >    FindMinMax(gSysMPF_Purity,min,max);
1527 >    purityhisto->SetStats(0);
1528 >    purityhisto->GetYaxis()->SetRangeUser(min,max);
1529 >    
1530 >    gSysMPF_Purity->Fit("pol1");
1531 >    TPolyLine *purity_fit_uncert = GetFitUncertaintyShape(gSysMPF_Purity, "pol1", min, max,0.0,1.0);
1532 >    TF1* pufit = (TF1*)gSysMPF_Purity->GetFunction("pol1");
1533 >    
1534 >    TF1* pufit_c = new TF1("pufit_c","[0]+[1]*x");
1535 >    pufit_c->SetParameters(pufit->GetParameters());
1536 >    pufit_c->SetLineColor(TColor::GetColor("#0101DF"));
1537 >    purity_fit_uncert->SetFillColor(TColor::GetColor("#5353FC"));
1538 >    pufit->SetLineColor(TColor::GetColor("#5353FC"));
1539 >    
1540 >    purityhisto->Draw();
1541 >    purity_fit_uncert->Draw("f");
1542 >    gSysMPF_Purity->Draw("P*");
1543 >    pufit_c->Draw("same");
1544 >    DrawPrelim();
1545 >    
1546 >    float ExtrapolatedResult=pufit->GetParameter(0)+1.0*pufit->GetParameter(1);
1547 >    float dExtrapolatedResult=sqrt(pufit->GetParError(0)*pufit->GetParError(0)+1.0*1.0*pufit->GetParError(1)*pufit->GetParError(1));
1548 >    
1549 >    stringstream summary;
1550 >    summary << "#splitline{Reference: " << std::setprecision(4) << Reference.MPF_Result_FaceValue.getValue() << "+/-" << std::setprecision(4) << Reference.MPF_Result_FaceValue.getError() << "}{";
1551 >    summary << "Extrapolation: " << std::setprecision(4) << ExtrapolatedResult << " +/- " << std::setprecision(4) << dExtrapolatedResult << "}";
1552 >    
1553 >    dout << "Have found an extrapolated result at 1.0 of " << ExtrapolatedResult << " ;  will use the difference between this result and the face value as ";
1554 >    dout << "purity-related systematic uncertainty" << endl;
1555 >    
1556 >    TText *infobox = write_title(summary.str());
1557 >    infobox->SetX(0.75);
1558 >    infobox->SetTextSize(0.03);
1559 >    infobox->SetY(0.75);
1560 >    infobox->Draw();
1561 >    
1562 >    CompleteSave(zbcanvas,"Systematics/Purity");
1563 >    
1564 >    
1565 >
1566 >    float SysMPF_Purity = abs(ExtrapolatedResult-Reference.MPF_Result_FaceValue.getValue());
1567 >    
1568 >    float sys_alpha    = abs(SysMPF_Alpha_down)>abs(SysMPF_Alpha_up)?abs(SysMPF_Alpha_down):abs(SysMPF_Alpha_up);
1569 >    float sys_eff      = abs(SysMPF_EFF_down)>abs(SysMPF_EFF_up)?abs(SysMPF_EFF_down):abs(SysMPF_EFF_up);
1570 >    float sys_purity   = abs(SysMPF_Purity);
1571 >    float sys_neutrino = abs(SysMPF_Neutrino_down)>abs(SysMPF_Neutrino_up)?abs(SysMPF_Neutrino_down):abs(SysMPF_Neutrino_up);
1572 >    float sys          = sqrt(sys_alpha*sys_alpha+sys_eff*sys_eff+sys_purity*sys_purity+sys_neutrino*sys_neutrino);
1573 >    
1574 >    dout << "MPF METHOD:" << endl;
1575 >    dout << "   Systematics :    down    up     (selected)" << endl;
1576 >    dout << "      Oversmearing         : " << SysMPF_Alpha_down << "   ,  " << SysMPF_Alpha_up << "   ( " << sys_alpha << ") " << endl;
1577 >    dout << "      eff/mistag variation : " << SysMPF_EFF_down << "   ,  " << SysMPF_EFF_up << "   ( " << sys_eff << ") " << endl;
1578 >    dout << "      purity               : " << SysMPF_Purity << "   ( " << sys_purity << " )" << endl;
1579 >    dout << "      neutrino             : " << SysMPF_Neutrino_down << "   ,  " << SysMPF_Neutrino_up << "   ( " << sys_neutrino << ") " << endl;
1580 >    dout << "      total                : " << sys << "     (note that this is the sys on C_{abs}^{b}, not yet on C_{corr} - there will also be an additional source!)" << endl;
1581 >    dout << endl << endl;
1582 >    dout << "   FINAL RESULTS : " << endl;
1583 >    dout << "     C_{abs}^{b}   = " << Reference.MPF_Result_FaceValue << " (stat) +/- " << sys << " (sys) " << endl;
1584 >    dout << "     C_{abs}^{inc} = " << inclusive.MPF_Result_FaceValue << " (stat) " << endl;
1585 >    
1586 >    float sysfinal = Reference.MPF_Result_FaceValue.getValue()/inclusive.MPF_Result_FaceValue.getValue();
1587 >    sysfinal *= sqrt(sys*sys/(Reference.MPF_Result_FaceValue.getValue()*Reference.MPF_Result_FaceValue.getValue()));
1588 >    //yes, this could be simplified but it would look like a bug.
1589 >    
1590 >
1591 >    
1592 >    
1593 >    dout << "Note: Varying the alpha cut (->face value) we get the following results: " << endl;
1594 >    dout << "-30%: " << AlphaThresholdVariationDown30.MPF_Result_FaceValue/AlphaThresholdVariationDown30inc.MPF_Result_FaceValue << " (stat)" << endl;
1595 >    dout << "-10%: " << AlphaThresholdVariationDown10.MPF_Result_FaceValue/AlphaThresholdVariationDown10inc.MPF_Result_FaceValue << " (stat)" << endl;
1596 >    dout << "  0%: " << Reference.MPF_Result_FaceValue/inclusive.MPF_Result_FaceValue << " (stat) +/- " << endl;
1597 >    dout << "+10%: " << AlphaThresholdVariationUp10.MPF_Result_FaceValue/AlphaThresholdVariationUp10inc.MPF_Result_FaceValue << " (stat)" << endl;
1598 >    dout << "+30%: " << AlphaThresholdVariationUp30.MPF_Result_FaceValue/AlphaThresholdVariationUp30inc.MPF_Result_FaceValue << " (stat)" << endl;
1599 >    
1600 >    dout << "          XCheck for Rabs, varying the alpha cut (->face value) we get the following results: " << endl;
1601 >    dout << "             -30%: " << AlphaThresholdVariationDown30.Rabs_Result_FaceValue/AlphaThresholdVariationDown30inc.Rabs_Result_FaceValue << " (stat)" << endl;
1602 >    dout << "             -10%: " << AlphaThresholdVariationDown10.Rabs_Result_FaceValue/AlphaThresholdVariationDown10inc.Rabs_Result_FaceValue << " (stat)" << endl;
1603 >    dout << "               0%: " << Reference.Rabs_Result_FaceValue/inclusive.Rabs_Result_FaceValue << " (stat) +/- " << endl;
1604 >    dout << "             +10%: " << AlphaThresholdVariationUp10.Rabs_Result_FaceValue/AlphaThresholdVariationUp10inc.Rabs_Result_FaceValue << " (stat)" << endl;
1605 >    dout << "             +30%: " << AlphaThresholdVariationUp30.Rabs_Result_FaceValue/AlphaThresholdVariationUp30inc.Rabs_Result_FaceValue << " (stat)" << endl;
1606 >
1607 >    
1608 >    float sys_alphavar_up   = abs(AlphaThresholdVariationUp30.MPF_Result_FaceValue.getValue()/AlphaThresholdVariationUp30inc.MPF_Result_FaceValue.getValue() - Reference.MPF_Result_FaceValue.getValue()/inclusive.MPF_Result_FaceValue.getValue());
1609 >    float sys_alphavar_down = abs(AlphaThresholdVariationDown30.MPF_Result_FaceValue.getValue()/AlphaThresholdVariationDown30inc.MPF_Result_FaceValue.getValue() - Reference.MPF_Result_FaceValue.getValue()/inclusive.MPF_Result_FaceValue.getValue());
1610 >    float sys_alphavar = sys_alphavar_up>sys_alphavar_down?sys_alphavar_up:sys_alphavar_down;
1611 >    
1612 >    dout << " -> Additional systematic (from 30%) : " << sys_alphavar << endl;
1613 >    sysfinal=sqrt(sysfinal*sysfinal+sys_alphavar*sys_alphavar);
1614 >    dout << "Final systematic is : " << sysfinal << endl;
1615 >    
1616 >    dout << "   ******************************************** <  FINAL RESULT > ******************************************** " << endl;
1617 >    dout << "     C_{corr}      = " << Reference.MPF_Result_FaceValue/inclusive.MPF_Result_FaceValue << " (stat) +/- " << sysfinal << " (sys) " << endl;
1618 >    dout << "   ******************************************** < /FINAL RESULT > ******************************************** " << endl;
1619 >    
1620 >    dout << "Note that if L5 corrections are applied, we get " << L5Reference.MPF_Result_FaceValue/L5inclusive.MPF_Result_FaceValue << " (stat) " << endl;
1621 >    
1622 >    LeadingB=TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898");
1623 >    
1624 >  } else {
1625 >    ZbResultContainer Reference = new_data_mc_agreement_2d(false,"ZbCHS3010_alpha","reference",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtT"),0.3);
1626 >  }
1627 >  
1628    
1629    delete zbcanvas;
1630   }
# Line 1177 | Line 1650 | void SpecialCutFlow(string identifier) {
1650    MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB) << ")";
1651    MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB) << ")";
1652    MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut) << ")";
1653 <  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000").c_str())) << ")";
1654 <  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.3").c_str())) << ")";
1655 <  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.2").c_str())) << ")";
1656 <  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.15").c_str())) << ")";
1657 <  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.1").c_str())) << ")";
1658 <  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>10&&pt<1000&&"+identifier+"_alpha<0.05").c_str())) << ")";
1653 >  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>30&&pt<1000").c_str())) << ")";
1654 >  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>30&&pt<1000&&"+identifier+"_alpha<0.3").c_str())) << ")";
1655 >  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>30&&pt<1000&&"+identifier+"_alpha<0.2").c_str())) << ")";
1656 >  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>30&&pt<1000&&"+identifier+"_alpha<0.15").c_str())) << ")";
1657 >  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>30&&pt<1000&&"+identifier+"_alpha<0.1").c_str())) << ")";
1658 >  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>30&&pt<1000&&"+identifier+"_alpha<0.05").c_str())) << ")";
1659    MegaCut << " ";
1660    
1661    TCanvas *can = new TCanvas("can","can");
# Line 1195 | Line 1668 | void SpecialCutFlow(string identifier) {
1668    
1669    
1670    TH1F *data  = allsamples.Draw("data",    MegaCut.str(),11,-0.5,10.5, "", "events", basecut,data,luminosity);
1671 <   THStack mcm = allsamples.DrawStack("mc", MegaCut.str(),11,-0.5,10.5, "", "events", basecut,mc,luminosity);
1671 >  THStack mcm = allsamples.DrawStack("mc", MegaCut.str(),11,-0.5,10.5, "", "events", basecut,mc,luminosity);
1672  
1673    float runningsum=0;
1674    for(int i=data->GetNbinsX();i>0;i--) {
# Line 1233 | Line 1706 | void SpecialCutFlow(string identifier) {
1706    data->GetXaxis()->SetBinLabel(3,"Leading jet is b-jet");
1707    data->GetXaxis()->SetBinLabel(4,"|#eta_{b}|<1.3 ");
1708    data->GetXaxis()->SetBinLabel(5,"|#delta#phi(b,Z)|>2.7");
1709 <  data->GetXaxis()->SetBinLabel(6,"10<p_{T}^{Z}<1000 GeV");
1709 >  data->GetXaxis()->SetBinLabel(6,"30<p_{T}^{Z}<1000 GeV");
1710    data->GetXaxis()->SetBinLabel(7,"#alpha < 0.3");
1711    data->GetXaxis()->SetBinLabel(8,"#alpha < 0.2");
1712    data->GetXaxis()->SetBinLabel(9,"#alpha < 0.15");
# Line 1251 | Line 1724 | void SpecialCutFlow(string identifier) {
1724    leg->AddEntry(crap,neventsinfo.str().c_str(),"l");
1725    
1726    data->Draw("e1");
1727 <  mcm.Draw("same");
1727 >  mcm.Draw("histo,same");
1728    data->Draw("same,e1,axis");
1729    data->Draw("same,e1");
1730   //  data->Draw("same,TEXT");
# Line 1264 | Line 1737 | void SpecialCutFlow(string identifier) {
1737    
1738    delete crap;
1739    delete data;
1740 <  
1740 >  CleanLegends();
1741 >  DeleteStack(mcm);
1742 >
1743    
1744   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines