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.19 by buchmann, Mon Jan 7 14:42:35 2013 UTC vs.
Revision 1.29 by buchmann, Fri Apr 26 07:18:47 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 16 | Line 17
17   #include <TProfile.h>
18   #include <TLegendEntry.h>
19  
20 + #include <RooRealVar.h>
21 + #include <RooArgSet.h>
22 + #include <RooDataSet.h>
23 + #include <RooDataHist.h>
24 + #include <RooMCStudy.h>
25 + #include <RooCategory.h>
26 +
27 + #include <RooPlot.h>
28 + #include <RooSimultaneous.h>
29 + #include <RooAddPdf.h>
30 + #include <RooFitResult.h>
31 + #include <RooVoigtian.h>
32 + #include <RooChebychev.h>
33 + #include <RooMsgService.h>
34 +
35 + #include <RooStats/ModelConfig.h>
36 + #include "RooStats/ProfileLikelihoodCalculator.h"
37 + #include "RooStats/LikelihoodInterval.h"
38 + #include "RooStats/HypoTestResult.h"
39 + #include "RooStats/SimpleLikelihoodRatioTestStat.h"
40 + #include "RooStats/ProfileLikelihoodTestStat.h"
41 + #include "RooStats/HybridCalculatorOriginal.h"
42 + #include "RooStats/HypoTestInverterOriginal.h"
43 +
44 +
45 +
46   using namespace std;
47  
48   using namespace PlottingSetup;
49  
50   namespace ZbData {
24 – vector<float> data_over_mc;
25 –
51   TCut ZplusBsel("pt1>20&&pt2>20&&mll>60&&mll<120&&id1==id2");
52   TCut LeadingB("ZbCHS3010_bTagProbCSVBP[0]>0.898");
53   TCut EtaB("abs(ZbCHS3010_pfJetGoodEta[0])<1.3");
54   TCut PhiZcut("abs(ZbCHS3010_pfJetDphiZ[0])>2.7");
55 <  
55 > float LowerHistoBoundary=0.5; // where to cut on MPF (or Rabs)
56 > float UpperHistoBoundary=1.5; // where to cut on MPF (or Rabs)
57 >
58   }
59  
60   TCut STDWEIGHT;
# Line 61 | Line 88 | class ZbResultContainer {
88   public:
89    string name;  
90    vector<ZbPointResult> MPF_Rabs_location;
91 +  
92    Value MPF_Result_FaceValue;
93    Value Rabs_Result_FaceValue;
94    Value MPF_Result_Extrapolated;
95    Value Rabs_Result_Extrapolated;
96    
97 +  vector<Value> MPF_Result_FaceValue_Bins;
98 +  vector<Value> Rabs_Result_FaceValue_Bins;
99 +  vector<Value> MPF_Result_Extrapolated_Bins;
100 +  vector<Value> Rabs_Result_Extrapolated_Bins;
101 +  
102 +  void Print();
103 +  
104    ZbResultContainer(string name);
70 –  ZbResultContainer(const ZbResultContainer &old);
105   };
106  
107  
108   ZbResultContainer::ZbResultContainer(string _name) {
109    name=_name;
110 <  cout<< "Result Container for \"" << name << "\" has been initialized" << endl;
77 <  
110 >  dout<< "Result Container for \"" << name << "\" has been initialized" << endl;
111   }
112  
113 < std::ostream &operator<<(std::ostream &ostr, ZbResultContainer &b)
113 > void ZbResultContainer::Print()
114   {
115 <  ostr << "Results for " << b.name << endl;
116 <  ostr << " MPF and Rabs locations: " << endl;
117 <  for(int i=0;i<b.MPF_Rabs_location.size();i++) ostr << "   " << b.MPF_Rabs_location[i] << endl;
118 <  ostr << " Results: " << endl;
119 <  ostr << "   Face Value: " << endl;
120 <  ostr << "     MPF " << b.MPF_Result_FaceValue << endl;
121 <  ostr << "     Rabs " << b.Rabs_Result_FaceValue << endl;
122 <  ostr << "   Extrapolated: " << endl;
123 <  ostr << "     MPF " << b.MPF_Result_Extrapolated << endl;
124 <  ostr << "             Rabs " << b.Rabs_Result_Extrapolated << endl;
125 <  
126 <  return ostr;
115 >  dout << "Results for " << this->name << endl;
116 >  dout << " MPF and Rabs locations: " << endl;
117 >  for(int i=0;i<this->MPF_Rabs_location.size();i++) dout << "   " << this->MPF_Rabs_location[i] << endl;
118 >  dout << " Results: " << endl;
119 >  dout << "   Face Value: " << endl;
120 >  dout << "     MPF " << this->MPF_Result_FaceValue << endl;
121 >  dout << "        Bins: [ ";
122 >  if(MPF_Result_FaceValue_Bins.size()>0) {
123 >    for(int i=0;i<MPF_Result_FaceValue_Bins.size()-1;i++) dout << MPF_Result_FaceValue_Bins[i] << " , ";
124 >    if(MPF_Result_FaceValue_Bins.size()>0) dout << MPF_Result_FaceValue_Bins[MPF_Result_FaceValue_Bins.size()-1] << " ]" << endl;                                                                                                  
125 >  }
126 >  
127 >  dout << "     Rabs " << this->Rabs_Result_FaceValue << endl;
128 >  if(Rabs_Result_FaceValue_Bins.size()>0) {
129 >    dout << "        Bins: [ ";                                                                                                                                                    
130 >    for(int i=0;i<Rabs_Result_FaceValue_Bins.size()-1;i++) dout << Rabs_Result_FaceValue_Bins[i] << " , ";                                                                          
131 >    if(Rabs_Result_FaceValue_Bins.size()>0) dout << Rabs_Result_FaceValue_Bins[Rabs_Result_FaceValue_Bins.size()-1] << " ]" << endl;                                                                                                
132 >  }
133 >  
134 >  dout << "   Extrapolated: " << endl;
135 >  dout << "     MPF " << this->MPF_Result_Extrapolated << endl;
136 >  if(MPF_Result_Extrapolated_Bins.size()>0) {
137 >    dout << "        Bins: [ ";
138 >    for(int i=0;i<MPF_Result_Extrapolated_Bins.size()-1;i++) dout << MPF_Result_Extrapolated_Bins[i] << " , ";
139 >    if(MPF_Result_Extrapolated_Bins.size()>0) dout << MPF_Result_Extrapolated_Bins[MPF_Result_Extrapolated_Bins.size()-1] << " ]" << endl;
140 >  }
141 >  
142 >  dout << "             Rabs " << this->Rabs_Result_Extrapolated << endl;
143 >  if(Rabs_Result_Extrapolated_Bins.size()>0) {
144 >    dout << "        Bins: [ ";
145 >    for(int i=0;i<Rabs_Result_Extrapolated_Bins.size()-1;i++) dout << Rabs_Result_Extrapolated_Bins[i] << " , ";
146 >    if(Rabs_Result_Extrapolated_Bins.size()>0) dout << Rabs_Result_Extrapolated_Bins[Rabs_Result_Extrapolated_Bins.size()-1] << " ]" << endl;
147 >  }
148   }
149  
150  
97 –
151   using namespace ZbData;
152  
153  
# Line 104 | Line 157 | void SpecialCutFlow(string identifier);
157   //*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*            AND THIS LINE             /*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*//
158  
159  
160 < void print_yield(TCut cut) {
160 > class CutYields{
161 > public:
162 >  float data;
163 >  float Zb;
164 >  float Zc;
165 >  float Zl;
166 >  float tt;
167 >  float SiTo;
168 >  float Dibosons;
169 >  float WJets;
170 >  float Other;
171 >  
172 >  CutYields();
173 > };
174 >
175 > CutYields::CutYields() {
176 >  this->data=0;
177 >  this->Zb=0;
178 >  this->Zc=0;
179 >  this->Zl=0;
180 >  this->tt=0;
181 >  this->SiTo=0;
182 >  this->Dibosons=0;
183 >  this->WJets=0;
184 >  this->Other=0;
185 > }
186 >
187 > std::ostream &operator<<(std::ostream &ostr, CutYields y)
188 > {
189 >  return ostr << y.data << ";" << y.Zb << ";" << y.Zc << ";" << y.Zl << ";" << y.tt << ";" << y.SiTo << ";" << y.Dibosons << ";" << y.WJets << ";" << endl;
190 > }
191 >
192 >
193 > CutYields print_yield(TCut cut) {
194 >  
195 >  CutYields Yield;
196    TH1F *data  = allsamples.Draw("data",    "mll",1,0,10000000, "m_{ll} [GeV]", "events", cut,data,luminosity);
197    dout << "  Yields for " << (const char*) cut << endl;
198 <  dout << "     data: " << data->Integral() << endl;
198 > //  dout << "     data: " << data->Integral() << endl;
199 >  
200 >  Yield.data=data->Integral();
201    THStack mcm = allsamples.DrawStack("mc", "mll",1,0,10000000, "m_{ll} [GeV]", "events", cut,mc,luminosity);
202    int nhists = mcm.GetHists()->GetEntries();
113 –  dout << "Going to loop over " << nhists << " histograms " << endl;
203    TFile *test = new TFile("test.root","RECREATE");
204    mcm.Write();
205    test->Close();
206    for (int i = 0; i < nhists; i++) {
207      TH1* hist = static_cast<TH1*>(mcm.GetHists()->At(i));
208 <    cout << "     " << hist->GetName() << ": " << hist->Integral() << endl;
208 >    string name=hist->GetName();
209 >    if(Contains(name,"t_bar_t")) Yield.tt+=hist->Integral();
210 >    if(Contains(name,"W_Jets")) Yield.WJets+=hist->Integral();
211 >    if(Contains(name,"Single_top")) Yield.SiTo+=hist->Integral();
212 >    if(Contains(name,"Dibosons")) Yield.Dibosons+=hist->Integral();
213 >    if(Contains(name,"Z_l_")) Yield.Zl+=hist->Integral();
214 >    if(Contains(name,"Z_c_")) Yield.Zc+=hist->Integral();
215 >    if(Contains(name,"Z_b_")) Yield.Zb+=hist->Integral();
216    }
217 <    
218 <  cout << endl << endl;
217 >  
218 >  dout << Yield << endl;
219 >  dout << endl << endl;
220    
221    delete data;
222 +  CleanLegends();
223 +  DeleteStack(mcm);
224 +
225 +  return Yield;
226   }
227  
228   void draw_kin_variable(string variable, TCut cut, int nbins, float min, float max, string xlabel, string saveas, bool dologscale, string speciallabel="") {
128 –  
229    TCanvas *ckin = new TCanvas("ckin","kin variable canvas");
230 < //  ckin->SetLogy(dologscale);
230 >  if(dologscale) ckin->SetLogy(1);
231    TH1F *datahisto = allsamples.Draw("datahisto",variable,nbins,min,max, xlabel, "events",cut,data,luminosity);
232    datahisto->SetMarkerSize(DataMarkerSize);
233    THStack mcstack = allsamples.DrawStack("mcstack",variable,nbins,min,max, xlabel, "events",cut,mc,luminosity);
234    if(dologscale) datahisto->SetMaximum(3*datahisto->GetMaximum());
235    else datahisto->SetMaximum(1.5*datahisto->GetMaximum());
236 <
236 >  
237 >  float SumMC=0.0;
238 >  TIter nextH(mcstack.GetHists());
239 >  TH1F* h;
240 >  while ( h = (TH1F*)nextH() ) {
241 >    SumMC += h->Integral();
242 >  }
243 >  float scale = datahisto->Integral()/SumMC;
244 >
245 >  // Re-scale stack to data integral
246 >  nextH = TIter(mcstack.GetHists());
247 >  
248 >  while ( h = (TH1F*)nextH() ) {
249 >    h->Scale(scale);
250 >  }
251 >  
252    TText *speciall;
253    datahisto->Draw("e1");
254    ckin->Update();
255    mcstack.Draw("histo,same");
256    datahisto->Draw("same,e1");
257    TLegend *kinleg = allsamples.allbglegend();
258 +
259    kinleg->Draw();
260  
261    TPad *kinpad = new TPad("kinpad","kinpad",0,0,1,1);
262    kinpad->cd();
263 <  kinpad->SetLogy(dologscale);
263 >  if(dologscale) kinpad->SetLogy(1);
264    datahisto->Draw("e1");
265    mcstack.Draw("histo,same");
266    datahisto->Draw("same,e1");
# Line 152 | Line 268 | void draw_kin_variable(string variable,
268    kinleg->Draw();
269    DrawPrelim();
270    saveas="kin/"+saveas;
271 +  kinpad->cd();
272 +  
273    if(speciallabel!="") {
274 <    speciall = new TText(0.2,0.8,speciallabel.c_str());
274 >    speciall = write_title(speciallabel);
275 >    speciall->SetX(0.18);
276 >    speciall->SetTextAlign(11);
277 >    speciall->SetTextSize(0.03);
278 >    speciall->SetY(0.85);
279      speciall->Draw();
280    }
159 –  save_with_ratio(datahisto,mcstack,kinpad->cd(),saveas);
281    
282 <  
162 <  ZbData::data_over_mc.push_back(datahisto->Integral()/CollapseStack(mcstack)->Integral());
282 >  Save_With_Ratio(datahisto,mcstack,kinpad->cd(),saveas);
283    
284    datahisto->Delete();
285 +  delete kinleg;
286    delete ckin;
287 +  CleanLegends();
288 +  DeleteStack(mcstack);
289   }
290 <  
290 >
291   void print_all_b_yields() {
292 <  cout << "Basic selection with a b jet" << endl;
293 <  print_yield(ZplusBsel&&"ZbCHS3010_pfJetGoodNumBtag>0");
294 <  cout << "Leading jet is a b-jet " << endl;
295 <  print_yield(ZplusBsel&&LeadingB);
296 <  cout << "|#eta_{b}|<1.3 " << endl;
297 <  print_yield(ZplusBsel&&LeadingB&&EtaB);
298 <  cout << "|#delta#phi(b<Z)|>2.7" << endl;
299 <  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut);
300 <  cout << "30<ptZ<1000 GeV" << endl;
301 <  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000");
302 <  cout << "#alpha < 0.35" << endl;
303 <  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"(ZbCHS3010_alpha)<0.35");
304 <  cout << "#alpha < 0.3" << endl;
305 <  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"(ZbCHS3010_alpha)<0.3");
306 <  cout << "#alpha < 0.2" << endl;
307 <  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"(ZbCHS3010_alpha)<0.2");
308 <  cout << "#alpha < 0.15" << endl;
309 <  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"(ZbCHS3010_alpha)<0.15");
310 <  cout << "#alpha < 0.1" << endl;
311 <  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"(ZbCHS3010_alpha)<0.1");
292 >  vector<CutYields> yields;
293 >  dout << "Basic selection with a b jet" << endl;
294 >  yields.push_back(print_yield(ZplusBsel&&"ZbCHS3010_pfJetGoodNumBtag>0"));
295 >  dout << "Leading jet is a b-jet " << endl;
296 >  yields.push_back(print_yield(ZplusBsel&&LeadingB));
297 >  dout << "|#eta_{b}|<1.3 " << endl;
298 >  yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB));
299 >  dout << "|#delta#phi(b<Z)|>2.7" << endl;
300 >  yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut));
301 >  dout << "30<ptZ<1000 GeV" << endl;
302 >  yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"));
303 >  dout << "#alpha < 0.35" << endl;
304 >  yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"(ZbCHS3010_alpha)<0.35"));
305 >  dout << "#alpha < 0.3" << endl;
306 >  yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"(ZbCHS3010_alpha)<0.3"));
307 > /*  dout << "#alpha < 0.2" << endl;
308 >  yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"(ZbCHS3010_alpha)<0.2"));
309 >  dout << "#alpha < 0.15" << endl;
310 >  yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"(ZbCHS3010_alpha)<0.15"));
311 >  dout << "#alpha < 0.1" << endl;
312 >  yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"(ZbCHS3010_alpha)<0.1"));*/
313 >
314 >  for(int iy=0;iy<yields.size();iy++) {
315 >    dout << yields[iy] << endl;
316 >  }
317   }  
318  
319   void DrawEvilCutFlow() {
# Line 213 | Line 341 | void DrawEvilCutFlow() {
341    
342    
343    TH1F *data  = allsamples.Draw("data",    MegaCut.str(),11,-0.5,10.5, "", "events", basecut,data,luminosity);
344 <   THStack mcm = allsamples.DrawStack("mc", MegaCut.str(),11,-0.5,10.5, "", "events", basecut,mc,luminosity);
344 >  THStack mcm = allsamples.DrawStack("mc", MegaCut.str(),11,-0.5,10.5, "", "events", basecut,mc,luminosity);
345  
346    float runningsum=0;
347    for(int i=data->GetNbinsX();i>0;i--) {
# Line 263 | Line 391 | void DrawEvilCutFlow() {
391    DrawPrelim();
392    
393    CompleteSave(can,"CutFlow");
394 +  CleanLegends();
395 +  DeleteStack(mcm);
396 +
397   }
398  
399   void draw_Zb_kin_vars() {
400 <   draw_kin_variable("pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000",25,0,200,"Z p_{T}","Official/Zpt",1);
401 <   draw_kin_variable("ZbCHS3010_pfJetGoodPt[1]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000",25,0,200,"Sub-Leading Jet Pt","Official/Jet2Pt",1);
402 <   draw_kin_variable("ZbCHS3010_pfJetGoodPt[0]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000",25,0,200,"Leading Jet Pt (B)","Official/LeadingJetPt",1);
403 <   draw_kin_variable("ZbCHS3010_pfJetGoodPt[1]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000",20,0,2,"p_{T}^{J2} / p_{T}^{Z}","Official/SubLeadingJetPt_Over_ZPt",1);
404 <   draw_kin_variable("ZbCHS3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<50",   40,0,2,"#alpha","Official/alpha_pt_30_to_50",1);
405 <   draw_kin_variable("ZbCHS3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>50&&pt<75",  40,0,2,"#alpha","Official/alpha_50_to_75",1);
406 <   draw_kin_variable("ZbCHS3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>75&&pt<125", 40,0,2,"#alpha","Official/alpha_pt_75_to_125",1);
407 <   draw_kin_variable("ZbCHS3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>125&&pt<1000",40,0,2,"#alpha","Official/alpha_pt_125_to_1000",1);
408 <    draw_kin_variable("ZbCHS3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000",40,0,2,"#alpha","Official/alpha_pt_30_to_1000",1);
409 <  
410 <  draw_kin_variable("ZbCHS3010_pfBJetDphiZ[0]",ZplusBsel&&LeadingB&&"pt>30&&pt<1000&&pfBJetDphiZ[0]>0",30,0,3.2,"#delta#phi (Z,b lead)","DeltaPhiZBlead",0);
411 <  
412 < /*
413 <  draw_kin_variable("pt1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000",20,0,200,"p_{T}^{l1}","Lep/pt1",1);
414 <  draw_kin_variable("pt1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000&&id1==0",20,0,200,"p_{T}^{e1}","Lep/pt1_e",1);
415 <  draw_kin_variable("pt1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000&&id1==1",20,0,200,"p_{T}^{#mu1}","Lep/pt1_m",1);
416 <  draw_kin_variable("pt2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000",20,0,200,"p_{T}^{l2}","Lep/pt2",1);
417 <  draw_kin_variable("pt2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000&&id2==0",20,0,200,"p_{T}^{e2}","Lep/pt2_e",1);
418 <  draw_kin_variable("pt2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000&&id2==1",20,0,200,"p_{T}^{#mu2}","Lep/pt2_m",1);
419 <  draw_kin_variable("eta1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000",20,-3.1,3.1,"#eta^{l1}","Lep/eta1",1);
420 <  draw_kin_variable("eta1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000&&id1==0",20,-3.1,3.1,"#eta^{e1}","Lep/eta1_e",1);
421 <  draw_kin_variable("eta1",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000&&id1==1",20,-3.1,3.1,"#eta^{m1}","Lep/eta1_m",1);
422 <  draw_kin_variable("eta2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000",20,-3.1,3.1,"#eta^{l1}","Lep/eta2",1);
423 <  draw_kin_variable("eta2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000&&id2==0",20,-3.1,3.1,"#eta^{e2}","Lep/eta2_e",1);
424 <  draw_kin_variable("eta2",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000&&id2==1",20,-3.1,3.1,"#eta^{#mu2}","Lep/eta2_m",1);
425 <  draw_kin_variable("numVtx",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000",30,0,30,"N(Vtx)","Lep/nVtx",1);
426 <  */
400 >  
401 >   draw_kin_variable("pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000"),25,0,200,"Z p_{T}","Official/Zpt",1);
402 >   draw_kin_variable("pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&id1==0"),25,0,200,"Z p_{T}","Official/Zpt_ee",1);
403 >   draw_kin_variable("pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&id1==1"),25,0,200,"Z p_{T}","Official/Zpt_mm",1);
404 >   draw_kin_variable("pt",ZplusBsel&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000"),25,0,200,"Z p_{T}","Official/Zpt__nonBtagged",1);
405 >   draw_kin_variable("pt",ZplusBsel&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&id1==0"),25,0,200,"Z p_{T}","Official/Zpt__nonBtagged_ee",1);
406 >   draw_kin_variable("pt",ZplusBsel&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&id1==1"),25,0,200,"Z p_{T}","Official/Zpt__nonBtagged_mm",1);
407 >  
408 >   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);
409 >   draw_kin_variable("ZbCHS3010_pfJetGoodPt[1]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&id1==0"),25,0,200,"Sub-Leading Jet Pt","Official/Jet2Pt_ee",1);
410 >   draw_kin_variable("ZbCHS3010_pfJetGoodPt[1]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&id1==1"),25,0,200,"Sub-Leading Jet Pt","Official/Jet2Pt_mm",1);
411 >   draw_kin_variable("ZbCHS3010_pfJetGoodPt[1]",ZplusBsel&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000"),25,0,200,"Sub-Leading Jet Pt","Official/Jet2Pt__nonBtagged",1);
412 >   draw_kin_variable("ZbCHS3010_pfJetGoodPt[1]",ZplusBsel&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&id1==0"),25,0,200,"Sub-Leading Jet Pt","Official/Jet2Pt__nonBtagged_ee",1);
413 >   draw_kin_variable("ZbCHS3010_pfJetGoodPt[1]",ZplusBsel&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&id1==1"),25,0,200,"Sub-Leading Jet Pt","Official/Jet2Pt__nonBtagged_mm",1);
414 >  
415 >   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);
416 >   draw_kin_variable("ZbCHS3010_pfJetGoodPt[0]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&id1==0"),25,0,200,"Leading Jet Pt (B)","Official/LeadingJetPt_ee",1);
417 >   draw_kin_variable("ZbCHS3010_pfJetGoodPt[0]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&id1==1"),25,0,200,"Leading Jet Pt (B)","Official/LeadingJetPt_mm",1);
418 >   draw_kin_variable("ZbCHS3010_pfJetGoodPt[0]",ZplusBsel&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000"),25,0,200,"Leading Jet Pt (B)","Official/LeadingJetPt__nonBtagged",1);
419 >   draw_kin_variable("ZbCHS3010_pfJetGoodPt[0]",ZplusBsel&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&id1==0"),25,0,200,"Leading Jet Pt (B)","Official/LeadingJetPt__nonBtagged_ee",1);
420 >   draw_kin_variable("ZbCHS3010_pfJetGoodPt[0]",ZplusBsel&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&id1==1"),25,0,200,"Leading Jet Pt (B)","Official/LeadingJetPt__nonBtagged_mm",1);
421 >  
422 >  
423 >   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);
424 >   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);
425 >   draw_kin_variable("ZbCHS3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>50&&pt<75"),  40,0,2,"#alpha","Official/alpha_50_to_75",1);
426 >   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);
427 >   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);
428 >   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);
429 >  
430 >   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);
431   }
432  
433   void draw_mpf_vars() {
434 <  ZbData::data_over_mc.clear();
435 <  draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"ZbCHS3010_alpha<0.3",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p3",0,"#alpha<0.3, 30 GeV < p_{T}^{Z} < 1000 GeV");
436 <  draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"ZbCHS3010_alpha<0.2",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p2",0,"#alpha<0.2, 30 GeV < p_{T}^{Z} < 1000 GeV");
437 <  draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"ZbCHS3010_alpha<0.15",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p15",0,"#alpha<0.15, 30 GeV < p_{T}^{Z} < 1000 GeV");
438 <  draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"ZbCHS3010_alpha<0.1",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p1",0,"#alpha<0.1, 30 GeV < p_{T}^{Z} < 1000 GeV");
439 <  draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"ZbCHS3010_alpha<0.05",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p05",0,"#alpha<0.05, 30 GeV < p_{T}^{Z} < 1000 GeV");
440 <  
441 <  draw_kin_variable("ZbCHS3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"ZbCHS3010_alpha<0.3",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p3",0,"#alpha<0.3, 30 GeV < p_{T}^{Z} < 1000 GeV");
442 <  draw_kin_variable("ZbCHS3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"ZbCHS3010_alpha<0.2",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p2",0,"#alpha<0.2, 30 GeV < p_{T}^{Z} < 1000 GeV");
443 <  draw_kin_variable("ZbCHS3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"ZbCHS3010_alpha<0.15",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p15",0,"#alpha<0.15, 30 GeV < p_{T}^{Z} < 1000 GeV");
444 <  draw_kin_variable("ZbCHS3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"ZbCHS3010_alpha<0.1",20,0,2,"p_{T} b-jet / p_{T} Z","ptb_over_ptz___alpha_smaller_0p1",0,"#alpha<0.1, 30 GeV < p_{T}^{Z} < 1000 GeV");
310 <  draw_kin_variable("ZbCHS3010_pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"ZbCHS3010_alpha<0.05",20,0,2,"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");
434 >  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");
435 >  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");
436 >  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");
437 >  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");
438 >  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");
439 >  
440 >  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");
441 >  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");
442 >  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");
443 >  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");
444 >  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");
445   }  
446  
447 < /*double GetMedian(TH1F *histo) {
448 <    int numBins = histo->GetXaxis()->GetNbins();
449 <    Double_t x[numBins];
450 <    Double_t y[numBins];
451 <    for (int i = 1; i <= numBins; i++) {
452 <        x[i] = histo->GetBinCenter(i);
453 <        y[i] = histo->GetBinContent(i);
454 <    }
455 <    return TMath::Median(numBins, &x, &y);
456 < }*/
447 > void make_RooFit(string ContainerName, TH1F *h, float &peak, float &error_peak, string saveas) {
448 >  TCanvas *fcan = new TCanvas("fcan","fcan");
449 >  RooRealVar mpf("mpf","mpf", -5, 5, "");
450 >  RooDataHist AllData("AllData", "CMS Real Dataset", mpf , h);
451 >
452 >  RooRealVar fraction("fraction", "fraction", 0.8, 0, 1);
453 >  RooRealVar mean("mean", "mean", 1, 0, 3);
454 >  RooRealVar sigma("sigma", "sigma", 1, 0, 5);
455 >  RooRealVar width("width", "width", 2.4, 0, 5);
456 >  RooRealVar chevvar1("chevvar1", "chevvar1", 0.2,-1.0, 1.0);
457 >  RooRealVar chevvar2("chevvar2", "chevvar2", 0.2,-1.0, 1.0);
458 >
459 >  RooVoigtian signal("signal", "signal", mpf, mean, width, sigma);
460 >  RooChebychev background("background_", "background", mpf, RooArgList(chevvar1, chevvar2));
461 >
462 >  RooAddPdf pdf("pdf",  "pdf", RooArgList(signal, background), fraction);
463 >
464 >  RooFitResult *result = pdf.fitTo(AllData, RooFit::Save());
465 >  RooPlot* frame = mpf.frame(RooFit::Bins(30), RooFit::Title("MPF")) ;
466 >  AllData.plotOn(frame);
467 >  pdf.plotOn(frame, RooFit::LineColor(kBlue)) ;
468 >  pdf.plotOn(frame, RooFit::Components(RooArgSet(background)), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed)) ;
469 >  pdf.plotOn(frame, RooFit::Components(RooArgSet(signal)), RooFit::LineStyle(kDashed), RooFit::LineColor(kGreen)) ;
470 >
471 >  RooRealVar* thePeak = (RooRealVar*) result->floatParsFinal().find("mean");
472 >  peak = thePeak->getVal();
473 >  error_peak = thePeak->getError();
474 >  
475 >  ///******************* just making it beautiful ....
476 >  RooPlot* frame2 = mpf.frame(RooFit::Bins(30), RooFit::Title("MPF")) ;
477 >  TH1F *hcopy = (TH1F*)h->Clone("hcopy");
478 >  hcopy->Rebin(40);
479 >  RooDataHist AllRData("AllRData", "CMS Dataset", mpf , hcopy);
480 >  
481 >  AllRData.plotOn(frame2);
482 >  pdf.plotOn(frame2, RooFit::LineColor(kBlue)) ;
483 >  pdf.plotOn(frame2, RooFit::Components(RooArgSet(background)), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed)) ;
484 >  pdf.plotOn(frame2, RooFit::Components(RooArgSet(signal)), RooFit::LineStyle(kDashed), RooFit::LineColor(kGreen)) ;
485 >
486 >  string saveasFinal = saveas + "_Fit";
487 >  fcan->cd() ; gPad->SetLeftMargin(0.15) ; frame->GetYaxis()->SetTitleOffset(1.4) ; frame2->Draw();
488 >  
489 >  if(Contains(saveas,"_data")) DrawPrelim();
490 >  else DrawMCPrelim();
491 >
492 >  CompleteSave(fcan, ContainerName+"/Fit/"+saveasFinal);
493 >  
494 >  delete hcopy;
495 >  delete fcan;
496 >  delete frame;
497 >  delete thePeak;
498 > }
499  
500 + float bother=0.1;
501 + float bother2=1;
502   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) {
325 – //  write_warning(__FUNCTION__,"Debugging this function, therefore always returning 3.1415+/-1.010101"); return Value(3.1415,1.010101);
503    TCanvas *cn = new TCanvas("cn","cn");
504    string varname="MPF";
505    
506 +  if(a2<15.0/pt2 && a1>0) dout << "this bin ( " << a1 << " < alpha < " << a2 <<" and " << pt1 << " < p_{T}^{Z} < " << pt2 <<" ) is not expected to do a lot." << endl;
507 +  //note that "15.0" here is the subleading jet cut.
508 +  
509    bool DoFit=false;
510    bool UseFit=false;
511 +  bool DoRooFit=false;
512    
513    
514    if(Contains(variable,"pfJetGood")) varname="R_{abs}";
515 <  TH1F *hdata  = allsamples.Draw("hdata",variable,1200,0,30, varname, "events",  cut,data,luminosity);
516 <  TH1F *hmc    = allsamples.Draw("hmc" , variable,1200,0,30, varname, "events", cut, mc,  luminosity);
515 > //  TH1F *hdata  = allsamples.Draw("hdata",variable,1200,0,30, varname, "events",  cut,data,luminosity);
516 > //  TH1F *hmc    = allsamples.Draw("hmc" , variable,1200,0,30, varname, "events", cut, mc,  luminosity);
517 >  TH1F *hdata, *hmc;
518 >  if(!DoRooFit) {
519 >    hdata = allsamples.Draw("hdata",variable,1200,ZbData::LowerHistoBoundary,ZbData::UpperHistoBoundary, varname, "events",  cut,data,luminosity); // making it more precise
520 >    hmc    = allsamples.Draw("hmc" , variable,1200,ZbData::LowerHistoBoundary,ZbData::UpperHistoBoundary, varname, "events", cut, mc,  luminosity); // making it more precise
521 >  } else {
522 >    hdata  = allsamples.Draw("hdata",variable,7200,0.0,6.0, varname, "events",  cut,data,luminosity); // making it more precise
523 >    hmc    = allsamples.Draw("hmc" , variable,7200,0.0,6.0, varname, "events", cut, mc,  luminosity); // making it more precise
524 >  }
525 >  
526    hmc->SetLineColor(kBlue);
527    
528    float a=hdata->GetMean();
# Line 355 | Line 545 | Value get_Zb_data_over_mc(string Contain
545      gaus_mc->SetParameter(2,hmc->GetRMS());
546      hmc->Fit(gaus_mc);
547    }
548 <  
548 >  if(DoRooFit) {
549 >    
550 >      if(hdata->Integral()>10) {
551 >        make_RooFit(ContainerName, hdata, a, da, saveas+"_data");
552 >        make_RooFit(ContainerName, hmc, b, db, saveas+"_mc");
553 >        factor = a / b;
554 >        error = TMath::Sqrt(  (1/(b*b)) * (da*da) + ((a*a)/(b*b))*db*db);
555 >      } else {
556 >        factor=NAN;//set to nan
557 >        error=NAN;//set to nan
558 >      }
559 >  }
560 >
561    cn->cd();
562 <  hdata->GetYaxis()->SetTitle("events / 0.1");
562 >  hdata->Rebin(int(0.05*1200/(ZbData::UpperHistoBoundary-ZbData::LowerHistoBoundary)));
563 >  hmc->Rebin(int(0.05*1200/(ZbData::UpperHistoBoundary-ZbData::LowerHistoBoundary)));
564 >  hdata->GetYaxis()->SetTitle("events / 0.05");
565    hdata->SetMarkerColor(kRed);
566    if(DoFit) {
567      gaus_mc->SetLineColor(kBlue);
# Line 367 | Line 571 | Value get_Zb_data_over_mc(string Contain
571    hdata->GetXaxis()->SetRangeUser(0,2);
572    hdata->GetYaxis()->SetTitleOffset(1.3);
573    hdata->Draw("e1");
574 +  
575 +  hmc->Scale(hdata->Integral()/hmc->Integral());
576 +  if(hmc->GetMaximum()>hdata->GetMaximum()) hdata->SetMaximum((hmc->GetMaximum()+sqrt((double)hmc->GetMaximum()))*1.2);
577 +  
578    hmc->Draw("histo,same");
579    hdata->Draw("e1,same");
580    if(DoFit) {
# Line 405 | Line 613 | Value get_Zb_data_over_mc(string Contain
613    CompleteSave(cn,ContainerName+"/ResponsePlots/"+saveas);
614    
615    res.MPF_Rabs_location.push_back(ZbPointResult(variable,pt1,pt2,a1,a2,Value(a,da),Value(b,db)));
616 <  cout << "Data;" << a << ";"<<da<<";MC;"<<b<<";"<<db<<endl;
616 >  dout << "Data;" << a << ";"<<da<<";MC;"<<b<<";"<<db<<endl;
617    delete cn;
618    delete hdata;
619    delete hmc;
# Line 414 | Line 622 | Value get_Zb_data_over_mc(string Contain
622   }
623    
624   int nFakes=1;
417 – ZbResultContainer Fake_new_data_mc_agreement_2d(string AlphaVariable, string ContainerName, TCut BTagCut, TCut BTagWgt) {
418 –  write_warning(__FUNCTION__,"You are using the fake way to extract the correction ... ");
419 –  nFakes++;
420 –  ZbResultContainer results(ContainerName);
421 –  
422 –  
423 –  results.MPF_Result_FaceValue=Value(1+nFakes*0.02,0.03);
424 –  results.Rabs_Result_FaceValue=Value(2,0.02);
425 –  results.MPF_Result_Extrapolated=Value(3,0.03);
426 –  results.Rabs_Result_Extrapolated=Value(4,0.04);
427 –  cout << "The fake result for " << ContainerName << " for the MPF face value is " << results.MPF_Result_FaceValue << endl;
625  
626 <  return results;
626 > float PurityInclusive=0.05;
627 > float PurityLoose=0.20;
628 > float PurityMedium=0.54;
629 > float PurityTight=0.82;
630 >
631 >
632 > void GetPtPointPosition(float AvPos[], const int nptcuts,string AlphaVariable, TCut basecut, float FaceValueToBeTakenAt, float ptcuts[]) {
633 >  for(int i=0;i<nptcuts-1;i++) {
634 >    stringstream specialcut;
635 >    specialcut << "((pt>30&& pt< 1000) && (" << AlphaVariable << "<" << FaceValueToBeTakenAt << "))";
636 >    TH1F *hdata  = allsamples.Draw("hada", "pt",100,ptcuts[i],ptcuts[i+1], "MPF", "events", TCut(basecut && specialcut.str().c_str()), data,  luminosity);
637 >    AvPos[i]=hdata->GetMean();
638 >    dout << AvPos[i] << " : ";
639 >    delete hdata;
640 >  }
641 >  dout << endl;
642   }
643 +    
644  
645 + void DeterminePurity(TCut basecut, int nptcuts, float ptcuts[], string AlphaVariable, string ContainerName, float FaceValueToBeTakenAt) {
646 +  
647 +  bool ComputePurity=false;
648 +  
649 +  if(ContainerName=="inclusive") ComputePurity=true;  
650 +  if(ContainerName=="loose") ComputePurity=true;
651 +  if(ContainerName=="medium") ComputePurity=true;
652 +  if(ContainerName=="reference") ComputePurity=true;
653 +  
654 +  if(!ComputePurity) return;
655 +  
656 +  stringstream specialcut;
657 +  specialcut << "((pt>30&& pt< 1000) && (" << AlphaVariable << "<" << FaceValueToBeTakenAt << "))";
658 +  
659 +  THStack mcm    = allsamples.DrawStack("mcm" , "mpf",1200,0,30, "MPF", "events", TCut(basecut && specialcut.str().c_str()), mc,  luminosity);
660 +  int nhists = mcm.GetHists()->GetEntries();
661 +  
662 +  float total=0;
663 +  float Zb=0;
664 +  
665 +  for (int i = 0; i < nhists; i++) {
666 +    TH1* hist = static_cast<TH1*>(mcm.GetHists()->At(i));
667 +    string name=hist->GetName();
668 +    if(Contains(name,"t_bar_t")) total+=hist->Integral();
669 +    if(Contains(name,"W_Jets")) total+=hist->Integral();
670 +    if(Contains(name,"Single_top")) total+=hist->Integral();
671 +    if(Contains(name,"Dibosons")) total+=hist->Integral();
672 +    if(Contains(name,"Z_l_")) total+=hist->Integral();
673 +    if(Contains(name,"Z_c_")) total+=hist->Integral();
674 +    if(Contains(name,"Z_b_")) {
675 +      total+=hist->Integral();
676 +      Zb+=hist->Integral();
677 +    }
678 +  }
679 +  
680 +  float purity=Zb/total;
681 +  dout << "Determined a purity of " << 100*purity << " % for the container \"" << ContainerName << "\"" << endl;
682 +  
683 +  if(ContainerName=="inclusive") PurityInclusive=purity;
684 +  if(ContainerName=="loose") PurityLoose=purity;
685 +  if(ContainerName=="medium") PurityMedium=purity;
686 +  if(ContainerName=="reference") PurityTight=purity;
687 +  
688 +  DeleteStack(mcm);
689 + }
690 +
691 + string NiceAlphaPrint(float alpha) {
692 +  //get in : 0.3, return 0p3 ...
693 +  stringstream res;
694 +  res << int(alpha) << "p" << int(100*alpha);
695 +  return res.str();
696 + }
697  
698   ZbResultContainer new_data_mc_agreement_2d(bool gofast, string AlphaVariable, string ContainerName, TCut BTagCut, TCut BTagWgt, float FaceValueToBeTakenAt) {
699    
# Line 438 | Line 703 | ZbResultContainer new_data_mc_agreement_
703    cutWeight=STDWEIGHT*BTagWgt;
704    LeadingB=BTagCut;
705    
706 +  
707    gStyle->SetOptFit(0);
708    TCut basecut(ZplusBsel&&LeadingB&&EtaB&&PhiZcut);
709    const int nalphacuts=6;
710    float alphacuts[nalphacuts] = {0.1,0.15,0.2,0.25,0.3,0.35};
711    
712 +  
713    const int nptcuts=5;
714    float ptcuts[nptcuts]={30,50,75,125,1000};
715 +  float ptPointPosition[nptcuts];
716 +  GetPtPointPosition(ptPointPosition,nptcuts,AlphaVariable,basecut,FaceValueToBeTakenAt,ptcuts);
717    
718    
719 +  DeterminePurity(basecut, nptcuts,ptcuts, AlphaVariable, ContainerName, FaceValueToBeTakenAt);
720    
721    float MPF_Results[nalphacuts][nptcuts];
722    float MPF_Errors[nalphacuts][nptcuts];
# Line 467 | Line 737 | ZbResultContainer new_data_mc_agreement_
737          alow=0;
738          ahigh=alphacuts[ia];
739        }
740 <      cout << specialcut.str() << endl;
741 <      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(alphacuts[ia]), ptcuts[ipt], ptcuts[ipt+1],alow,ahigh);
742 <      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(alphacuts[ia]), ptcuts[ipt], ptcuts[ipt+1],alow,ahigh);
740 >      dout << specialcut.str() << endl;
741 >      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);
742 >      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);
743        
744        MPF_Results[ia][ipt]=MPF_data_over_mc.getValue();
745        MPF_Errors[ia][ipt]=MPF_data_over_mc.getError();
# Line 477 | Line 747 | ZbResultContainer new_data_mc_agreement_
747        RABS_Results[ia][ipt]=RABS_data_over_mc.getValue();
748        RABS_Errors[ia][ipt]=RABS_data_over_mc.getError();
749        
750 <      cout << "alpha cut at " << alphacuts[ia+1] << " and pt range " << ptcuts[ipt] << " , " << ptcuts[ipt+1] << " \t : \t" << specialcut.str() << " \t \t --> " << MPF_data_over_mc << " (MPF) " << RABS_data_over_mc << " (RABS)" << endl;
750 >      dout << "alpha in [" << alow << " , " << ahigh << " and pt range [" << ptcuts[ipt] << " , " << ptcuts[ipt+1] << "] \t : \t" << specialcut.str() << " \t \t --> " << MPF_data_over_mc << " (MPF) " << RABS_data_over_mc << " (RABS)" << endl;
751      }
752    }
753    
# Line 496 | Line 766 | ZbResultContainer new_data_mc_agreement_
766      
767      stringstream filename;
768      filename << "Extrapolation/Zplusb_data_over_mc___" << ptcuts[ipt] << "_to_" << ptcuts[ipt+1];
769 <    cout << "Going to create histo for " << filename.str() << endl;
769 >    dout << "Going to create histo for " << filename.str() << endl;
770      TGraphErrors *mpf_gr =new TGraphErrors();
771      TGraphErrors *rabs_gr =new TGraphErrors();
772 +    mpf_gr->SetName("MPF_Graph");
773 +    rabs_gr->SetName("Rabs_Graph");
774      
775      int rabs_counter=0;
776      int mpf_counter=0;
777      
778      for(int ia=0;ia<nalphacuts && !gofast;ia++) {
779 <      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;
779 >      dout << "Setting a point for an alpha cut at " << alphacuts[ia] << " with value " << MPF_Results[ia][ipt] << " (MPF) and " << RABS_Results[ia][ipt] << " (RABS)" << endl;
780        if(MPF_Results[ia][ipt]==MPF_Results[ia][ipt] && MPF_Errors[ia][ipt]==MPF_Errors[ia][ipt]) { // remove nan's
781         mpf_gr->SetPoint(mpf_counter,alphacuts[ia],MPF_Results[ia][ipt]);
782         mpf_gr->SetPointError(mpf_counter,0,MPF_Errors[ia][ipt]);
# Line 512 | Line 784 | ZbResultContainer new_data_mc_agreement_
784        } else {
785          dout << "Detected BS for MPF method (bin ignored) : " << endl;
786          dout << "       alpha: " << alphacuts[ia] << " , pt range: " << ptcuts[ipt] << " < pt < " << ptcuts[ipt+1] << endl;
787 <        cout << "       value (MPF) : " << MPF_Results[ia][ipt] << " +/- " << MPF_Errors[ia][ipt] << endl;
787 >        dout << "       value (MPF) : " << MPF_Results[ia][ipt] << " +/- " << MPF_Errors[ia][ipt] << endl;
788        }
789        if(RABS_Results[ia][ipt]==RABS_Results[ia][ipt] && RABS_Errors[ia][ipt]==RABS_Errors[ia][ipt]) {
790         rabs_gr->SetPoint(rabs_counter,alphacuts[ia],RABS_Results[ia][ipt]);
# Line 521 | Line 793 | ZbResultContainer new_data_mc_agreement_
793        } else {
794          dout << "Detected BS for RABS method (bin ignored) : " << endl;
795          dout << "       alpha: " << alphacuts[ia] << " , pt range: " << ptcuts[ipt] << " < pt < " << ptcuts[ipt+1] << endl;
796 <        cout << "       value (RABS) : " << RABS_Results[ia][ipt] << " +/- " << RABS_Errors[ia][ipt] << endl;
796 >        dout << "       value (RABS) : " << RABS_Results[ia][ipt] << " +/- " << RABS_Errors[ia][ipt] << endl;
797        }
798      }
799      
# Line 529 | Line 801 | ZbResultContainer new_data_mc_agreement_
801      specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (" << AlphaVariable << "<" << FaceValueToBeTakenAt << "))";
802      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);
803      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);
804 <    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;
804 >    dout << "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;
805      MPF_FaceValueAtPoint3->SetPoint(ipt,(ptcuts[ipt]+ptcuts[ipt+1])/2,MPF_data_over_mc.getValue());
806      MPF_FaceValueAtPoint3->SetPointError(ipt,0,MPF_data_over_mc.getError());
807      RABS_FaceValueAtPoint3->SetPoint(ipt,(ptcuts[ipt]+ptcuts[ipt+1])/2,RABS_data_over_mc.getValue());
808      RABS_FaceValueAtPoint3->SetPointError(ipt,0,RABS_data_over_mc.getError());
809      
810      
811 <    
811 >    TH1D *mpf_band, *rabs_band;
812      
813      if(!gofast) {
814        can->cd();
815        mpf_gr->SetMarkerStyle(21);
816        mpf_gr->SetMarkerColor(kBlue);
817        mpf_gr->Draw("AP");
818 <      mpf_gr->Fit("pol1","");
818 >      
819 >      TF1 *mpf_pol = new TF1("mpf_pol","[0]+[1]*x",0.0,0.4);
820 >      
821 >      mpf_gr->Fit(mpf_pol); // better error estimation using MINOS (E), quiet (Q) and using loglikelihood method
822 >      mpf_band = getBand(mpf_pol,"MPF_band");
823        mpf_gr->GetXaxis()->SetTitle("#alpha");
824        mpf_gr->GetXaxis()->CenterTitle();
825        mpf_gr->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
826        mpf_gr->GetYaxis()->CenterTitle();
827 <      TF1 *mpf_pol=(TF1*)mpf_gr->GetFunction("pol1");
827 >      mpf_band->SetMarkerSize(0);
828        mpf_pol->SetLineWidth(0);
829 +      mpf_pol->SetLineColor(mpf_band->GetFillColor());
830        TF1 *mpf_pol_complete = new TF1("mpf_pol_complete","[0]+[1]*x");
831        mpf_pol_complete->SetParameters(mpf_pol->GetParameters());
832 +      
833        mpf_pol_complete->SetLineWidth(1);
834        mpf_pol_complete->SetLineColor(kBlue);
835        
# Line 568 | Line 846 | ZbResultContainer new_data_mc_agreement_
846        histo->GetYaxis()->SetRangeUser(min,max);
847        mpf_gr->GetYaxis()->SetRangeUser(min,max);
848        
849 <      TPolyLine *mpf_fit_uncert = GetFitUncertaintyShape(mpf_gr, "pol1", min, max,0.0,0.4);
849 >      //TPolyLine *mpf_fit_uncert = GetFitUncertaintyShape(mpf_gr, "pol1", min, max,0.0,0.4);
850 >      
851        mpf_pol->SetLineColor(TColor::GetColor("#0101DF"));
852 <      mpf_fit_uncert->SetFillColor(TColor::GetColor("#5353FC"));
852 >      //mpf_fit_uncert->SetFillColor(TColor::GetColor("#A2A2FA"));
853 >      mpf_band->SetFillColor(TColor::GetColor("#A2A2FA"));
854        histo->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
855        histo->Draw();
856 <      mpf_fit_uncert->Draw("F");
857 <      mpf_fit_uncert->Draw("");
856 >      //mpf_fit_uncert->Draw("F");
857 >      //mpf_fit_uncert->Draw("");
858 >      mpf_gr->Draw("P");
859 >      mpf_band->Draw("C E3 SAME");
860        mpf_gr->Draw("P");
861        mpf_pol_complete->Draw("same");
862        
863        float mpf_a=mpf_pol->GetParameter(0);
864        float mpf_b=mpf_pol->GetParameter(1);
865 <      MPF_FinalGraph->SetPoint(ipt,0.5*(ptcuts[ipt]+ptcuts[ipt+1]),mpf_a);
866 <      MPF_FinalGraph->SetPointError(ipt,0,mpf_pol->GetParError(0));
865 >      MPF_FinalGraph->SetPoint(ipt,0.5*(ptcuts[ipt]+ptcuts[ipt+1]),mpf_pol->GetParameter(0));
866 >      MPF_FinalGraph->SetPointError(ipt,0,mpf_band->GetBinError(mpf_band->FindBin(0)));
867        
868        MPF_ExtraPolatedResults[ipt]=mpf_a;
869        
# Line 598 | Line 880 | ZbResultContainer new_data_mc_agreement_
880        filename << "__MPF";
881        DrawPrelim();
882        CompleteSave(can,ContainerName+"/"+filename.str());
883 <      cout << "MPF : " << mpf_a << " + " << mpf_b << " * alpha " << endl;
883 >      dout << "MPF : " << mpf_a << " + " << mpf_b << " * alpha " << endl;
884        
885        filename.str("");
886        
887 +      TF1 *rabs_pol = new TF1("rabs_pol","[0]+[1]*x",0.0,0.4);
888 +      rabs_gr->Fit(rabs_pol,"E");
889 +      rabs_band = getBand(rabs_pol,"Rabs_band");
890 +      rabs_pol->SetLineColor(rabs_band->GetFillColor());
891        rabs_gr->SetMarkerStyle(21);
892        rabs_gr->SetMarkerColor(kRed);
893        rabs_gr->Draw("AP");
894 <      rabs_gr->Fit("pol1","");
894 >      //rabs_gr->Fit("pol1","E");
895 >      rabs_band->SetMarkerSize(0);
896 >      
897 >      
898        rabs_gr->GetXaxis()->SetTitle("#alpha");
899        rabs_gr->GetXaxis()->CenterTitle();
900        rabs_gr->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
901        rabs_gr->GetYaxis()->CenterTitle();
902 <      TF1 *rabs_pol=(TF1*)rabs_gr->GetFunction("pol1");
902 >      //TF1 *rabs_pol=(TF1*)rabs_gr->GetFunction("pol1");
903        rabs_pol->SetLineWidth(0);
904        TF1 *rabs_pol_complete = new TF1("rabs_pol_complete","[0]+[1]*x");
905        rabs_pol_complete->SetParameters(rabs_pol->GetParameters());
# Line 619 | Line 908 | ZbResultContainer new_data_mc_agreement_
908        FindMinMax(rabs_gr,min,max);
909        rabs_gr->GetYaxis()->SetRangeUser(min,max);
910        histo->GetYaxis()->SetRangeUser(min,max);
911 <      TPolyLine *rabs_fit_uncert = GetFitUncertaintyShape(rabs_gr, "pol1", min, max,0.0,0.4);
912 <      rabs_pol->SetLineColor(TColor::GetColor("#FA5858"));
913 <      rabs_pol->SetFillColor(TColor::GetColor("#FA5858"));
914 <      rabs_fit_uncert->SetLineColor(TColor::GetColor("#FA5858"));
915 <      rabs_fit_uncert->SetFillColor(TColor::GetColor("#FA5858"));
911 >      //TPolyLine *rabs_fit_uncert = GetFitUncertaintyShape(rabs_gr, "pol1", min, max,0.0,0.4);
912 >      rabs_pol->SetLineColor(TColor::GetColor("#F99C9C"));
913 >      rabs_pol->SetFillColor(TColor::GetColor("#F99C9C"));
914 > //      rabs_band->SetLineColor(TColor::GetColor("#F99C9C"));
915 >      rabs_band->SetFillColor(TColor::GetColor("#F99C9C"));
916 >      //rabs_fit_uncert->SetLineColor(TColor::GetColor("#F99C9C"));
917 >      //rabs_fit_uncert->SetFillColor(TColor::GetColor("#F99C9C"));
918        histo->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
919        histo->Draw();
920 <      rabs_fit_uncert->Draw("F");
921 <      rabs_fit_uncert->Draw("");
920 >      rabs_band->Draw("C E3 SAME");
921 >      //rabs_fit_uncert->Draw("F");
922 >      //rabs_fit_uncert->Draw("");
923        rabs_gr->Draw("P");
924        rabs_pol_complete->Draw("same");
925        
634 –      
926        float rabs_a=rabs_pol->GetParameter(0);
927        float rabs_b=rabs_pol->GetParameter(1);
928 +      
929        RABS_FinalGraph->SetPoint(ipt,0.5*(ptcuts[ipt]+ptcuts[ipt+1]),rabs_a);
930 <      RABS_FinalGraph->SetPointError(ipt,0,rabs_pol->GetParError(0));
930 >      RABS_FinalGraph->SetPointError(ipt,0,rabs_band->GetBinError(rabs_band->FindBin(0)));
931        
932 <      cout << "!+!+!+!+!+!+!+!+!+ Just added a point to the final plot : " << rabs_a << " +/- " << rabs_pol->GetParError(0) << endl;
932 >      dout << "!+!+!+!+!+!+!+!+!+ Just added a point to the final plot : " << rabs_a << " +/- " << rabs_band->GetBinError(rabs_band->FindBin(0)) << endl;
933        
934        stringstream rabs_info;
935        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() << "}";
# Line 651 | Line 943 | ZbResultContainer new_data_mc_agreement_
943        filename << filenamebkp << "__RABS";
944        DrawPrelim();
945        CompleteSave(can,ContainerName+"/"+filename.str());
946 <      cout << "RABS : " << rabs_a << " + " << rabs_b << " * alpha " << endl;
946 >      dout << "RABS : " << rabs_a << " + " << rabs_b << " * alpha " << endl;
947 >      
948 >      delete mpf_band;
949 >      delete rabs_band;
950      }
951    }
952    
# Line 787 | Line 1082 | ZbResultContainer new_data_mc_agreement_
1082    filename2 << "Extrapolation/FINAL_RESULT_RABS_FaceValp3";
1083    CompleteSave(can,ContainerName+"/"+filename2.str());
1084    
1085 <  cout << "FINAL RESULTS: " << endl;
1086 <  if(!gofast) cout << "MPF : " << MPFResult  << " +/- " << MPFResultError  << endl;
1087 <  if(!gofast) cout << "Rabs: " << RabsResult << " +/- " << RabsResultError << endl;
1088 <  cout << "Face value at " << FaceValueToBeTakenAt << " :" << FaceValResult << " +/- " << FaceValResultError << "  (MPF) " << endl;
1089 <  cout << "Face value at " << FaceValueToBeTakenAt << " :" << RABSFaceValResult << " +/- " << RABSFaceValResultError << "  (RABS) " << endl;
1085 >  dout << "FINAL RESULTS: " << endl;
1086 >  if(!gofast) dout << "MPF : " << MPFResult  << " +/- " << MPFResultError  << endl;
1087 >  if(!gofast) dout << "Rabs: " << RabsResult << " +/- " << RabsResultError << endl;
1088 >  dout << "Face value at " << FaceValueToBeTakenAt << " :" << FaceValResult << " +/- " << FaceValResultError << "  (MPF) " << endl;
1089 >  dout << "Face value at " << FaceValueToBeTakenAt << " :" << RABSFaceValResult << " +/- " << RABSFaceValResultError << "  (RABS) " << endl;
1090    
1091    delete can;
1092    
1093 +  
1094 +  
1095 +  for(unsigned int i=0;i<MPF_FinalGraph->GetN();i++) {
1096 +    double x,y,dy;
1097 +    MPF_FinalGraph->GetPoint(i,x,y);
1098 +    dy = MPF_FinalGraph->GetErrorY(i);
1099 +    results.MPF_Result_Extrapolated_Bins.push_back(Value(y,dy));
1100 +  }
1101 +  
1102 +  for(unsigned int i=0;i<MPF_FaceValueAtPoint3->GetN();i++) {
1103 +    double x,y,dy;
1104 +    MPF_FaceValueAtPoint3->GetPoint(i,x,y);
1105 +    dy = MPF_FaceValueAtPoint3->GetErrorY(i);
1106 +    results.MPF_Result_FaceValue_Bins.push_back(Value(y,dy));
1107 +  }
1108 +
1109 +  for(unsigned int i=0;i<RABS_FinalGraph->GetN();i++) {
1110 +    double x,y,dy;
1111 +    RABS_FinalGraph->GetPoint(i,x,y);
1112 +    dy = RABS_FinalGraph->GetErrorY(i);
1113 +    results.Rabs_Result_Extrapolated_Bins.push_back(Value(y,dy));
1114 +  }
1115 +  
1116 +  for(unsigned int i=0;i<RABS_FaceValueAtPoint3->GetN();i++) {
1117 +    double x,y,dy;
1118 +    RABS_FaceValueAtPoint3->GetPoint(i,x,y);
1119 +    dy = RABS_FaceValueAtPoint3->GetErrorY(i);
1120 +    results.Rabs_Result_FaceValue_Bins.push_back(Value(y,dy));
1121 +  }
1122 +  
1123 +  
1124    results.MPF_Result_FaceValue=Value(FaceValResult,FaceValResultError);
1125    results.Rabs_Result_FaceValue=Value(RABSFaceValResult,RABSFaceValResultError);
1126    if(!gofast) results.MPF_Result_Extrapolated=Value(MPFResult,MPFResultError);
1127    if(!gofast) results.Rabs_Result_Extrapolated=Value(RabsResult,RabsResultError);
1128 <  cout << results << endl;
1128 >  results.Print();
1129    return results;
1130   }
1131  
# Line 811 | Line 1137 | void DoPUStudy(string identifier) {
1137    TCanvas *can = new TCanvas("can","can");
1138    TH1F *malpha[8];
1139    TH1F *dalpha[8];
1140 <  cout << identifier << ";";
1141 <  cout << endl;
1140 >  dout << identifier << ";";
1141 >  dout << endl;
1142    for(int i=1;i<5;i++) {
1143      stringstream sSpecialCut;
1144      if(i<4) {
1145        sSpecialCut << "numVtx>" << numVtxcuts[i-1] << " && numVtx<" << numVtxcuts[i];
1146 <      cout << numVtxcuts[i-1] << "  < nVtx < " << numVtxcuts[i] << ";";
1146 >      dout << numVtxcuts[i-1] << "  < nVtx < " << numVtxcuts[i] << ";";
1147      } else {
1148        sSpecialCut << "numVtx>" << numVtxcuts[i-1];
1149 <      cout << numVtxcuts[i-1] << ";";
1149 >      dout << numVtxcuts[i-1] << ";";
1150      }
1151      
1152      sSpecialCut << " && " << identifier << "_alpha<0.3";
# Line 838 | Line 1164 | void DoPUStudy(string identifier) {
1164      float factor = a / b;
1165      float error = TMath::Sqrt(  (1/(b*b)) * (da*da) + ((a*a)/(b*b))*db*db);
1166    
1167 <    cout << malpha[i]->GetMean() << ";" << malpha[i]->GetMeanError() << ";" << dalpha[i]->GetMean() << ";" << dalpha[i]->GetMeanError() << ";" << malpha[i]->Integral()<<";"<<dalpha[i]->Integral() << ";"<<factor << ";" << error <<  endl;
1167 >    dout << malpha[i]->GetMean() << ";" << malpha[i]->GetMeanError() << ";" << dalpha[i]->GetMean() << ";" << dalpha[i]->GetMeanError() << ";" << malpha[i]->Integral()<<";"<<dalpha[i]->Integral() << ";"<<factor << ";" << error <<  endl;
1168      delete malpha[i];
1169      delete dalpha[i];
1170    }
# Line 919 | Line 1245 | void ScenarioComparison() {
1245      }
1246      
1247      
1248 +    gStyle->SetOptFit(0);
1249      gr[i]->SetTitle();
1250      gr[i]->GetXaxis()->SetTitle("N(Vtx)");
1251      gr[i]->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
# Line 1051 | Line 1378 | void ScenarioComparisonInclusive() {
1378   void compare_selection(string identifier) {
1379    bool recognized_scenario=false;
1380    
1381 <  cout << "Running with identifier " << identifier << endl;
1381 >  dout << "Running with identifier " << identifier << endl;
1382    if(identifier=="Zb30") {
1383        LeadingB=TCut ("Zb30_bTagProbCSVBP[0]>0.898");
1384        EtaB=TCut("abs(Zb30_pfJetGoodEta[0])<1.3");
# Line 1163 | Line 1490 | void compare_selection(string identifier
1490    
1491    assert(recognized_scenario);
1492    
1493 <  cout << "Cuts: " << endl;
1494 <  cout << "   " << (const char*) LeadingB << endl;
1495 <  cout << "   " << (const char*) EtaB << endl;
1496 <  cout << "   " << (const char*) PhiZcut << endl;
1493 >  dout << "Cuts: " << endl;
1494 >  dout << "   " << (const char*) LeadingB << endl;
1495 >  dout << "   " << (const char*) EtaB << endl;
1496 >  dout << "   " << (const char*) PhiZcut << endl;
1497    
1498 <  cout << endl << endl;
1498 >  dout << endl << endl;
1499  
1500   //  ScenarioComparison();
1501   //  ScenarioComparisonInclusive();
# Line 1223 | Line 1550 | void TEST_GetNumberEventsInsideOutsideAl
1550   //   GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==1"), "mm");
1551   //   GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2"), "ee/mm");
1552    
1553 <  cout << "BASE SELECTION" << endl;
1553 >  dout << "BASE SELECTION" << endl;
1554    GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==0"), "ee");
1555    GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==1"), "mm");
1556    GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2"), "ee/mm");
1557  
1558 <  cout << "BASE SELECTION: Z+b" << endl;
1558 >  dout << "BASE SELECTION: Z+b" << endl;
1559    GetNumberEventsInsideOutsideAlphaWindow(ZplusBsel&&TCut("id1==id2&&id1==0"), "ee");
1560    GetNumberEventsInsideOutsideAlphaWindow(ZplusBsel&&TCut("id1==id2&&id1==1"), "mm");
1561    GetNumberEventsInsideOutsideAlphaWindow(ZplusBsel&&TCut("id1==id2"), "ee/mm");
1562    
1563 <  cout << "Leading B" << endl;
1563 >  dout << "Leading B" << endl;
1564    GetNumberEventsInsideOutsideAlphaWindow(LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==0"), "ee");
1565    GetNumberEventsInsideOutsideAlphaWindow(LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==1"), "mm");
1566    GetNumberEventsInsideOutsideAlphaWindow(LeadingB&&ZplusBsel&&TCut("id1==id2"), "ee/mm");
1567    
1568 <  cout << "PhiZcut" << endl;
1568 >  dout << "PhiZcut" << endl;
1569    GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==0"), "ee");
1570    GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==1"), "mm");
1571    GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2"), "ee/mm");
1572    
1573 <  cout << "alpha cut" << endl;
1573 >  dout << "alpha cut" << endl;
1574    GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==0")&&TCut("ZbCHS3010_alpha<0.3"), "ee");
1575    GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==1")&&TCut("ZbCHS3010_alpha<0.3"), "mm");
1576    GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2")&&TCut("ZbCHS3010_alpha<0.3"), "ee/mm");
# Line 1254 | Line 1581 | void TEST_GetNumberEventsInsideOutsideAl
1581   }
1582    
1583  
1584 <
1584 > void MetSpectrumDepletionIllustration(int isdata) {
1585 >  
1586 >  TCanvas *can = new TCanvas("can","can");
1587 >  TCut basecut(ZplusBsel&&LeadingB&&EtaB&&PhiZcut);
1588 >  
1589 >  cutWeight=STDWEIGHT*TCut("ZbCHS3010_BTagWgtTUp*(1.0)*(ZbCHS3010_bTagProbCSVBP[0]>0.898)");
1590 >  TH1F *standard  = allsamples.Draw("standard","met[1]",80,0,400, "MET [GeV]", "events", basecut,isdata,luminosity);
1591 >  
1592 >  cutWeight=STDWEIGHT*TCut("ZbCHS3010_BTagWgtTUp*(1.0*softMuon)*(ZbCHS3010_bTagProbCSVBP[0]>0.898)");
1593 >  TH1F *Neutrinos  = allsamples.Draw("Neutrinos","met[1]",80,0,400, "MET [GeV]", "events", basecut,isdata,luminosity);
1594 >  
1595 >  cutWeight=STDWEIGHT*TCut("ZbCHS3010_BTagWgtTUp*(1.0*(1.0-softMuon))*(ZbCHS3010_bTagProbCSVBP[0]>0.898)");
1596 >  TH1F *NoNeutrinos  = allsamples.Draw("NoNeutrinos","met[1]",80,0,400, "MET [GeV]", "events", basecut,isdata,luminosity);
1597 >  
1598 >  standard->Scale(1.0/standard->Integral());
1599 >  Neutrinos->Scale(1.0/Neutrinos->Integral());
1600 >  NoNeutrinos->Scale(1.0/NoNeutrinos->Integral());
1601 >  
1602 >  Neutrinos->SetLineColor(kGreen);
1603 >  NoNeutrinos->SetLineColor(kRed);
1604 >  
1605 >  can->SetLogy(1);
1606 >  TLegend *leg = make_legend();
1607 >  leg->AddEntry(standard,"Standard sample","ep");
1608 >  leg->AddEntry(Neutrinos,"#nu enriched sample","l");
1609 >  leg->AddEntry(NoNeutrinos,"#nu depleted sample","l");
1610 >  leg->SetX1(0.9*leg->GetX1());
1611 >  standard->Draw("e1");
1612 >  Neutrinos->Draw("histo,same");
1613 >  NoNeutrinos->Draw("histo,same");
1614 >  if(isdata==data) leg->SetHeader("Data:");
1615 >  else leg->SetHeader("MC:");
1616 >  leg->Draw();
1617 >  DrawPrelim();
1618 >  
1619 >  if(isdata==data) CompleteSave(can,"NeutrinoMetSpectrum__Data");
1620 >  else CompleteSave(can,"NeutrinoMetSpectrum__MC");
1621 >  delete standard;
1622 >  delete Neutrinos;
1623 >  delete NoNeutrinos;
1624 >  delete can;
1625 >  CleanLegends();
1626 > }
1627 >    
1628 > string PrintVector(vector<float> Systematics) {
1629 >  stringstream fullanswer;
1630 >  fullanswer << " [ ";
1631 >  if(Systematics.size()>0) fullanswer << Systematics[0];
1632 >  for(unsigned int is=1;is<Systematics.size();is++) fullanswer << " , " << Systematics[is];
1633 >  fullanswer << " ] ";
1634 >  return fullanswer.str();
1635 > }
1636 >  
1637      
1638   void do_basic_ZB_analysis(bool DoCompleteAnalysis) {
1639    STDWEIGHT=TCut(cutWeight);
1640    cutWeight=TCut(STDWEIGHT*TCut("ZbCHS3010_BTagWgtT"));
1641 <  cout << "The lepton requirement is " << (const char*) leptoncut << endl;
1641 >  dout << "The lepton requirement is " << (const char*) leptoncut << endl;
1642    bool doquick=false;
1643    
1644     TCanvas *zbcanvas = new TCanvas("zbcanvas","zbcanvas");
1645    
1646 < //   ScenarioComparison();
1647 < //   ScenarioComparisonInclusive();
1648 <
1646 >  
1647 > //   MetSpectrumDepletionIllustration(data);
1648 > //   MetSpectrumDepletionIllustration(mc);
1649 >  
1650   //   compare_selection("ZbCHS3010_p5Clean");
1651   //   compare_selection("ZbCHS3010");
1652   //   compare_selection("Zb30");
# Line 1285 | Line 1665 | void do_basic_ZB_analysis(bool DoComplet
1665   //   compare_selection("Zb40");
1666  
1667    if(!doquick) {
1668 <    print_all_b_yields();
1669 <    draw_mpf_vars();
1670 <    DrawEvilCutFlow();
1668 > //    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.)
1669 > //    print_all_b_yields();
1670 > //     draw_mpf_vars();
1671 > //    DrawEvilCutFlow();
1672          
1673 <    draw_Zb_kin_vars();
1673 > //     draw_Zb_kin_vars();
1674          
1675     }
1676    
# Line 1306 | Line 1687 | void do_basic_ZB_analysis(bool DoComplet
1687      
1688      bool fast=true;bool slow=false;//don't change these
1689      
1690 <    ZbResultContainer inclusive = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","inclusive",TCut("mll>0"),TCut("1.0"),0.3);
1690 >    ZbResultContainer inclusive = new_data_mc_agreement_2d(slow,"ZbCHS3010_alpha","inclusive",TCut("mll>0"),TCut("1.0"),0.3);
1691      ZbResultContainer Reference = new_data_mc_agreement_2d(slow,"ZbCHS3010_alpha","reference",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtT"),0.3);
1692 +    
1693 +    
1694 +    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);
1695 +    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);
1696 +    
1697      ZbResultContainer effmistagUp = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","effmistagUp",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtTUp"),0.3);
1698      ZbResultContainer effmistagDown = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","effmistagDown",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtTDown"),0.3);
1699      ZbResultContainer medium = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","medium",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.679"),TCut("ZbCHS3010_BTagWgtM"),0.3);
# Line 1315 | Line 1701 | void do_basic_ZB_analysis(bool DoComplet
1701      ZbResultContainer AlphaUp = new_data_mc_agreement_2d(fast,"ZbCHS3010_alphaUp","AlphaUp",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtT"),0.3);
1702      ZbResultContainer AlphaDown = new_data_mc_agreement_2d(fast,"ZbCHS3010_alphaDown","AlphaDown",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtT"),0.3);
1703      
1704 <    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);
1705 <    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);
1706 <    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);
1707 <    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);
1704 >    ZbResultContainer L5inclusive = new_data_mc_agreement_2d(fast,"ZbCHS3010_alphaL5","L5inclusive",TCut("mll>0"),TCut("1.0"),0.3);
1705 >    ZbResultContainer L5Reference = new_data_mc_agreement_2d(fast,"ZbCHS3010_alphaL5","L5reference",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtT"),0.3);
1706 >    
1707 >    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);
1708 >    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);
1709 >    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);
1710 >    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);
1711 >    ZbResultContainer AlphaThresholdVariationDown10inc = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","AlphaThresholdVariationDown10inc",TCut("mll>0"),                              TCut("1.0"),0.3-0.1*0.3);
1712 >    ZbResultContainer AlphaThresholdVariationUp10inc =   new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","AlphaThresholdVariationUp10inc",  TCut("mll>0"),                              TCut("1.0"),0.3+0.1*0.3);
1713 >    ZbResultContainer AlphaThresholdVariationDown30inc = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","AlphaThresholdVariationDown30inc",TCut("mll>0"),                              TCut("1.0"),0.3-0.3*0.3);
1714 >    ZbResultContainer AlphaThresholdVariationUp30inc =   new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","AlphaThresholdVariationUp30inc",  TCut("mll>0"),                              TCut("1.0"),0.3+0.3*0.3);
1715 >    
1716 >    
1717 >    ZbResultContainer OnlyOneJet   = new_data_mc_agreement_2d(slow,"ZbCHS3010_alpha","OnlyOneJet",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898&&ZbCHS3010_pfJetGoodNum==1"),TCut("ZbCHS3010_BTagWgtT"),0.3);
1718 >    ZbResultContainer MultipleJets = new_data_mc_agreement_2d(slow,"ZbCHS3010_alpha","MultipleJets",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898&&ZbCHS3010_pfJetGoodNum>1"),TCut("ZbCHS3010_BTagWgtT"),0.3);
1719 >
1720 >    ZbResultContainer OnlyOneJetinc   = new_data_mc_agreement_2d(slow,"ZbCHS3010_alpha","OnlyOneJetinc",TCut("ZbCHS3010_pfJetGoodNum==1"),TCut("1.0"),0.3);
1721 >    ZbResultContainer MultipleJetsinc = new_data_mc_agreement_2d(slow,"ZbCHS3010_alpha","MultipleJetsinc",TCut("ZbCHS3010_pfJetGoodNum>1"),TCut("1.0"),0.3);
1722  
1723      //and now let's compute some systematics!
1724      //1 alpha variation
1725      float SysMPF_Alpha_down = Reference.MPF_Result_FaceValue.getValue()-AlphaDown.MPF_Result_FaceValue.getValue();
1726      float SysMPF_Alpha_up = Reference.MPF_Result_FaceValue.getValue()-AlphaUp.MPF_Result_FaceValue.getValue();
1727      
1728 +    vector<float> SysMPF_Alpha_down_Bins;
1729 +    vector<float> SysMPF_Alpha_up_Bins;
1730 +    for(unsigned int iFV=0;iFV<Reference.MPF_Result_FaceValue_Bins.size();iFV++) SysMPF_Alpha_down_Bins.push_back(Reference.MPF_Result_FaceValue_Bins[iFV].getValue()-AlphaDown.MPF_Result_FaceValue_Bins[iFV].getValue());
1731 +    for(unsigned int iFV=0;iFV<Reference.MPF_Result_FaceValue_Bins.size();iFV++) SysMPF_Alpha_up_Bins.push_back(Reference.MPF_Result_FaceValue_Bins[iFV].getValue()-AlphaUp.MPF_Result_FaceValue_Bins[iFV].getValue());
1732 +    
1733 +    
1734      //2 btagging efficiency/mistag correction
1735      float SysMPF_EFF_down = Reference.MPF_Result_FaceValue.getValue()-effmistagDown.MPF_Result_FaceValue.getValue();
1736      float SysMPF_EFF_up = Reference.MPF_Result_FaceValue.getValue()-effmistagUp.MPF_Result_FaceValue.getValue();
1737      
1738 <    //3 purity
1738 >    vector<float> SysMPF_EFF_down_Bins;
1739 >    vector<float> SysMPF_EFF_up_Bins;
1740 >    for(unsigned int iFV=0;iFV<Reference.MPF_Result_FaceValue_Bins.size();iFV++) SysMPF_EFF_down_Bins.push_back(Reference.MPF_Result_FaceValue_Bins[iFV].getValue()-effmistagDown.MPF_Result_FaceValue_Bins[iFV].getValue());
1741 >    for(unsigned int iFV=0;iFV<Reference.MPF_Result_FaceValue_Bins.size();iFV++) SysMPF_EFF_up_Bins.push_back(Reference.MPF_Result_FaceValue_Bins[iFV].getValue()-effmistagUp.MPF_Result_FaceValue_Bins[iFV].getValue());
1742 >    
1743 >    //3 Presence of soft muons (->neutrino question)
1744 >    float SysMPF_Neutrino_up = Reference.MPF_Result_FaceValue.getValue()-NeutrinoQ.MPF_Result_FaceValue.getValue();
1745 >    float SysMPF_Neutrino_down = Reference.MPF_Result_FaceValue.getValue()-ANeutrinoQ.MPF_Result_FaceValue.getValue();
1746 >    
1747 >    vector<float> SysMPF_Neutrino_up_Bins;
1748 >    vector<float> SysMPF_Neutrino_down_Bins;
1749 >    for(unsigned int iFV=0;iFV<Reference.MPF_Result_FaceValue_Bins.size();iFV++) SysMPF_Neutrino_up_Bins.push_back(Reference.MPF_Result_FaceValue_Bins[iFV].getValue()-NeutrinoQ.MPF_Result_FaceValue_Bins[iFV].getValue());
1750 >    for(unsigned int iFV=0;iFV<Reference.MPF_Result_FaceValue_Bins.size();iFV++) SysMPF_Neutrino_down_Bins.push_back(Reference.MPF_Result_FaceValue_Bins[iFV].getValue()-ANeutrinoQ.MPF_Result_FaceValue_Bins[iFV].getValue());
1751 >    
1752 >    
1753 >    //4 purity
1754      zbcanvas->cd();
1755      gStyle->SetOptFit(0);
1756      write_info(__FUNCTION__,"Don't know the exact purity at the moment, still needs to be determined!");
# Line 1340 | Line 1761 | void do_basic_ZB_analysis(bool DoComplet
1761      purityhisto->GetYaxis()->CenterTitle();
1762      
1763      TGraphErrors *gSysMPF_Purity = new TGraphErrors();
1764 <    gSysMPF_Purity->SetPoint(0,0.05,inclusive.MPF_Result_FaceValue.getValue());
1764 >    
1765 >    gSysMPF_Purity->SetPoint(0,PurityInclusive,inclusive.MPF_Result_FaceValue.getValue());
1766      gSysMPF_Purity->SetPointError(0,0.0,inclusive.MPF_Result_FaceValue.getError());
1767      
1768 <    gSysMPF_Purity->SetPoint(1,0.20,loose.MPF_Result_FaceValue.getValue());//x value?
1768 >    gSysMPF_Purity->SetPoint(1,PurityLoose,loose.MPF_Result_FaceValue.getValue());//x value?
1769      gSysMPF_Purity->SetPointError(1,0.0,loose.MPF_Result_FaceValue.getError());
1770 <    gSysMPF_Purity->SetPoint(2,0.54,medium.MPF_Result_FaceValue.getValue());//x value?
1770 >    gSysMPF_Purity->SetPoint(2,PurityMedium,medium.MPF_Result_FaceValue.getValue());//x value?
1771      gSysMPF_Purity->SetPointError(2,0.0,medium.MPF_Result_FaceValue.getError());
1772 <    gSysMPF_Purity->SetPoint(3,0.82,Reference.MPF_Result_FaceValue.getValue());//x value?
1772 >    gSysMPF_Purity->SetPoint(3,PurityTight,Reference.MPF_Result_FaceValue.getValue());//x value?
1773      gSysMPF_Purity->SetPointError(3,0.0,Reference.MPF_Result_FaceValue.getError());
1774      
1775      
# Line 1390 | Line 1812 | void do_basic_ZB_analysis(bool DoComplet
1812      
1813      CompleteSave(zbcanvas,"Systematics/Purity");
1814      
1815 +    dout << "Bin-by-bin systematics : " << endl;
1816 +    dout << "    Oversmearing : " << endl;
1817 +    dout << "      down: " << PrintVector(SysMPF_Alpha_down_Bins) << endl;
1818 +    dout << "      up  : " << PrintVector(SysMPF_Alpha_up_Bins) << endl;
1819 +    dout << "    eff/mistag variation: " << endl;
1820 +    dout << "      down: " << PrintVector(SysMPF_EFF_down_Bins) << endl;
1821 +    dout << "      up  : " << PrintVector(SysMPF_EFF_up_Bins) << endl;
1822 +    dout << "    neutrinos: " << endl;
1823 +    dout << "      down: " << PrintVector(SysMPF_Neutrino_down_Bins) << endl;
1824 +    dout << "      up  : " << PrintVector(SysMPF_Neutrino_up_Bins) << endl;
1825      
1826  
1827      float SysMPF_Purity = abs(ExtrapolatedResult-Reference.MPF_Result_FaceValue.getValue());
1828      
1829 <    float sys_alpha  = abs(SysMPF_Alpha_down)>abs(SysMPF_Alpha_up)?abs(SysMPF_Alpha_down):abs(SysMPF_Alpha_up);
1830 <    float sys_eff    = abs(SysMPF_EFF_down)>abs(SysMPF_EFF_up)?abs(SysMPF_EFF_down):abs(SysMPF_EFF_up);
1831 <    float sys_purity = abs(SysMPF_Purity);
1832 <    float sys        = sqrt(sys_alpha*sys_alpha+sys_eff*sys_eff+sys_purity*sys_purity);
1829 >    float sys_alpha    = abs(SysMPF_Alpha_down)>abs(SysMPF_Alpha_up)?abs(SysMPF_Alpha_down):abs(SysMPF_Alpha_up);
1830 >    float sys_eff      = abs(SysMPF_EFF_down)>abs(SysMPF_EFF_up)?abs(SysMPF_EFF_down):abs(SysMPF_EFF_up);
1831 >    float sys_purity   = abs(SysMPF_Purity);
1832 >    float sys_neutrino = abs(SysMPF_Neutrino_down)>abs(SysMPF_Neutrino_up)?abs(SysMPF_Neutrino_down):abs(SysMPF_Neutrino_up);
1833 >    float sys          = sqrt(sys_alpha*sys_alpha+sys_eff*sys_eff+sys_purity*sys_purity+sys_neutrino*sys_neutrino);
1834      
1835      dout << "MPF METHOD:" << endl;
1836      dout << "   Systematics :    down    up     (selected)" << endl;
1837 <    dout << "      Alpha variation      : " << SysMPF_Alpha_down << "   ,  " << SysMPF_Alpha_up << "   ( " << sys_alpha << ") " << endl;
1837 >    dout << "      Oversmearing         : " << SysMPF_Alpha_down << "   ,  " << SysMPF_Alpha_up << "   ( " << sys_alpha << ") " << endl;
1838      dout << "      eff/mistag variation : " << SysMPF_EFF_down << "   ,  " << SysMPF_EFF_up << "   ( " << sys_eff << ") " << endl;
1839      dout << "      purity               : " << SysMPF_Purity << "   ( " << sys_purity << " )" << endl;
1840 <    dout << "      total                : " << sys << endl;
1840 >    dout << "      neutrino             : " << SysMPF_Neutrino_down << "   ,  " << SysMPF_Neutrino_up << "   ( " << sys_neutrino << ") " << endl;
1841 >    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;
1842      dout << endl << endl;
1843      dout << "   FINAL RESULTS : " << endl;
1844      dout << "     C_{abs}^{b}   = " << Reference.MPF_Result_FaceValue << " (stat) +/- " << sys << " (sys) " << endl;
# Line 1414 | Line 1848 | void do_basic_ZB_analysis(bool DoComplet
1848      sysfinal *= sqrt(sys*sys/(Reference.MPF_Result_FaceValue.getValue()*Reference.MPF_Result_FaceValue.getValue()));
1849      //yes, this could be simplified but it would look like a bug.
1850      
1851 +
1852 +    
1853 +    
1854 +    dout << "Note: Varying the alpha cut (->face value) we get the following results: " << endl;
1855 +    dout << "-30%: " << AlphaThresholdVariationDown30.MPF_Result_FaceValue/AlphaThresholdVariationDown30inc.MPF_Result_FaceValue << " (stat)" << endl;
1856 +    dout << "-10%: " << AlphaThresholdVariationDown10.MPF_Result_FaceValue/AlphaThresholdVariationDown10inc.MPF_Result_FaceValue << " (stat)" << endl;
1857 +    dout << "  0%: " << Reference.MPF_Result_FaceValue/inclusive.MPF_Result_FaceValue << " (stat) +/- " << endl;
1858 +    dout << "+10%: " << AlphaThresholdVariationUp10.MPF_Result_FaceValue/AlphaThresholdVariationUp10inc.MPF_Result_FaceValue << " (stat)" << endl;
1859 +    dout << "+30%: " << AlphaThresholdVariationUp30.MPF_Result_FaceValue/AlphaThresholdVariationUp30inc.MPF_Result_FaceValue << " (stat)" << endl;
1860 +    
1861 +    dout << "          XCheck for Rabs, varying the alpha cut (->face value) we get the following results: " << endl;
1862 +    dout << "             -30%: " << AlphaThresholdVariationDown30.Rabs_Result_FaceValue/AlphaThresholdVariationDown30inc.Rabs_Result_FaceValue << " (stat)" << endl;
1863 +    dout << "             -10%: " << AlphaThresholdVariationDown10.Rabs_Result_FaceValue/AlphaThresholdVariationDown10inc.Rabs_Result_FaceValue << " (stat)" << endl;
1864 +    dout << "               0%: " << Reference.Rabs_Result_FaceValue/inclusive.Rabs_Result_FaceValue << " (stat) +/- " << endl;
1865 +    dout << "             +10%: " << AlphaThresholdVariationUp10.Rabs_Result_FaceValue/AlphaThresholdVariationUp10inc.Rabs_Result_FaceValue << " (stat)" << endl;
1866 +    dout << "             +30%: " << AlphaThresholdVariationUp30.Rabs_Result_FaceValue/AlphaThresholdVariationUp30inc.Rabs_Result_FaceValue << " (stat)" << endl;
1867 +
1868 +    
1869 +    dout << "  Also checked what happens when you only use events with one jet and only with more than one jets: " << endl;
1870 +    dout << "      Only 1      : " << OnlyOneJet.MPF_Result_FaceValue/OnlyOneJetinc.MPF_Result_FaceValue << " (stat) " << endl;
1871 +    dout << "      More than 1 : " << MultipleJets.MPF_Result_FaceValue/MultipleJetsinc.MPF_Result_FaceValue << " (stat) " << endl;
1872 +    
1873 +    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());
1874 +    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());
1875 +    float sys_alphavar = sys_alphavar_up>sys_alphavar_down?sys_alphavar_up:sys_alphavar_down;
1876 +    
1877 +    dout << " -> Additional systematic (from 30%) : " << sys_alphavar << endl;
1878 +    sysfinal=sqrt(sysfinal*sysfinal+sys_alphavar*sys_alphavar);
1879 +    dout << "Final systematic is : " << sysfinal << endl;
1880 +    
1881 +    dout << "   ******************************************** <  FINAL RESULT > ******************************************** " << endl;
1882      dout << "     C_{corr}      = " << Reference.MPF_Result_FaceValue/inclusive.MPF_Result_FaceValue << " (stat) +/- " << sysfinal << " (sys) " << endl;
1883 +    dout << "   ******************************************** < /FINAL RESULT > ******************************************** " << endl;
1884 +    
1885 +    dout << "Note that if L5 corrections are applied, we get " << L5Reference.MPF_Result_FaceValue/L5inclusive.MPF_Result_FaceValue << " (stat) " << endl;
1886      
1887      LeadingB=TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898");
1888      
# Line 1465 | Line 1933 | void SpecialCutFlow(string identifier) {
1933    
1934    
1935    TH1F *data  = allsamples.Draw("data",    MegaCut.str(),11,-0.5,10.5, "", "events", basecut,data,luminosity);
1936 <   THStack mcm = allsamples.DrawStack("mc", MegaCut.str(),11,-0.5,10.5, "", "events", basecut,mc,luminosity);
1936 >  THStack mcm = allsamples.DrawStack("mc", MegaCut.str(),11,-0.5,10.5, "", "events", basecut,mc,luminosity);
1937  
1938    float runningsum=0;
1939    for(int i=data->GetNbinsX();i>0;i--) {
# Line 1534 | Line 2002 | void SpecialCutFlow(string identifier) {
2002    
2003    delete crap;
2004    delete data;
2005 <  
2005 >  CleanLegends();
2006 >  DeleteStack(mcm);
2007 >
2008    
2009   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines