[ViewVC] Diff of: cvsroot/UserCode/cbrown/Development/Plotting/Modules/ZbTools.C

Comparing UserCode/cbrown/Development/Plotting/Modules/ZbTools.C (file contents):
Revision 1.3 by buchmann, Thu Jul 12 12:53:41 2012 UTC vs.
Revision 1.28 by buchmann, Fri Apr 26 06:52:09 2013 UTC

+#include <iostream>
+#include <vector>
+#include <sys/stat.h>
-+
+#include <math.h>
+#include <TCut.h>
+#include <TROOT.h>
+#include <TProfile.h>
+#include <TLegendEntry.h>
-+
+#include <RooRealVar.h>
-+
+#include <RooArgSet.h>
-+
+#include <RooDataSet.h>
-+
+#include <RooDataHist.h>
-+
+#include <RooMCStudy.h>
-+
+#include <RooCategory.h>
-+
-+
+#include <RooPlot.h>
-+
+#include <RooSimultaneous.h>
-+
+#include <RooAddPdf.h>
-+
+#include <RooFitResult.h>
-+
+#include <RooVoigtian.h>
-+
+#include <RooChebychev.h>
-+
+#include <RooMsgService.h>
-+
-+
+#include <RooStats/ModelConfig.h>
-+
+#include "RooStats/ProfileLikelihoodCalculator.h"
-+
+#include "RooStats/LikelihoodInterval.h"
-+
+#include "RooStats/HypoTestResult.h"
-+
+#include "RooStats/SimpleLikelihoodRatioTestStat.h"
-+
+#include "RooStats/ProfileLikelihoodTestStat.h"
-+
+#include "RooStats/HybridCalculatorOriginal.h"
-+
+#include "RooStats/HypoTestInverterOriginal.h"
-+
-+
-+
+using namespace std;
+using namespace PlottingSetup;
+namespace ZbData {
-<
+ vector<float> data_over_mc;
->
+ TCut ZplusBsel("pt1>20&&pt2>20&&mll>60&&mll<120&&id1==id2");
->
+ TCut LeadingB("ZbCHS3010_bTagProbCSVBP[0]>0.898");
->
+ TCut EtaB("abs(ZbCHS3010_pfJetGoodEta[0])<1.3");
->
+ TCut PhiZcut("abs(ZbCHS3010_pfJetDphiZ[0])>2.7");
->
+ float LowerHistoBoundary=0.5; // where to cut on MPF (or Rabs)
->
+ float UpperHistoBoundary=1.5; // where to cut on MPF (or Rabs)
-–
+ TCut ZplusBsel("pt1>20&&pt2>20&&mll>60&&mll<120");
-–
+ TCut LeadingB("bTagProbCSVBP[0]>0.679");
-–
+ TCut EtaB("abs(pfJetGoodEta[0])<1.3");
-–
+ TCut PhiZcut("pfJetDphiZ[0]>2.7");
-–
-+
+TCut STDWEIGHT;
-+
-+
+class ZbPointResult {
-+
+public:
-+
+  ZbPointResult(string var, float ptlow, float pthigh, float alphalow, float alphahigh, Value Data, Value MC);
-+
+  float PtLow;
-+
+  float PtHigh;
-+
+  float AlphaLow;
-+
+  float AlphaHigh;
-+
+  string variable;
-+
+  Value Data;
-+
+  Value MC;
-+
+  Value Ratio;
-+
+};
-+
-+
+ZbPointResult::ZbPointResult(string var, float ptlow, float pthigh, float alphalow, float alphahigh, Value data, Value mc) :
-+
+  variable(var), PtLow(ptlow), PtHigh(pthigh), AlphaLow(alphalow), AlphaHigh(alphahigh), Data(data), MC(mc) {
-+
+    Ratio = data/mc;
-+
+}
-+
-+
+std::ostream &operator<<(std::ostream &ostr, ZbPointResult &b)
-+
+{
-+
+  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();
-+
+  return ostr;
-+
+}
-+
-+
-+
+class ZbResultContainer {
-+
+public:
-+
+  string name;
-+
+  vector<ZbPointResult> MPF_Rabs_location;
-+
-+
+  Value MPF_Result_FaceValue;
-+
+  Value Rabs_Result_FaceValue;
-+
+  Value MPF_Result_Extrapolated;
-+
+  Value Rabs_Result_Extrapolated;
-+
-+
+  vector<Value> MPF_Result_FaceValue_Bins;
-+
+  vector<Value> Rabs_Result_FaceValue_Bins;
-+
+  vector<Value> MPF_Result_Extrapolated_Bins;
-+
+  vector<Value> Rabs_Result_Extrapolated_Bins;
-+
-+
+  void Print();
-+
-+
+  ZbResultContainer(string name);
-+
+};
-+
-+
-+
+ZbResultContainer::ZbResultContainer(string _name) {
-+
+  name=_name;
-+
+  dout<< "Result Container for \"" << name << "\" has been initialized" << endl;
-+
+}
-+
-+
+void ZbResultContainer::Print()
-+
+{
-+
+  dout << "Results for " << this->name << endl;
-+
+  dout << " MPF and Rabs locations: " << endl;
-+
+  for(int i=0;i<this->MPF_Rabs_location.size();i++) dout << "   " << this->MPF_Rabs_location[i] << endl;
-+
+  dout << " Results: " << endl;
-+
+  dout << "   Face Value: " << endl;
-+
+  dout << "     MPF " << this->MPF_Result_FaceValue << endl;
-+
+  dout << "        Bins: [ ";
-+
+  if(MPF_Result_FaceValue_Bins.size()>0) {
-+
+    for(int i=0;i<MPF_Result_FaceValue_Bins.size()-1;i++) dout << MPF_Result_FaceValue_Bins[i] << " , ";
-+
+    if(MPF_Result_FaceValue_Bins.size()>0) dout << MPF_Result_FaceValue_Bins[MPF_Result_FaceValue_Bins.size()-1] << " ]" << endl;
-+
+  }
-+
-+
+  dout << "     Rabs " << this->Rabs_Result_FaceValue << endl;
-+
+  if(Rabs_Result_FaceValue_Bins.size()>0) {
-+
+    dout << "        Bins: [ ";
-+
+    for(int i=0;i<Rabs_Result_FaceValue_Bins.size()-1;i++) dout << Rabs_Result_FaceValue_Bins[i] << " , ";
-+
+    if(Rabs_Result_FaceValue_Bins.size()>0) dout << Rabs_Result_FaceValue_Bins[Rabs_Result_FaceValue_Bins.size()-1] << " ]" << endl;
-+
+  }
-+
-+
+  dout << "   Extrapolated: " << endl;
-+
+  dout << "     MPF " << this->MPF_Result_Extrapolated << endl;
-+
+  if(MPF_Result_Extrapolated_Bins.size()>0) {
-+
+    dout << "        Bins: [ ";
-+
+    for(int i=0;i<MPF_Result_Extrapolated_Bins.size()-1;i++) dout << MPF_Result_Extrapolated_Bins[i] << " , ";
-+
+    if(MPF_Result_Extrapolated_Bins.size()>0) dout << MPF_Result_Extrapolated_Bins[MPF_Result_Extrapolated_Bins.size()-1] << " ]" << endl;
-+
+  }
-+
-+
+  dout << "             Rabs " << this->Rabs_Result_Extrapolated << endl;
-+
+  if(Rabs_Result_Extrapolated_Bins.size()>0) {
-+
+    dout << "        Bins: [ ";
-+
+    for(int i=0;i<Rabs_Result_Extrapolated_Bins.size()-1;i++) dout << Rabs_Result_Extrapolated_Bins[i] << " , ";
-+
+    if(Rabs_Result_Extrapolated_Bins.size()>0) dout << Rabs_Result_Extrapolated_Bins[Rabs_Result_Extrapolated_Bins.size()-1] << " ]" << endl;
-+
+  }
-+
+}
-+
-+
+using namespace ZbData;
-<
+void print_yield(TCut cut) {
->
->
+//*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*  DELETE EVERYTHING BETWEEN THIS LINE /*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*//
->
+void SpecialCutFlow(string identifier);
->
->
+//*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*            AND THIS LINE             /*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*//
->
->
->
+class CutYields{
->
+public:
->
+  float data;
->
+  float Zb;
->
+  float Zc;
->
+  float Zl;
->
+  float tt;
->
+  float SiTo;
->
+  float Dibosons;
->
+  float WJets;
->
+  float Other;
->
->
+  CutYields();
->
+};
->
->
+CutYields::CutYields() {
->
+  this->data=0;
->
+  this->Zb=0;
->
+  this->Zc=0;
->
+  this->Zl=0;
->
+  this->tt=0;
->
+  this->SiTo=0;
->
+  this->Dibosons=0;
->
+  this->WJets=0;
->
+  this->Other=0;
->
+}
->
->
+std::ostream &operator<<(std::ostream &ostr, CutYields y)
->
+{
->
+  return ostr << y.data << ";" << y.Zb << ";" << y.Zc << ";" << y.Zl << ";" << y.tt << ";" << y.SiTo << ";" << y.Dibosons << ";" << y.WJets << ";" << endl;
->
+}
->
->
->
+CutYields print_yield(TCut cut) {
->
->
+  CutYields Yield;
+  TH1F *data  = allsamples.Draw("data",    "mll",1,0,10000000, "m_{ll} [GeV]", "events", cut,data,luminosity);
-<
+  cout << "  Yields for " << (const char*) cut << endl;
-<
+  cout << "     data: " << data->Integral() << endl;
->
+  dout << "  Yields for " << (const char*) cut << endl;
->
+//  dout << "     data: " << data->Integral() << endl;
->
->
+  Yield.data=data->Integral();
+  THStack mcm = allsamples.DrawStack("mc", "mll",1,0,10000000, "m_{ll} [GeV]", "events", cut,mc,luminosity);
+  int nhists = mcm.GetHists()->GetEntries();
-–
+  cout << "Going to loop over " << nhists << " histograms " << endl;
+  TFile *test = new TFile("test.root","RECREATE");
+  mcm.Write();
+  test->Close();
+  for (int i = 0; i < nhists; i++) {
+    TH1* hist = static_cast<TH1*>(mcm.GetHists()->At(i));
-<
+    cout << "     " << hist->GetName() << ": " << hist->Integral() << endl;
->
+    string name=hist->GetName();
->
+    if(Contains(name,"t_bar_t")) Yield.tt+=hist->Integral();
->
+    if(Contains(name,"W_Jets")) Yield.WJets+=hist->Integral();
->
+    if(Contains(name,"Single_top")) Yield.SiTo+=hist->Integral();
->
+    if(Contains(name,"Dibosons")) Yield.Dibosons+=hist->Integral();
->
+    if(Contains(name,"Z_l_")) Yield.Zl+=hist->Integral();
->
+    if(Contains(name,"Z_c_")) Yield.Zc+=hist->Integral();
->
+    if(Contains(name,"Z_b_")) Yield.Zb+=hist->Integral();
-<
+  cout << endl << endl;
->
->
+  dout << Yield << endl;
->
+  dout << endl << endl;
->
+  delete data;
-+
+  CleanLegends();
-+
+  DeleteStack(mcm);
-+
-+
+  return Yield;
+void draw_kin_variable(string variable, TCut cut, int nbins, float min, float max, string xlabel, string saveas, bool dologscale, string speciallabel="") {
-–
+  TCanvas *ckin = new TCanvas("ckin","kin variable canvas");
-<
+  ckin->SetLogy(dologscale);
->
+  if(dologscale) ckin->SetLogy(1);
+  TH1F *datahisto = allsamples.Draw("datahisto",variable,nbins,min,max, xlabel, "events",cut,data,luminosity);
+  datahisto->SetMarkerSize(DataMarkerSize);
+  THStack mcstack = allsamples.DrawStack("mcstack",variable,nbins,min,max, xlabel, "events",cut,mc,luminosity);
+  if(dologscale) datahisto->SetMaximum(3*datahisto->GetMaximum());
+  else datahisto->SetMaximum(1.5*datahisto->GetMaximum());
-<
+  TText *speciall;
-<
+  if(speciallabel!="") {
-<
+    speciall = new TText(0.2,0.8,speciallabel.c_str());
-<
+    speciall->Draw();
->
->
+  float SumMC=0.0;
->
+  TIter nextH(mcstack.GetHists());
->
+  TH1F* h;
->
+  while ( h = (TH1F*)nextH() ) {
->
+    SumMC += h->Integral();
->
+  }
->
+  float scale = datahisto->Integral()/SumMC;
->
->
+  // Re-scale stack to data integral
->
+  nextH = TIter(mcstack.GetHists());
->
->
+  while ( h = (TH1F*)nextH() ) {
->
+    h->Scale(scale);
-+
-+
+  TText *speciall;
+  datahisto->Draw("e1");
+  ckin->Update();
-<
+  mcstack.Draw("same");
->
+  mcstack.Draw("histo,same");
+  datahisto->Draw("same,e1");
+  TLegend *kinleg = allsamples.allbglegend();
-+
+  kinleg->Draw();
+  TPad *kinpad = new TPad("kinpad","kinpad",0,0,1,1);
+  kinpad->cd();
-<
+  kinpad->SetLogy(dologscale);
->
+  if(dologscale) kinpad->SetLogy(1);
+  datahisto->Draw("e1");
-<
+  mcstack.Draw("same");
->
+  mcstack.Draw("histo,same");
+  datahisto->Draw("same,e1");
+  datahisto->Draw("same,axis");
+  kinleg->Draw();
+  DrawPrelim();
+  saveas="kin/"+saveas;
-<
+  save_with_ratio(datahisto,mcstack,kinpad->cd(),saveas);
-<
-<
->
+  kinpad->cd();
-+
+  if(speciallabel!="") {
-+
+    speciall = write_title(speciallabel);
-+
+    speciall->SetX(0.18);
-+
+    speciall->SetTextAlign(11);
-+
+    speciall->SetTextSize(0.03);
-+
+    speciall->SetY(0.85);
-+
+    speciall->Draw();
-+
+  }
-<
+  ZbData::data_over_mc.push_back(datahisto->Integral()/CollapseStack(mcstack)->Integral());
->
+  Save_With_Ratio(datahisto,mcstack,kinpad->cd(),saveas);
+  datahisto->Delete();
-+
+  delete kinleg;
+  delete ckin;
-+
+  CleanLegends();
-+
+  DeleteStack(mcstack);
-+
-+
+void print_all_b_yields() {
-+
+  vector<CutYields> yields;
-+
+  dout << "Basic selection with a b jet" << endl;
-+
+  yields.push_back(print_yield(ZplusBsel&&"ZbCHS3010_pfJetGoodNumBtag>0"));
-+
+  dout << "Leading jet is a b-jet " << endl;
-+
+  yields.push_back(print_yield(ZplusBsel&&LeadingB));
-+
+  dout << "|#eta_{b}|<1.3 " << endl;
-+
+  yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB));
-+
+  dout << "|#delta#phi(b<Z)|>2.7" << endl;
-+
+  yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut));
-+
+  dout << "30<ptZ<1000 GeV" << endl;
-+
+  yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"));
-+
+  dout << "#alpha < 0.35" << endl;
-+
+  yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"(ZbCHS3010_alpha)<0.35"));
-+
+  dout << "#alpha < 0.3" << endl;
-+
+  yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"(ZbCHS3010_alpha)<0.3"));
-+
+/*  dout << "#alpha < 0.2" << endl;
-+
+  yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"(ZbCHS3010_alpha)<0.2"));
-+
+  dout << "#alpha < 0.15" << endl;
-+
+  yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"(ZbCHS3010_alpha)<0.15"));
-+
+  dout << "#alpha < 0.1" << endl;
-+
+  yields.push_back(print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000"&&"(ZbCHS3010_alpha)<0.1"));*/
-+
-+
+  for(int iy=0;iy<yields.size();iy++) {
-+
+    dout << yields[iy] << endl;
-+
+  }
-+
+}
-+
-+
+void DrawEvilCutFlow() {
-+
+  stringstream MegaCut;
-<
+Value getfrom2Dmap(TH2F *map, int ixbin, int iybin) {
-<
+  Value sum(0,0);
-<
+  for(int iy=1;iy<=iybin;iy++) {
-<
+    sum=sum+Value(map->GetBinContent(ixbin,iy),map->GetBinError(ixbin,iy));
-<
+  }
-<
+  return sum;
-<
+}
->
+  MegaCut << "   1*(" << (const char*) (ZplusBsel&&TCut("ZbCHS3010_pfJetGoodNumBtag>0")) << ")";
->
+  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB) << ")";
->
+  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB) << ")";
->
+  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut) << ")";
->
+  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000")) << ")";
->
+  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.3")) << ")";
->
+  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.2")) << ")";
->
+  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.15")) << ")";
->
+  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.1")) << ")";
->
+  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000&&ZbCHS3010_alpha<0.05")) << ")";
->
+  MegaCut << " ";
-+
+  TCanvas *can = new TCanvas("can","can");
-+
+  can->SetLogy(1);
-+
+  TCut basecut("mll>6||mll<10");
-<
+void data_mc_agreement_2d() {
-<
-<
+  vector<int> emptyvector;
-<
+  float binX[5]={10,50,100,200,1000};
-<
+  float binY[5]={0.0,0.1,0.15,0.2,0.3};
-<
+  vector<float>binningx;
-<
+  vector<float>binningy;
-<
+  binningx.assign(binX,binX+5);
-<
+  binningy.assign(binY,binY+5);
-<
-<
+  TCut cut(ZplusBsel&&LeadingB&&EtaB&&PhiZcut);
-<
-<
+  TCanvas *agrcan = new TCanvas("agrcan","agrcan");
-<
-<
+  cout << "Drawing 2d for data " << endl;
-<
+  TH2F *alphazptplane_data = allsamples.Draw("alphazptplane_data","alpha:pt",binningx,binningy,"Z p_{T}","#alpha",cut,data,luminosity,emptyvector,false);
-<
+  cout << "Drawing 2d for MC " << endl;
-<
+  TH2F *alphazptplane_mc   = allsamples.Draw("alphazptplane_mc",  "alpha:pt",binningx,binningy,"Z p_{T}","#alpha",cut,mc,  luminosity,emptyvector,false);
-<
-<
-<
+  TFile *f = new TFile("test.root","RECREATE");
-<
+  alphazptplane_data->Write();
-<
+  f->Close();
-<
-<
+  for(int ix=1;ix<binningx.size();ix++) {
-<
+    cout << "pt in " << binningx[ix-1] << " , " << binningx[ix] << endl;
-<
+    TGraphErrors gr;
-<
+    for(int iy=1;iy<binningy.size();iy++) {
-<
+      Value mc   = Value((getfrom2Dmap(alphazptplane_mc,ix,iy)).getValue(),0.0);// not taking MC statistics into account
-<
+      Value data = getfrom2Dmap(alphazptplane_data,ix,iy);
-<
+      cout << " For an alpha cut at " << binningy[iy] << " we get : mc (" << mc << ") , data (" << data << ") --> data/mc = " << data/mc << endl;
-<
+      gr.SetPoint(iy-1,binningy[iy],(data/mc).getValue());
-<
+      gr.SetPointError(iy-1,0,(data/mc).getError());
-<
+    }
-<
-<
+    gr.Draw("AP*");
-<
+    gr.GetYaxis()->SetTitle("data / mc");
-<
+    gr.GetYaxis()->CenterTitle();
-<
+    gr.GetXaxis()->SetTitle("#alpha");
-<
+    gr.GetXaxis()->CenterTitle();
-<
+    gr.GetXaxis()->SetRangeUser(0,0.4);
-<
+    gr.GetYaxis()->SetRangeUser(0.5,1.5);
-<
+    gr.Draw("AP*");
-<
+    DrawPrelim();
-<
+    stringstream ptbin;
-<
+    ptbin << binningx[ix-1] << " < p_{T}^{Z} < " << binningx[ix] << " GeV";
-<
+    TText *ptt = write_title(ptbin.str());
-<
+    ptt->SetX(0.7);
-<
+    ptt->SetY(0.7);
-<
+    ptt->SetTextSize(0.03);
-<
+    ptt->Draw();
-<
+    TLine *l = new TLine(0.0,1.0,0.4,1.0);
-<
+    l->SetLineStyle(2);
-<
+    l->SetLineColor(kBlue);
-<
+    l->Draw();
-<
-<
+    stringstream saveas;
-<
+    saveas << "Zplusb_data_over_mc___" << binningx[ix-1] << "to" << binningx[ix];
-<
+    DrawPrelim();
-<
+    CompleteSave(agrcan,saveas.str());
-<
+//    agrcan->SaveAs("tester.png");
->
->
->
+  TLegend *leg = allsamples.allbglegend();
->
->
->
+  TH1F *data  = allsamples.Draw("data",    MegaCut.str(),11,-0.5,10.5, "", "events", basecut,data,luminosity);
->
+  THStack mcm = allsamples.DrawStack("mc", MegaCut.str(),11,-0.5,10.5, "", "events", basecut,mc,luminosity);
->
->
+  float runningsum=0;
->
+  for(int i=data->GetNbinsX();i>0;i--) {
->
+    runningsum+=data->GetBinContent(i);
->
+    data->SetBinContent(i,runningsum);
->
+  }
->
->
+  float minimum=data->GetMaximum();
->
->
+  TH1F* h;
->
+  TIter nextOF(mcm.GetHists());
->
+  while ( h = (TH1F*)nextOF() ) {
->
+    float runningsum=0;
->
+    minimum=data->GetMaximum();//done deliberately at each step so get the minimum of the last (!) contribution
->
+    for(int i=h->GetNbinsX();i>0;i--) {
->
+      runningsum+=h->GetBinContent(i);
->
+      h->SetBinContent(i,runningsum);
->
+      if(runningsum<minimum&&runningsum>0) minimum=runningsum;
->
+    }
-–
+  delete agrcan;
-–
+}
-<
+void print_all_b_yields() {
-<
+  cout << "Basic selection with a b jet" << endl;
-<
+  print_yield(ZplusBsel&&"pfJetGoodNumBtag>0");
-<
+  cout << "Leading jet is a b-jet " << endl;
-<
+  print_yield(ZplusBsel&&LeadingB);
-<
+  cout << "|#eta_{b}|<1.3 " << endl;
-<
+  print_yield(ZplusBsel&&LeadingB&&EtaB);
-<
+  cout << "|#delta#phi(b<Z)|>2.7" << endl;
-<
+  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut);
-<
+  cout << "10<ptZ<1000 GeV" << endl;
-<
+  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000");
-<
+  cout << "#alpha < 0.3" << endl;
-<
+  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"alpha<0.3");
-<
+  cout << "#alpha < 0.2" << endl;
-<
+  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"alpha<0.2");
-<
+  cout << "#alpha < 0.15" << endl;
-<
+  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"alpha<0.15");
-<
+  cout << "#alpha < 0.1" << endl;
-<
+  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"alpha<0.1");
-<
+  cout << "#alpha < 0.05" << endl;
-<
+  print_yield(ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"alpha<0.05");
-<
+}
->
+  data->SetMinimum(0.05*minimum);
->
+  can->cd(1)->SetBottomMargin(0.25);
->
+  data->GetXaxis()->SetBinLabel(1,"Pre-selection");
->
+  data->GetXaxis()->SetBinLabel(2,"Contains b jet");
->
+  data->GetXaxis()->SetBinLabel(3,"Leading jet is b-jet");
->
+  data->GetXaxis()->SetBinLabel(4,"|#eta_{b}|<1.3 ");
->
+  data->GetXaxis()->SetBinLabel(5,"|#delta#phi(b,Z)|>2.7");
->
+  data->GetXaxis()->SetBinLabel(6,"30<p_{T}^{Z}<1000 GeV");
->
+  data->GetXaxis()->SetBinLabel(7,"#alpha < 0.3");
->
+  data->GetXaxis()->SetBinLabel(8,"#alpha < 0.2");
->
+  data->GetXaxis()->SetBinLabel(9,"#alpha < 0.15");
->
+  data->GetXaxis()->SetBinLabel(10,"#alpha < 0.1");
->
+  data->GetXaxis()->SetBinLabel(11,"#alpha < 0.05");
->
+  data->GetXaxis()->LabelsOption("v");
->
->
+  data->Draw("e1");
->
+  mcm.Draw("histo,same");
->
+  data->Draw("same,e1,axis");
->
+  data->Draw("same,e1");
->
+//  data->Draw("same,TEXT");
->
->
+  leg->Draw();
->
->
+  DrawPrelim();
->
->
+  CompleteSave(can,"CutFlow");
->
+  CleanLegends();
->
+  DeleteStack(mcm);
->
->
+}
+void draw_Zb_kin_vars() {
-<
+  draw_kin_variable("pfBJetDphiZ",ZplusBsel&&LeadingB,100,0,3.2,"#delta#phi (Z,b lead)","DeltaPhiZBlead",1);
-<
+  draw_kin_variable("pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",50,0,200,"Z p_{T}","Zpt",1);
-<
+  draw_kin_variable("pfJetGoodPt[1]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",50,0,200,"Sub-Leading Jet Pt","Jet2Pt",1);
-<
+  draw_kin_variable("pfJetGoodPtlBtag[0]",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",50,0,200,"Leading Jet Pt (B)","LeadingJetPt",1);
-<
+  draw_kin_variable("alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000",20,0,2,"#alpha","alpha",1);
->
->
+   draw_kin_variable("pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000"),25,0,200,"Z p_{T}","Official/Zpt",1);
->
+   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);
->
+   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);
->
+   draw_kin_variable("pt",ZplusBsel&&EtaB&&PhiZcut&&TCut("pt>30&&pt<1000"),25,0,200,"Z p_{T}","Official/Zpt__nonBtagged",1);
->
+   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);
->
+   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);
->
->
+   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);
->
+   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);
->
+   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);
->
+   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);
->
+   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);
->
+   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);
->
->
+   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);
->
+   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);
->
+   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);
->
+   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);
->
+   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);
->
+   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);
->
->
->
+   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);
->
+   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);
->
+   draw_kin_variable("ZbCHS3010_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut("pt>50&&pt<75"),  40,0,2,"#alpha","Official/alpha_50_to_75",1);
->
+   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);
->
+   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);
->
+   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);
->
->
+   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);
+void draw_mpf_vars() {
-<
+  ZbData::data_over_mc.clear();
-<
+  draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"alpha<0.3",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p3",0,"#alpha<0.3, 10 GeV < p_{T}^{Z} < 1000 GeV");
-<
+  draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"alpha<0.2",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p2",0,"#alpha<0.2, 10 GeV < p_{T}^{Z} < 1000 GeV");
-<
+  draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"alpha<0.15",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p15",0,"#alpha<0.15, 10 GeV < p_{T}^{Z} < 1000 GeV");
-<
+  draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"alpha<0.1",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p1",0,"#alpha<0.1, 10 GeV < p_{T}^{Z} < 1000 GeV");
-<
+  draw_kin_variable("mpf",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"alpha<0.05",20,0,2,"Z+b MPF","mpf_alpha_smaller_0p05",0,"#alpha<0.05, 10 GeV < p_{T}^{Z} < 1000 GeV");
-<
-<
+  draw_kin_variable("pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"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");
-<
+  draw_kin_variable("pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"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");
-<
+  draw_kin_variable("pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"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");
-<
+  draw_kin_variable("pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"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");
-<
+  draw_kin_variable("pfJetGoodPt[0]/pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>10&&pt<1000"&&"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");
->
+  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");
->
+  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");
->
+  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");
->
+  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");
->
+  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");
->
->
+  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");
->
+  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");
->
+  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");
->
+  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");
->
+  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");
-<
+Value get_Zb_data_over_mc(string variable, TCut cut) {
-<
+//  write_warning(__FUNCTION__,"Debugging this function, therefore always returning 3.1415+/-1.010101"); return Value(3.1415,1.010101);
-<
+  TH1F *hdata  = allsamples.Draw("hdata",variable,1000,0,30, "#alpha [GeV]", "events", cut,data,luminosity);
-<
+  TH1F *hmc    = allsamples.Draw("hmc" , variable,1000,0,30, "#alpha [GeV]", "events", cut,mc,  luminosity);
->
+void make_RooFit(string ContainerName, TH1F *h, float &peak, float &error_peak, string saveas) {
->
+  TCanvas *fcan = new TCanvas("fcan","fcan");
->
+  RooRealVar mpf("mpf","mpf", -5, 5, "");
->
+  RooDataHist AllData("AllData", "CMS Real Dataset", mpf , h);
->
->
+  RooRealVar fraction("fraction", "fraction", 0.8, 0, 1);
->
+  RooRealVar mean("mean", "mean", 1, 0, 3);
->
+  RooRealVar sigma("sigma", "sigma", 1, 0, 5);
->
+  RooRealVar width("width", "width", 2.4, 0, 5);
->
+  RooRealVar chevvar1("chevvar1", "chevvar1", 0.2,-1.0, 1.0);
->
+  RooRealVar chevvar2("chevvar2", "chevvar2", 0.2,-1.0, 1.0);
->
->
+  RooVoigtian signal("signal", "signal", mpf, mean, width, sigma);
->
+  RooChebychev background("background_", "background", mpf, RooArgList(chevvar1, chevvar2));
->
->
+  RooAddPdf pdf("pdf",  "pdf", RooArgList(signal, background), fraction);
->
->
+  RooFitResult *result = pdf.fitTo(AllData, RooFit::Save());
->
+  RooPlot* frame = mpf.frame(RooFit::Bins(30), RooFit::Title("MPF")) ;
->
+  AllData.plotOn(frame);
->
+  pdf.plotOn(frame, RooFit::LineColor(kBlue)) ;
->
+  pdf.plotOn(frame, RooFit::Components(RooArgSet(background)), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed)) ;
->
+  pdf.plotOn(frame, RooFit::Components(RooArgSet(signal)), RooFit::LineStyle(kDashed), RooFit::LineColor(kGreen)) ;
->
->
+  RooRealVar* thePeak = (RooRealVar*) result->floatParsFinal().find("mean");
->
+  peak = thePeak->getVal();
->
+  error_peak = thePeak->getError();
->
->
+  ///******************* just making it beautiful ....
->
+  RooPlot* frame2 = mpf.frame(RooFit::Bins(30), RooFit::Title("MPF")) ;
->
+  TH1F *hcopy = (TH1F*)h->Clone("hcopy");
->
+  hcopy->Rebin(40);
->
+  RooDataHist AllRData("AllRData", "CMS Dataset", mpf , hcopy);
->
->
+  AllRData.plotOn(frame2);
->
+  pdf.plotOn(frame2, RooFit::LineColor(kBlue)) ;
->
+  pdf.plotOn(frame2, RooFit::Components(RooArgSet(background)), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed)) ;
->
+  pdf.plotOn(frame2, RooFit::Components(RooArgSet(signal)), RooFit::LineStyle(kDashed), RooFit::LineColor(kGreen)) ;
->
->
+  string saveasFinal = saveas + "_Fit";
->
+  fcan->cd() ; gPad->SetLeftMargin(0.15) ; frame->GetYaxis()->SetTitleOffset(1.4) ; frame2->Draw();
->
->
+  if(Contains(saveas,"_data")) DrawPrelim();
->
+  else DrawMCPrelim();
->
->
+  CompleteSave(fcan, ContainerName+"/Fit/"+saveasFinal);
->
->
+  delete hcopy;
->
+  delete fcan;
->
+  delete frame;
->
+  delete thePeak;
->
+}
->
->
+float bother=0.1;
->
+float bother2=1;
->
+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) {
->
+  TCanvas *cn = new TCanvas("cn","cn");
->
+  string varname="MPF";
->
->
+  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;
->
+  //note that "15.0" here is the subleading jet cut.
->
->
+  bool DoFit=false;
->
+  bool UseFit=false;
->
+  bool DoRooFit=false;
->
->
->
+  if(Contains(variable,"pfJetGood")) varname="R_{abs}";
->
+//  TH1F *hdata  = allsamples.Draw("hdata",variable,1200,0,30, varname, "events",  cut,data,luminosity);
->
+//  TH1F *hmc    = allsamples.Draw("hmc" , variable,1200,0,30, varname, "events", cut, mc,  luminosity);
->
+  TH1F *hdata, *hmc;
->
+  if(!DoRooFit) {
->
+    hdata = allsamples.Draw("hdata",variable,1200,ZbData::LowerHistoBoundary,ZbData::UpperHistoBoundary, varname, "events",  cut,data,luminosity); // making it more precise
->
+    hmc    = allsamples.Draw("hmc" , variable,1200,ZbData::LowerHistoBoundary,ZbData::UpperHistoBoundary, varname, "events", cut, mc,  luminosity); // making it more precise
->
+  } else {
->
+    hdata  = allsamples.Draw("hdata",variable,7200,0.0,6.0, varname, "events",  cut,data,luminosity); // making it more precise
->
+    hmc    = allsamples.Draw("hmc" , variable,7200,0.0,6.0, varname, "events", cut, mc,  luminosity); // making it more precise
->
+  }
->
->
+  hmc->SetLineColor(kBlue);
+  float a=hdata->GetMean();
+  float b=hmc->GetMean();
+  float db=hmc->GetMeanError();
+  float factor = a / b;
+  float error = TMath::Sqrt(  (1/(b*b)) * (da*da) + ((a*a)/(b*b))*db*db);
-+
-+
+  TF1 *gaus_data = new TF1("gaus_data","gaus",0.7,1.3);
-+
+  TF1 *gaus_mc = new TF1("gaus_mc","gaus",0.7,1.3);
-+
+  if(DoFit) {
-+
+    gaus_data->SetParameter(0,hdata->GetMaximum());
-+
+    gaus_data->SetParameter(2,hdata->GetMean());
-+
+    gaus_data->SetParameter(1,hdata->GetRMS());
-+
+    hdata->Fit(gaus_data);
-+
-+
+    gaus_mc->SetParameter(0,hmc->GetMaximum());
-+
+    gaus_mc->SetParameter(1,hmc->GetMean());
-+
+    gaus_mc->SetParameter(2,hmc->GetRMS());
-+
+    hmc->Fit(gaus_mc);
-+
+  }
-+
+  if(DoRooFit) {
-+
-+
+      if(hdata->Integral()>10) {
-+
+        make_RooFit(ContainerName, hdata, a, da, saveas+"_data");
-+
+        make_RooFit(ContainerName, hmc, b, db, saveas+"_mc");
-+
+        factor = a / b;
-+
+        error = TMath::Sqrt(  (1/(b*b)) * (da*da) + ((a*a)/(b*b))*db*db);
-+
+      } else {
-+
+        factor=NAN;//set to nan
-+
+        error=NAN;//set to nan
-+
+      }
-+
+  }
-+
-+
+  cn->cd();
-+
+  hdata->Rebin(int(0.05*1200/(ZbData::UpperHistoBoundary-ZbData::LowerHistoBoundary)));
-+
+  hmc->Rebin(int(0.05*1200/(ZbData::UpperHistoBoundary-ZbData::LowerHistoBoundary)));
-+
+  hdata->GetYaxis()->SetTitle("events / 0.05");
-+
+  hdata->SetMarkerColor(kRed);
-+
+  if(DoFit) {
-+
+    gaus_mc->SetLineColor(kBlue);
-+
+    gaus_data->SetLineColor(kRed);
-+
+  }
-+
-+
+  hdata->GetXaxis()->SetRangeUser(0,2);
-+
+  hdata->GetYaxis()->SetTitleOffset(1.3);
-+
+  hdata->Draw("e1");
-+
-+
+  hmc->Scale(hdata->Integral()/hmc->Integral());
-+
+  if(hmc->GetMaximum()>hdata->GetMaximum()) hdata->SetMaximum((hmc->GetMaximum()+sqrt((double)hmc->GetMaximum()))*1.2);
-+
-+
+  hmc->Draw("histo,same");
-+
+  hdata->Draw("e1,same");
-+
+  if(DoFit) {
-+
+    gaus_data->Draw("same");
-+
+    gaus_mc->Draw("same");
-+
+  }
-+
-+
+  stringstream summary;
-+
+  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 << "}{";
-+
+  summary << "#splitline{mc:   " << std::setprecision(4) << b << " #pm " << std::setprecision(4) << db << "}{#splitline{ratio: " << factor;
-+
+  summary << " #pm " << error;
-+
+  if(DoFit) {
-+
+    summary << "}{#splitline{}{#splitline{data fit: " << std::setprecision(4) << gaus_data->GetParameter(1) << "+/-";
-+
+    summary << std::setprecision(4) << gaus_data->GetParError(1) << "}{#splitline{";
-+
+    summary << "mc fit:" << std::setprecision(4) << gaus_mc->GetParameter(1) << "#pm" << std::setprecision(4) << gaus_mc->GetParError(1) << "}{ratio: ";
-+
+    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)) ) << "}}}}";
-+
+    summary << "}}}}";
-+
+  } else {
-+
+    summary << "}{}}}}}";
-+
+  }
-+
-+
+  if(UseFit) {
-+
+    write_error(__FUNCTION__,"Requested using fit result for final result - not implemented yet!");
-+
+    assert(0);
-+
+  }
-+
-+
-+
+  TText *infobox = write_title(summary.str());
-+
+  infobox->SetX(0.75);
-+
+  infobox->SetTextSize(0.03);
-+
+  infobox->SetY(0.75);
-+
+  infobox->Draw();
-+
-+
-+
+  DrawPrelim();
-+
+  CompleteSave(cn,ContainerName+"/ResponsePlots/"+saveas);
-+
-+
+  res.MPF_Rabs_location.push_back(ZbPointResult(variable,pt1,pt2,a1,a2,Value(a,da),Value(b,db)));
-+
+  dout << "Data;" << a << ";"<<da<<";MC;"<<b<<";"<<db<<endl;
-+
+  delete cn;
+  delete hdata;
+  delete hmc;
-+
+  return Value(factor,error);
-+
+int nFakes=1;
-+
-+
+float PurityInclusive=0.05;
-+
+float PurityLoose=0.20;
-+
+float PurityMedium=0.54;
-+
+float PurityTight=0.82;
-+
-+
-+
+void GetPtPointPosition(float AvPos[], const int nptcuts,string AlphaVariable, TCut basecut, float FaceValueToBeTakenAt, float ptcuts[]) {
-+
+  for(int i=0;i<nptcuts-1;i++) {
-+
+    stringstream specialcut;
-+
+    specialcut << "((pt>30&& pt< 1000) && (" << AlphaVariable << "<" << FaceValueToBeTakenAt << "))";
-+
+    TH1F *hdata  = allsamples.Draw("hada", "pt",100,ptcuts[i],ptcuts[i+1], "MPF", "events", TCut(basecut && specialcut.str().c_str()), data,  luminosity);
-+
+    AvPos[i]=hdata->GetMean();
-+
+    dout << AvPos[i] << " : ";
-+
+    delete hdata;
-+
+  }
-+
+  dout << endl;
-+
+}
-+
-+
-+
+void DeterminePurity(TCut basecut, int nptcuts, float ptcuts[], string AlphaVariable, string ContainerName, float FaceValueToBeTakenAt) {
-+
-+
+  bool ComputePurity=false;
-+
-+
+  if(ContainerName=="inclusive") ComputePurity=true;
-+
+  if(ContainerName=="loose") ComputePurity=true;
-+
+  if(ContainerName=="medium") ComputePurity=true;
-+
+  if(ContainerName=="reference") ComputePurity=true;
-+
-+
+  if(!ComputePurity) return;
-+
-+
+  stringstream specialcut;
-+
+  specialcut << "((pt>30&& pt< 1000) && (" << AlphaVariable << "<" << FaceValueToBeTakenAt << "))";
-+
-+
+  THStack mcm    = allsamples.DrawStack("mcm" , "mpf",1200,0,30, "MPF", "events", TCut(basecut && specialcut.str().c_str()), mc,  luminosity);
-+
+  int nhists = mcm.GetHists()->GetEntries();
-+
-+
+  float total=0;
-+
+  float Zb=0;
-+
-+
+  for (int i = 0; i < nhists; i++) {
-+
+    TH1* hist = static_cast<TH1*>(mcm.GetHists()->At(i));
-+
+    string name=hist->GetName();
-+
+    if(Contains(name,"t_bar_t")) total+=hist->Integral();
-+
+    if(Contains(name,"W_Jets")) total+=hist->Integral();
-+
+    if(Contains(name,"Single_top")) total+=hist->Integral();
-+
+    if(Contains(name,"Dibosons")) total+=hist->Integral();
-+
+    if(Contains(name,"Z_l_")) total+=hist->Integral();
-+
+    if(Contains(name,"Z_c_")) total+=hist->Integral();
-+
+    if(Contains(name,"Z_b_")) {
-+
+      total+=hist->Integral();
-+
+      Zb+=hist->Integral();
-+
+    }
-+
+  }
-+
+  float purity=Zb/total;
-+
+  dout << "Determined a purity of " << 100*purity << " % for the container \"" << ContainerName << "\"" << endl;
-+
-+
+  if(ContainerName=="inclusive") PurityInclusive=purity;
-+
+  if(ContainerName=="loose") PurityLoose=purity;
-+
+  if(ContainerName=="medium") PurityMedium=purity;
-+
+  if(ContainerName=="reference") PurityTight=purity;
-+
-+
+  DeleteStack(mcm);
-+
+}
-<
+void new_data_mc_agreement_2d() {
->
+string NiceAlphaPrint(float alpha) {
->
+  //get in : 0.3, return 0p3 ...
->
+  stringstream res;
->
+  res << int(alpha) << "p" << int(100*alpha);
->
+  return res.str();
->
+}
->
->
+ZbResultContainer new_data_mc_agreement_2d(bool gofast, string AlphaVariable, string ContainerName, TCut BTagCut, TCut BTagWgt, float FaceValueToBeTakenAt) {
->
->
+  if(gofast) write_info(__FUNCTION__,"Fast mode activated for "+ContainerName);
->
+  ZbResultContainer results(ContainerName);
->
->
+  cutWeight=STDWEIGHT*BTagWgt;
->
+  LeadingB=BTagCut;
->
->
+  gStyle->SetOptFit(0);
+  TCut basecut(ZplusBsel&&LeadingB&&EtaB&&PhiZcut);
-+
+  const int nalphacuts=6;
-+
+  float alphacuts[nalphacuts] = {0.1,0.15,0.2,0.25,0.3,0.35};
-+
-–
+  const int nalphacuts=5;
-–
+  float alphacuts[nalphacuts] = {0.1,0.15,0.2,0.25,0.3};
+  const int nptcuts=5;
-<
+  float ptcuts[nptcuts]={10,50,100,200,1000};
->
+  float ptcuts[nptcuts]={30,50,75,125,1000};
->
+  float ptPointPosition[nptcuts];
->
+  GetPtPointPosition(ptPointPosition,nptcuts,AlphaVariable,basecut,FaceValueToBeTakenAt,ptcuts);
->
->
->
+  DeterminePurity(basecut, nptcuts,ptcuts, AlphaVariable, ContainerName, FaceValueToBeTakenAt);
+  float MPF_Results[nalphacuts][nptcuts];
+  float MPF_Errors[nalphacuts][nptcuts];
+  float RABS_Errors[nalphacuts][nptcuts];
-<
+  for(int ia=0;ia<nalphacuts;ia++) {
->
+  for(int ia=0;ia<nalphacuts && !gofast;ia++) {
+    for(int ipt=0;ipt<nptcuts-1;ipt++) {
+      stringstream specialcut;
-<
+      specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (alpha<" << alphacuts[ia] << "))";
-<
+      Value MPF_data_over_mc = get_Zb_data_over_mc("mpf",TCut(basecut && specialcut.str().c_str()));
-<
+      Value RABS_data_over_mc = get_Zb_data_over_mc("pfJetGoodPt[0]/pt",TCut(basecut && specialcut.str().c_str()));
->
+      float alow,ahigh;
->
+      if(ia>0) {
->
+        specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (" << AlphaVariable << "<" << alphacuts[ia] << " && " << AlphaVariable << ">" << alphacuts[ia-1] << "))";
->
+        alow=alphacuts[ia-1];
->
+        ahigh=alphacuts[ia];
->
+      } else {
->
+        specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (" << AlphaVariable << "<" << alphacuts[ia] << "))";
->
+        alow=0;
->
+        ahigh=alphacuts[ia];
->
+      }
->
+      dout << specialcut.str() << endl;
->
+      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);
->
+      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);
+      MPF_Results[ia][ipt]=MPF_data_over_mc.getValue();
+      MPF_Errors[ia][ipt]=MPF_data_over_mc.getError();
+      RABS_Results[ia][ipt]=RABS_data_over_mc.getValue();
+      RABS_Errors[ia][ipt]=RABS_data_over_mc.getError();
-<
+      cout << "alpha cut at " << alphacuts[ia] << " 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;
->
+      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;
+  float RABS_ExtraPolatedResults[nptcuts];
+  TCanvas *can = new TCanvas("can");
-+
+  can->cd(1)->SetLeftMargin(0.15);
+  TGraphErrors *MPF_FinalGraph = new TGraphErrors();
+  TGraphErrors *RABS_FinalGraph = new TGraphErrors();
-+
+  TGraphErrors *MPF_FaceValueAtPoint3 = new TGraphErrors();
-+
+  TGraphErrors *RABS_FaceValueAtPoint3 = new TGraphErrors();
+  for(int ipt=0;ipt<nptcuts-1;ipt++) {
+    stringstream filename;
+    filename << "Extrapolation/Zplusb_data_over_mc___" << ptcuts[ipt] << "_to_" << ptcuts[ipt+1];
-<
+    cout << "Going to create histo for " << filename.str() << endl;
->
+    dout << "Going to create histo for " << filename.str() << endl;
+    TGraphErrors *mpf_gr =new TGraphErrors();
+    TGraphErrors *rabs_gr =new TGraphErrors();
-<
+    for(int ia=0;ia<nalphacuts;ia++) {
-<
+      mpf_gr->SetPoint(ia,alphacuts[ia],MPF_Results[ia][ipt]);
-<
+      mpf_gr->SetPointError(ia,0,MPF_Errors[ia][ipt]);
-<
+      rabs_gr->SetPoint(ia,alphacuts[ia],RABS_Results[ia][ipt]);
-<
+      rabs_gr->SetPointError(ia,0,RABS_Errors[ia][ipt]);
-<
+    }
-<
-<
+    mpf_gr->SetMarkerStyle(21);
-<
+    mpf_gr->SetMarkerColor(kBlue);
-<
+    mpf_gr->Draw("AP");
-<
+    mpf_gr->Fit("pol1","QE"); // quiet, use minos
-<
+    mpf_gr->GetXaxis()->SetTitle("#alpha");
-<
+    mpf_gr->GetXaxis()->CenterTitle();
-<
+    mpf_gr->GetYaxis()->SetTitle("data/mc");
-<
+    mpf_gr->GetYaxis()->CenterTitle();
-<
+    TF1 *mpf_pol=(TF1*)mpf_gr->GetFunction("pol1");
-<
+    float mpf_a=mpf_pol->GetParameter(0);
-<
+    float mpf_b=mpf_pol->GetParameter(1);
-<
+    MPF_FinalGraph->SetPoint(ipt,0.5*(ptcuts[ipt]+ptcuts[ipt+1]),mpf_a);
-<
+    MPF_FinalGraph->SetPointError(ipt,0,mpf_pol->GetParError(0));
-<
-<
+    MPF_ExtraPolatedResults[ipt]=mpf_a;
-<
-<
+    stringstream mpf_info;
-<
+    mpf_info <<  "#splitline{#splitline{#splitline{" << ptcuts[ipt] << " GeV < p_{T}^{Z} < " << ptcuts[ipt+1] << " GeV }{ (MPF) }}{Extrapolated value: " << std::setprecision(3) << mpf_a << "}}{#chi^{2}/NDF : " << mpf_pol->GetChisquare() << " / " << mpf_pol->GetNDF() << "}";
-<
-<
+    TText *mpf_mark = write_title(mpf_info.str());
-<
+    mpf_mark->SetX(0.75);
-<
+    mpf_mark->SetTextSize(0.03);
-<
+    mpf_mark->SetY(0.75);
-<
+    mpf_mark->Draw();
->
+    mpf_gr->SetName("MPF_Graph");
->
+    rabs_gr->SetName("Rabs_Graph");
-<
+    string filenamebkp=filename.str();
-<
+    filename << "__MPF";
-<
+    DrawPrelim();
-<
+    CompleteSave(can,filename.str());
-<
+    cout << "MPF : " << mpf_a << " + " << mpf_b << " * alpha " << endl;
->
+    int rabs_counter=0;
->
+    int mpf_counter=0;
-<
+    filename.str("");
->
+    for(int ia=0;ia<nalphacuts && !gofast;ia++) {
->
+      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;
->
+      if(MPF_Results[ia][ipt]==MPF_Results[ia][ipt] && MPF_Errors[ia][ipt]==MPF_Errors[ia][ipt]) { // remove nan's
->
+       mpf_gr->SetPoint(mpf_counter,alphacuts[ia],MPF_Results[ia][ipt]);
->
+       mpf_gr->SetPointError(mpf_counter,0,MPF_Errors[ia][ipt]);
->
+       mpf_counter++;
->
+      } else {
->
+        dout << "Detected BS for MPF method (bin ignored) : " << endl;
->
+        dout << "       alpha: " << alphacuts[ia] << " , pt range: " << ptcuts[ipt] << " < pt < " << ptcuts[ipt+1] << endl;
->
+        dout << "       value (MPF) : " << MPF_Results[ia][ipt] << " +/- " << MPF_Errors[ia][ipt] << endl;
->
+      }
->
+      if(RABS_Results[ia][ipt]==RABS_Results[ia][ipt] && RABS_Errors[ia][ipt]==RABS_Errors[ia][ipt]) {
->
+       rabs_gr->SetPoint(rabs_counter,alphacuts[ia],RABS_Results[ia][ipt]);
->
+       rabs_gr->SetPointError(rabs_counter,0,RABS_Errors[ia][ipt]);
->
+       rabs_counter++;
->
+      } else {
->
+        dout << "Detected BS for RABS method (bin ignored) : " << endl;
->
+        dout << "       alpha: " << alphacuts[ia] << " , pt range: " << ptcuts[ipt] << " < pt < " << ptcuts[ipt+1] << endl;
->
+        dout << "       value (RABS) : " << RABS_Results[ia][ipt] << " +/- " << RABS_Errors[ia][ipt] << endl;
->
+      }
->
+    }
->
->
+    stringstream specialcut;
->
+    specialcut << "((pt>" << ptcuts[ipt] << " && pt< " << ptcuts[ipt+1] << ") && (" << AlphaVariable << "<" << FaceValueToBeTakenAt << "))";
->
+    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);
->
+    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);
->
+    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;
->
+    MPF_FaceValueAtPoint3->SetPoint(ipt,(ptcuts[ipt]+ptcuts[ipt+1])/2,MPF_data_over_mc.getValue());
->
+    MPF_FaceValueAtPoint3->SetPointError(ipt,0,MPF_data_over_mc.getError());
->
+    RABS_FaceValueAtPoint3->SetPoint(ipt,(ptcuts[ipt]+ptcuts[ipt+1])/2,RABS_data_over_mc.getValue());
->
+    RABS_FaceValueAtPoint3->SetPointError(ipt,0,RABS_data_over_mc.getError());
-–
+    rabs_gr->SetMarkerStyle(21);
-–
+    rabs_gr->SetMarkerColor(kBlue);
-–
+    rabs_gr->Fit("pol1","QE"); // quiet, use minos
-–
+    rabs_gr->Draw("AP");
-–
+    rabs_gr->GetXaxis()->SetTitle("#alpha");
-–
+    rabs_gr->GetXaxis()->CenterTitle();
-–
+    rabs_gr->GetYaxis()->SetTitle("data/mc");
-–
+    rabs_gr->GetYaxis()->CenterTitle();
-<
+    TF1 *rabs_pol=(TF1*)rabs_gr->GetFunction("pol1");
-<
+    float rabs_a=rabs_pol->GetParameter(0);
-<
+    float rabs_b=rabs_pol->GetParameter(1);
-<
+    RABS_FinalGraph->SetPoint(ipt,0.5*(ptcuts[ipt]+ptcuts[ipt+1]),rabs_a);
-<
+    RABS_FinalGraph->SetPointError(ipt,0,rabs_pol->GetParError(0));
->
+    TH1D *mpf_band, *rabs_band;
-<
+    cout << "!+!+!+!+!+!+!+!+!+ Just added a point to the final plot : " << rabs_a << " +/- " << rabs_pol->GetParError(0) << endl;
->
+    if(!gofast) {
->
+      can->cd();
->
+      mpf_gr->SetMarkerStyle(21);
->
+      mpf_gr->SetMarkerColor(kBlue);
->
+      mpf_gr->Draw("AP");
->
->
+      TF1 *mpf_pol = new TF1("mpf_pol","[0]+[1]*x",0.0,0.4);
->
->
+      mpf_gr->Fit(mpf_pol); // better error estimation using MINOS (E), quiet (Q) and using loglikelihood method
->
+      mpf_band = getBand(mpf_pol,"MPF_band");
->
+      mpf_gr->GetXaxis()->SetTitle("#alpha");
->
+      mpf_gr->GetXaxis()->CenterTitle();
->
+      mpf_gr->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
->
+      mpf_gr->GetYaxis()->CenterTitle();
->
+      mpf_band->SetMarkerSize(0);
->
+      mpf_pol->SetLineWidth(0);
->
+      mpf_pol->SetLineColor(mpf_band->GetFillColor());
->
+      TF1 *mpf_pol_complete = new TF1("mpf_pol_complete","[0]+[1]*x");
->
+      mpf_pol_complete->SetParameters(mpf_pol->GetParameters());
->
->
+      mpf_pol_complete->SetLineWidth(1);
->
+      mpf_pol_complete->SetLineColor(kBlue);
->
->
+      float min,max;
->
+      FindMinMax(mpf_gr,min,max);
->
+      TH1F *histo = new TH1F("histo","histo",1,0,0.4);
->
+      histo->GetXaxis()->SetTitle("#alpha");
->
+      histo->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
->
+      histo->GetXaxis()->CenterTitle();
->
+      histo->GetYaxis()->CenterTitle();
->
+      histo->GetYaxis()->SetTitleOffset(1.3);
->
+      histo->SetStats(0);
->
+      histo->GetXaxis()->SetRangeUser(0,0.4);
->
+      histo->GetYaxis()->SetRangeUser(min,max);
->
+      mpf_gr->GetYaxis()->SetRangeUser(min,max);
->
->
+      //TPolyLine *mpf_fit_uncert = GetFitUncertaintyShape(mpf_gr, "pol1", min, max,0.0,0.4);
->
->
+      mpf_pol->SetLineColor(TColor::GetColor("#0101DF"));
->
+      //mpf_fit_uncert->SetFillColor(TColor::GetColor("#A2A2FA"));
->
+      mpf_band->SetFillColor(TColor::GetColor("#A2A2FA"));
->
+      histo->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
->
+      histo->Draw();
->
+      //mpf_fit_uncert->Draw("F");
->
+      //mpf_fit_uncert->Draw("");
->
+      mpf_gr->Draw("P");
->
+      mpf_band->Draw("C E3 SAME");
->
+      mpf_gr->Draw("P");
->
+      mpf_pol_complete->Draw("same");
->
->
+      float mpf_a=mpf_pol->GetParameter(0);
->
+      float mpf_b=mpf_pol->GetParameter(1);
->
+      MPF_FinalGraph->SetPoint(ipt,0.5*(ptcuts[ipt]+ptcuts[ipt+1]),mpf_pol->GetParameter(0));
->
+      MPF_FinalGraph->SetPointError(ipt,0,mpf_band->GetBinError(mpf_band->FindBin(0)));
->
->
+      MPF_ExtraPolatedResults[ipt]=mpf_a;
->
->
+      stringstream mpf_info;
->
+      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() << "}";
->
->
+      TText *mpf_mark = write_title(mpf_info.str());
->
+      mpf_mark->SetX(0.75);
->
+      mpf_mark->SetTextSize(0.03);
->
+      mpf_mark->SetY(0.75);
->
+      mpf_mark->Draw();
->
->
+      string filenamebkp=filename.str();
->
+      filename << "__MPF";
->
+      DrawPrelim();
->
+      CompleteSave(can,ContainerName+"/"+filename.str());
->
+      dout << "MPF : " << mpf_a << " + " << mpf_b << " * alpha " << endl;
->
->
+      filename.str("");
->
->
+      TF1 *rabs_pol = new TF1("rabs_pol","[0]+[1]*x",0.0,0.4);
->
+      rabs_gr->Fit(rabs_pol,"E");
->
+      rabs_band = getBand(rabs_pol,"Rabs_band");
->
+      rabs_pol->SetLineColor(rabs_band->GetFillColor());
->
+      rabs_gr->SetMarkerStyle(21);
->
+      rabs_gr->SetMarkerColor(kRed);
->
+      rabs_gr->Draw("AP");
->
+      //rabs_gr->Fit("pol1","E");
->
+      rabs_band->SetMarkerSize(0);
->
->
->
+      rabs_gr->GetXaxis()->SetTitle("#alpha");
->
+      rabs_gr->GetXaxis()->CenterTitle();
->
+      rabs_gr->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
->
+      rabs_gr->GetYaxis()->CenterTitle();
->
+      //TF1 *rabs_pol=(TF1*)rabs_gr->GetFunction("pol1");
->
+      rabs_pol->SetLineWidth(0);
->
+      TF1 *rabs_pol_complete = new TF1("rabs_pol_complete","[0]+[1]*x");
->
+      rabs_pol_complete->SetParameters(rabs_pol->GetParameters());
->
+      rabs_pol_complete->SetLineColor(kRed);
->
+      rabs_pol_complete->SetLineWidth(1);
->
+      FindMinMax(rabs_gr,min,max);
->
+      rabs_gr->GetYaxis()->SetRangeUser(min,max);
->
+      histo->GetYaxis()->SetRangeUser(min,max);
->
+      //TPolyLine *rabs_fit_uncert = GetFitUncertaintyShape(rabs_gr, "pol1", min, max,0.0,0.4);
->
+      rabs_pol->SetLineColor(TColor::GetColor("#F99C9C"));
->
+      rabs_pol->SetFillColor(TColor::GetColor("#F99C9C"));
->
+//      rabs_band->SetLineColor(TColor::GetColor("#F99C9C"));
->
+      rabs_band->SetFillColor(TColor::GetColor("#F99C9C"));
->
+      //rabs_fit_uncert->SetLineColor(TColor::GetColor("#F99C9C"));
->
+      //rabs_fit_uncert->SetFillColor(TColor::GetColor("#F99C9C"));
->
+      histo->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
->
+      histo->Draw();
->
+      rabs_band->Draw("C E3 SAME");
->
+      //rabs_fit_uncert->Draw("F");
->
+      //rabs_fit_uncert->Draw("");
->
+      rabs_gr->Draw("P");
->
+      rabs_pol_complete->Draw("same");
->
->
+      float rabs_a=rabs_pol->GetParameter(0);
->
+      float rabs_b=rabs_pol->GetParameter(1);
->
->
+      RABS_FinalGraph->SetPoint(ipt,0.5*(ptcuts[ipt]+ptcuts[ipt+1]),rabs_a);
->
+      RABS_FinalGraph->SetPointError(ipt,0,rabs_band->GetBinError(rabs_band->FindBin(0)));
->
->
+      dout << "!+!+!+!+!+!+!+!+!+ Just added a point to the final plot : " << rabs_a << " +/- " << rabs_band->GetBinError(rabs_band->FindBin(0)) << endl;
->
->
+      stringstream rabs_info;
->
+      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() << "}";
->
+      TText *rabs_mark = write_title(rabs_info.str());
->
+      rabs_mark->SetX(0.75);
->
+      rabs_mark->SetTextSize(0.03);
->
+      rabs_mark->SetY(0.75);
->
+      rabs_mark->Draw();
->
->
+      RABS_ExtraPolatedResults[ipt]=rabs_a;
->
+      filename << filenamebkp << "__RABS";
->
+      DrawPrelim();
->
+      CompleteSave(can,ContainerName+"/"+filename.str());
->
+      dout << "RABS : " << rabs_a << " + " << rabs_b << " * alpha " << endl;
->
->
+      delete mpf_band;
->
+      delete rabs_band;
->
+    }
->
+  }
->
->
+  float MPFResult;
->
+  float MPFResultError;
->
->
+  float RabsResult;
->
+  float RabsResultError;
->
->
+  if(!gofast) {
->
+    MPF_FinalGraph->SetMarkerStyle(21);
->
+    MPF_FinalGraph->SetMarkerColor(kBlue);
->
+    MPF_FinalGraph->Fit("pol0",""); // quiet, use minos
-+
+    TF1 *mpf_pol0=(TF1*)MPF_FinalGraph->GetFunction("pol0");
-+
+    TF1 *mpf_pol0_complete = new TF1("mpf_pol0_complete","[0]+[1]*x");
-+
+    mpf_pol0_complete->SetLineWidth(1);
-+
+    mpf_pol0_complete->SetLineColor(kBlue);
-+
+    mpf_pol0->SetLineColor(kBlue);
-+
+    MPF_FinalGraph->Draw("AP");
-+
+    MPF_FinalGraph->GetXaxis()->SetTitle("p_{T}^{Z}");
-+
+    MPF_FinalGraph->GetXaxis()->CenterTitle();
-+
+    MPF_FinalGraph->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
-+
+    MPF_FinalGraph->GetYaxis()->CenterTitle();
-+
+    MPF_FinalGraph->Draw("AP");
-+
+    mpf_pol0_complete->Draw("same");
-<
+    stringstream rabs_info;
-<
+    rabs_info << "#splitline{#splitline{#splitline{" << ptcuts[ipt] << " GeV < p_{T}^{Z} < " << ptcuts[ipt+1] << " GeV }{ (R_{abs}) }}{Extrapolated value: " << std::setprecision(3) << rabs_a << "}}{#chi^{2}/NDF : " << rabs_pol->GetChisquare() << " / " << rabs_pol->GetNDF() << "}";
-<
+    TText *rabs_mark = write_title(rabs_info.str());
-<
+    rabs_mark->SetX(0.75);
-<
+    rabs_mark->SetTextSize(0.03);
-<
+    rabs_mark->SetY(0.75);
-<
+    rabs_mark->Draw();
->
+    stringstream mpf_result;
->
+    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)}";
->
+    TText *rmpf_final_mark = write_title(mpf_result.str());
->
+    rmpf_final_mark->SetX(0.75);
->
+    rmpf_final_mark->SetY(0.75);
->
+    rmpf_final_mark->SetTextSize(0.03);
->
+    rmpf_final_mark->Draw();
->
+    DrawPrelim();
->
+    MPFResult=1/mpf_pol0->GetParameter(0);
->
+    MPFResultError=mpf_pol0->GetParError(0)/(mpf_pol0->GetParameter(0)*mpf_pol0->GetParameter(0));
-+
+    stringstream filename2;
-+
+    filename2 << "Extrapolation/FINAL_RESULT_MPF";
-+
+    CompleteSave(can,ContainerName+"/"+filename2.str());
-<
-<
+    RABS_ExtraPolatedResults[ipt]=rabs_a;
-<
+    filename << filenamebkp << "__RABS";
->
+    RABS_FinalGraph->SetMarkerStyle(21);
->
+    RABS_FinalGraph->SetMarkerStyle(21);
->
+    RABS_FinalGraph->SetMarkerColor(kRed);
->
+    RABS_FinalGraph->Fit("pol0","QE"); // quiet, use minos
->
+    RABS_FinalGraph->Draw("AP");
->
+    RABS_FinalGraph->GetXaxis()->SetTitle("p_{T}^{Z}");
->
+    RABS_FinalGraph->GetXaxis()->CenterTitle();
->
+    RABS_FinalGraph->GetYaxis()->SetTitle("<R_{abs}>_{data}/<R_{abs}>_{mc}");
->
+    RABS_FinalGraph->GetYaxis()->CenterTitle();
->
+    RABS_FinalGraph->Draw("AP");
->
+    TF1 *rabs_pol0=(TF1*)RABS_FinalGraph->GetFunction("pol0");
->
+    stringstream rabs_result;
->
+    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})}";;
->
+    TText *rabs_final_mark = write_title(rabs_result.str());
->
+    rabs_final_mark->SetX(0.75);
->
+    rabs_final_mark->SetY(0.75);
->
+    rabs_final_mark->SetTextSize(0.03);
->
+    rabs_final_mark->Draw();
+    DrawPrelim();
-<
+    CompleteSave(can,filename.str());
-<
+    cout << "RABS : " << rabs_a << " + " << rabs_b << " * alpha " << endl;
->
->
+    RabsResult=1/rabs_pol0->GetParameter(0);
->
+    RabsResultError=rabs_pol0->GetParError(0)/(rabs_pol0->GetParameter(0)*rabs_pol0->GetParameter(0));
->
->
+    filename2.str("");
->
+    filename2 << "Extrapolation/FINAL_RESULT_RABS";
->
+    CompleteSave(can,ContainerName+"/"+filename2.str());
-<
-<
+  MPF_FinalGraph->SetMarkerStyle(21);
-<
+  MPF_FinalGraph->SetMarkerColor(kRed);
-<
+  MPF_FinalGraph->Fit("pol0","QE"); // quiet, use minos
-<
-<
+  TF1 *mpf_pol0=(TF1*)MPF_FinalGraph->GetFunction("pol0");
-<
+  MPF_FinalGraph->Draw("AP");
-<
+  MPF_FinalGraph->GetXaxis()->SetTitle("p_{T}^{Z}");
-<
+  MPF_FinalGraph->GetXaxis()->CenterTitle();
-<
+  MPF_FinalGraph->GetYaxis()->SetTitle("data/mc");
-<
+  MPF_FinalGraph->GetYaxis()->CenterTitle();
-<
+  MPF_FinalGraph->Draw("AP");
-<
-<
+  stringstream mpf_result;
-<
+  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)}";
-<
+  TText *rmpf_final_mark = write_title(mpf_result.str());
-<
+  rmpf_final_mark->SetX(0.75);
-<
+  rmpf_final_mark->SetY(0.75);
-<
+  rmpf_final_mark->SetTextSize(0.03);
-<
+  rmpf_final_mark->Draw();
->
->
+  can->cd();
->
+  MPF_FaceValueAtPoint3->SetMarkerStyle(21);
->
+  MPF_FaceValueAtPoint3->SetMarkerColor(kBlue);
->
+  MPF_FaceValueAtPoint3->Fit("pol0","QE"); // quiet, use minos
->
+  MPF_FaceValueAtPoint3->Draw("AP");
->
+  MPF_FaceValueAtPoint3->GetXaxis()->SetTitle("p_{T}^{Z}");
->
+  MPF_FaceValueAtPoint3->GetXaxis()->CenterTitle();
->
+  MPF_FaceValueAtPoint3->GetYaxis()->SetTitle(((string)"<MPF>__{data}/<MPF>_{mc} for #alpha<"+any2string(FaceValueToBeTakenAt)).c_str());
->
+  MPF_FaceValueAtPoint3->GetYaxis()->CenterTitle();
->
+  MPF_FaceValueAtPoint3->Draw("AP");
->
+  TF1 *faceval_pol0=(TF1*)MPF_FaceValueAtPoint3->GetFunction("pol0");
->
+  faceval_pol0->SetLineColor(kBlue);
->
+  faceval_pol0->SetLineWidth(1);
->
+  faceval_pol0->Draw("same");
->
+  TF1 *faceval_pol0_complete = new TF1("faceval_pol0_complete","[0]+[1]*x");
->
+  faceval_pol0_complete->SetParameters(faceval_pol0->GetParameters());
->
+  stringstream faceval_result;
->
+  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 << ")}";;
->
+  TText *faceval_final_mark = write_title(faceval_result.str());
->
+  faceval_final_mark->SetX(0.75);
->
+  faceval_final_mark->SetY(0.75);
->
+  faceval_final_mark->SetTextSize(0.03);
->
+  faceval_final_mark->Draw();
+  DrawPrelim();
-<
+  stringstream filename2;
-<
+  filename2 << "Extrapolation/FINAL_RESULT_MPF";
-<
+  CompleteSave(can,filename2.str());
->
+  float FaceValResult=1/faceval_pol0->GetParameter(0);
->
+  float FaceValResultError=faceval_pol0->GetParError(0)/(faceval_pol0->GetParameter(0)*faceval_pol0->GetParameter(0));
-+
+  stringstream filename2;
-+
+  filename2.str("");
-+
+  filename2 << "Extrapolation/FINAL_RESULT_MPF_FaceValp3";
-+
+  CompleteSave(can,ContainerName+"/"+filename2.str());
-<
+  RABS_FinalGraph->SetMarkerStyle(21);
-<
+  RABS_FinalGraph->SetMarkerStyle(21);
-<
+  RABS_FinalGraph->SetMarkerColor(kRed);
-<
+  RABS_FinalGraph->Fit("pol0","QE"); // quiet, use minos
-<
+  RABS_FinalGraph->Draw("AP");
-<
+  RABS_FinalGraph->GetXaxis()->SetTitle("p_{T}^{Z}");
-<
+  RABS_FinalGraph->GetXaxis()->CenterTitle();
-<
+  RABS_FinalGraph->GetYaxis()->SetTitle("data/mc");
-<
+  RABS_FinalGraph->GetYaxis()->CenterTitle();
-<
+  RABS_FinalGraph->Draw("AP");
-<
+  TF1 *rabs_pol0=(TF1*)RABS_FinalGraph->GetFunction("pol0");
-<
+  stringstream rabs_result;
-<
+  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})}";;
-<
+  TText *rabs_final_mark = write_title(rabs_result.str());
-<
+  rabs_final_mark->SetX(0.75);
-<
+  rabs_final_mark->SetY(0.75);
-<
+  rabs_final_mark->SetTextSize(0.03);
-<
+  rabs_final_mark->Draw();
->
+  can->cd();
->
+  RABS_FaceValueAtPoint3->SetMarkerStyle(21);
->
+  RABS_FaceValueAtPoint3->SetMarkerColor(kRed);
->
+  RABS_FaceValueAtPoint3->Fit("pol0","QE"); // quiet, use minos
->
+  RABS_FaceValueAtPoint3->Draw("AP");
->
+  RABS_FaceValueAtPoint3->GetXaxis()->SetTitle("p_{T}^{Z}");
->
+  RABS_FaceValueAtPoint3->GetXaxis()->CenterTitle();
->
+  RABS_FaceValueAtPoint3->GetYaxis()->SetTitle(((string)"<R_{abs}>_{data}/<R_{abs}>_{mc} for #alpha<"+any2string(FaceValueToBeTakenAt)).c_str());
->
+  RABS_FaceValueAtPoint3->GetYaxis()->CenterTitle();
->
+  RABS_FaceValueAtPoint3->Draw("AP");
->
+  TF1 *faceval_pol0RABS=(TF1*)RABS_FaceValueAtPoint3->GetFunction("pol0");
->
+  TF1 *faceval_pol0RABS_complete = new TF1("faceval_pol0_complete","[0]+[1]*x");
->
+  faceval_pol0RABS_complete->SetParameters(faceval_pol0RABS->GetParameters());
->
+  faceval_pol0RABS->SetLineWidth(1);
->
+  faceval_pol0RABS->SetLineColor(kRed);
->
+  faceval_pol0RABS_complete->Draw("same");
->
+  stringstream RABSfaceval_result;
->
+  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 << ")}";;
->
+  TText *RABSfaceval_final_mark = write_title(RABSfaceval_result.str());
->
+  RABSfaceval_final_mark->SetX(0.75);
->
+  RABSfaceval_final_mark->SetY(0.75);
->
+  RABSfaceval_final_mark->SetTextSize(0.03);
->
+  RABSfaceval_final_mark->Draw();
+  DrawPrelim();
-+
+  float RABSFaceValResult=1/faceval_pol0RABS->GetParameter(0);
-+
+  float RABSFaceValResultError=faceval_pol0RABS->GetParError(0)/(faceval_pol0RABS->GetParameter(0)*faceval_pol0RABS->GetParameter(0));
-+
+  filename2.str("");
-<
+  filename2 << "Extrapolation/FINAL_RESULT_RABS";
-<
+  CompleteSave(can,filename2.str());
->
+  filename2 << "Extrapolation/FINAL_RESULT_RABS_FaceValp3";
->
+  CompleteSave(can,ContainerName+"/"+filename2.str());
-+
+  dout << "FINAL RESULTS: " << endl;
-+
+  if(!gofast) dout << "MPF : " << MPFResult  << " +/- " << MPFResultError  << endl;
-+
+  if(!gofast) dout << "Rabs: " << RabsResult << " +/- " << RabsResultError << endl;
-+
+  dout << "Face value at " << FaceValueToBeTakenAt << " :" << FaceValResult << " +/- " << FaceValResultError << "  (MPF) " << endl;
-+
+  dout << "Face value at " << FaceValueToBeTakenAt << " :" << RABSFaceValResult << " +/- " << RABSFaceValResultError << "  (RABS) " << endl;
+  delete can;
-+
-+
+  for(unsigned int i=0;i<MPF_FinalGraph->GetN();i++) {
-+
+    double x,y,dy;
-+
+    MPF_FinalGraph->GetPoint(i,x,y);
-+
+    dy = MPF_FinalGraph->GetErrorY(i);
-+
+    results.MPF_Result_Extrapolated_Bins.push_back(Value(y,dy));
-+
+  }
-+
-+
+  for(unsigned int i=0;i<MPF_FaceValueAtPoint3->GetN();i++) {
-+
+    double x,y,dy;
-+
+    MPF_FaceValueAtPoint3->GetPoint(i,x,y);
-+
+    dy = MPF_FaceValueAtPoint3->GetErrorY(i);
-+
+    results.MPF_Result_FaceValue_Bins.push_back(Value(y,dy));
-+
+  }
-+
-+
+  for(unsigned int i=0;i<RABS_FinalGraph->GetN();i++) {
-+
+    double x,y,dy;
-+
+    RABS_FinalGraph->GetPoint(i,x,y);
-+
+    dy = RABS_FinalGraph->GetErrorY(i);
-+
+    results.Rabs_Result_Extrapolated_Bins.push_back(Value(y,dy));
-+
+  }
-+
-+
+  for(unsigned int i=0;i<RABS_FaceValueAtPoint3->GetN();i++) {
-+
+    double x,y,dy;
-+
+    RABS_FaceValueAtPoint3->GetPoint(i,x,y);
-+
+    dy = RABS_FaceValueAtPoint3->GetErrorY(i);
-+
+    results.Rabs_Result_FaceValue_Bins.push_back(Value(y,dy));
-+
+  }
-+
-+
-+
+  results.MPF_Result_FaceValue=Value(FaceValResult,FaceValResultError);
-+
+  results.Rabs_Result_FaceValue=Value(RABSFaceValResult,RABSFaceValResultError);
-+
+  if(!gofast) results.MPF_Result_Extrapolated=Value(MPFResult,MPFResultError);
-+
+  if(!gofast) results.Rabs_Result_Extrapolated=Value(RabsResult,RabsResultError);
-+
+  results.Print();
-+
+  return results;
-+
+void DoPUStudy(string identifier) {
-+
-+
+  float numVtxcuts[7]={0,10,15,20};
-+
-+
+  TCanvas *can = new TCanvas("can","can");
-+
+  TH1F *malpha[8];
-+
+  TH1F *dalpha[8];
-+
+  dout << identifier << ";";
-+
+  dout << endl;
-+
+  for(int i=1;i<5;i++) {
-+
+    stringstream sSpecialCut;
-+
+    if(i<4) {
-+
+      sSpecialCut << "numVtx>" << numVtxcuts[i-1] << " && numVtx<" << numVtxcuts[i];
-+
+      dout << numVtxcuts[i-1] << "  < nVtx < " << numVtxcuts[i] << ";";
-+
+    } else {
-+
+      sSpecialCut << "numVtx>" << numVtxcuts[i-1];
-+
+      dout << numVtxcuts[i-1] << ";";
-+
+    }
-+
-+
+    sSpecialCut << " && " << identifier << "_alpha<0.3";
-+
+    TCut SpecialCut(sSpecialCut.str().c_str());
-+
-+
-+
+    malpha[i] = allsamples.Draw("malpha"+any2string(i),"mpf",1000,0,2,"MPF","events",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&SpecialCut,mc,luminosity);
-+
+    dalpha[i] = allsamples.Draw("dalpha"+any2string(i),"mpf",1000,0,2,"MPF","events",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&SpecialCut,data,luminosity);
-+
-+
-+
+    float a=dalpha[i]->GetMean();
-+
+    float b=malpha[i]->GetMean();
-+
+    float da=dalpha[i]->GetMeanError();
-+
+    float db=malpha[i]->GetMeanError();
-+
+    float factor = a / b;
-+
+    float error = TMath::Sqrt(  (1/(b*b)) * (da*da) + ((a*a)/(b*b))*db*db);
-+
-+
+    dout << malpha[i]->GetMean() << ";" << malpha[i]->GetMeanError() << ";" << dalpha[i]->GetMean() << ";" << dalpha[i]->GetMeanError() << ";" << malpha[i]->Integral()<<";"<<dalpha[i]->Integral() << ";"<<factor << ";" << error <<  endl;
-+
+    delete malpha[i];
-+
+    delete dalpha[i];
-+
+  }
-+
-+
+  delete can;
-+
-+
+}
-+
-+
+void ScenarioComparison() {
-+
+//very dumb way to do this.
-+
-+
+  TGraphAsymmErrors *gr[5];
-+
+  TF1 *fit[5];
-+
+  bool dofit=true;
-+
-+
+  TCanvas *can = new TCanvas("can","can",500,500);
-+
-+
+  TLegend *leg = new TLegend(0.2,0.2,0.5,0.45);
-+
+  leg->SetBorderSize(0);
-+
+  leg->SetFillColor(kWhite);
-+
-+
+  for(int i=0;i<6;i++) {
-+
+    gr[i] = new TGraphAsymmErrors();
-+
+    double x[4]   = {5,12.5,17.5,25};
-+
+    double dxl[4] = {5,2.5,2.5,5};
-+
+    double dxh[4] = {5,2.5,2.5,10};
-+
-+
+    string name="";
-+
+    string color="";
-+
-+
+    if(i==0) {
-+
+      double y[4]   = {1.01935,1.01558,0.996438,1.01901};
-+
+      double dyl[4] = {0.0132973,0.0131661,0.017629,0.0319767};
-+
+      double dyh[4] = {0.0132973,0.0131661,0.017629,0.0319767};
-+
+      name="30/30";
-+
+      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
-+
+      color="#a01d99";
-+
+    }
-+
+    if(i==1) {
-+
+      double y[4]   = {1.00577,1.00706,0.991029,0.954098};
-+
+      double dyl[4] = {0.0135575,0.0136548,0.0170778,0.0263919};
-+
+      double dyh[4] = {0.0135575,0.0136548,0.0170778,0.0263919};
-+
+      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
-+
+      name="30/15";
-+
+      color="#2464bc";
-+
+    }
-+
+    if(i==2) {
-+
+      double y[4]   = {1.023,1.00648,1.02217,1.03181};
-+
+      double dyl[4] = {0.0157372,0.0187778,0.0261512,0.0422295};
-+
+      double dyh[4] = {0.0157372,0.0187778,0.0261512,0.0422295};
-+
+      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
-+
+      name="15/15";
-+
+      color="#f49230";
-+
+    }
-+
+    if(i==3) {
-+
+      double y[4]   = {1.02029,1.01338,0.999769,1.00357};
-+
+      double dyl[4] = {0.0129509,0.0120127,0.0160537,0.0268972};
-+
+      double dyh[4] = {0.0129509,0.0120127,0.0160537,0.0268972};
-+
+      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
-+
+      name="CHS 30/30";
-+
+      color="#c72f3c";
-+
+    }
-+
+    if(i==4) {
-+
+      double y[4]   = {1.00274,1.01056,0.993214,0.997963};
-+
+      double dyl[4] = {0.0127179,0.0140344,0.0162814,0.0278914};
-+
+      double dyh[4] = {0.0127179,0.0140344,0.0162814,0.0278914};
-+
+      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
-+
+      name="CHS 30/15";
-+
+      color="#7cb433";
-+
+    }
-+
+    if(i==5) {
-+
+      double y[4]   = {1.04103,1.03726,1.04649,1.05383};
-+
+      double dyl[4] = {0.0156465,0.0168194,0.0231982,0.0376114};
-+
+      double dyh[4] = {0.0156465,0.0168194,0.0231982,0.0376114};
-+
+      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
-+
+      name="CHS 15/15";
-+
+      color="#DF013A";
-+
+    }
-+
-+
-+
+    gStyle->SetOptFit(0);
-+
+    gr[i]->SetTitle();
-+
+    gr[i]->GetXaxis()->SetTitle("N(Vtx)");
-+
+    gr[i]->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
-+
+    gr[i]->GetXaxis()->CenterTitle();
-+
+    gr[i]->GetYaxis()->CenterTitle();
-+
+    gr[i]->SetName(name.c_str());
-+
+    gr[i]->SetMarkerSize(0.00001);
-+
+    gr[i]->GetYaxis()->SetRangeUser(0.90,1.10);
-+
+    gr[i]->SetLineColor(TColor::GetColor(color.c_str()));
-+
+    gr[i]->SetMarkerColor(TColor::GetColor(color.c_str()));
-+
+    if(dofit) {
-+
+      gr[i]->Fit("pol1");
-+
+      fit[i] = (TF1*)gr[i]->GetFunction("pol1");
-+
+      fit[i]->SetLineColor(TColor::GetColor(color.c_str()));
-+
+    }
-+
+    leg->AddEntry(gr[i],name.c_str(),"l");
-+
+    if(i==0) gr[i]->Draw("AP*");
-+
+    else gr[i]->Draw("P");
-+
+  }
-+
+  leg->Draw();
-+
+  DrawPrelim();
-+
-+
+  CompleteSave(can,"ScenarioComparison");
-+
+  delete can;
-+
-+
+}
-+
-+
+void ScenarioComparisonInclusive() {
-+
+  //dumbest way ever to do this. but ok we only need to do it once.
-+
+  TGraphAsymmErrors *gr[5];
-+
+  TF1 *fit[5];
-+
+  bool dofit=true;
-+
-+
+  TCanvas *can = new TCanvas("can","can",500,500);
-+
-+
+  TLegend *leg = new TLegend(0.2,0.2,0.65,0.45);
-+
+  leg->SetBorderSize(0);
-+
+  leg->SetFillColor(kWhite);
-+
-+
+  for(int i=0;i<6;i++) {
-+
+    gr[i] = new TGraphAsymmErrors();
-+
+    double x[4]   = {5,12.5,17.5,25};
-+
+    double dxl[4] = {5,2.5,2.5,5};
-+
+    double dxh[4] = {5,2.5,2.5,5};
-+
-+
+    string name="";
-+
+    string color="";
-+
-+
+    if(i==0) {
-+
+      double y[4]   = {1.01093,1.01379,1.0147,1.02766};
-+
+      double dyl[4] = {0.00243596,0.00242534,0.00341941,0.00600663};
-+
+      double dyh[4] = {0.00243596,0.00242534,0.00341941,0.00600663};
-+
+      name="30/30";
-+
+      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
-+
+      color="#a01d99";
-+
+    }
-+
+    if(i==1) {
-+
+      double y[4]   = {1.0068,1.00779,0.999692,0.994484};
-+
+      double dyl[4] = {0.00249233,0.00257169,0.00371399,0.00631138};
-+
+      double dyh[4] = {0.00249233,0.00257169,0.00371399,0.00631138};
-+
+      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
-+
+      name="30/15";
-+
+      color="#2464bc";
-+
+    }
-+
+    if(i==2) {
-+
+      double y[4]   = {1.01208,1.01991,1.02185,1.03479};
-+
+      double dyl[4] = {0.00366745,0.00400001,0.00584054,0.0104661};
-+
+      double dyh[4] = {0.00366745,0.00400001,0.00584054,0.0104661};
-+
+      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
-+
+      name="15/15";
-+
+      color="#f49230";
-+
+    }
-+
+    if(i==3) {
-+
+      double y[4]   = {1.01177,1.01568,1.01616,1.03211};
-+
+      double dyl[4] = {0.00217317,0.00203014,0.00276838,0.00465643};
-+
+      double dyh[4] = {0.00217317,0.00203014,0.00276838,0.00465643};
-+
+      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
-+
+      name="CHS 30/30";
-+
+      color="#c72f3c";
-+
+    }
-+
+    if(i==4) {
-+
+      double y[4]   = {1.00751,1.00964,1.00804,1.01073};
-+
+      double dyl[4] = {0.00224851,0.00215563,0.0029592,0.00499805};
-+
+      double dyh[4] = {0.00224851,0.00215563,0.0029592,0.00499805};
-+
+      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
-+
+      name="CHS 30/15";
-+
+      color="#7cb433";
-+
+    }
-+
+    if(i==5) {
-+
+      double y[4]   = {1.01557,1.0318,1.03535,1.06331};
-+
+      double dyl[4] = {0.00300783,0.00302486,0.00421064,0.00739416};
-+
+      double dyh[4] = {0.00300783,0.00302486,0.00421064,0.00739416};
-+
+      gr[i] = new TGraphAsymmErrors(4,x,y,dxl,dxh,dyl,dyh);
-+
+      name="CHS 15/15";
-+
+      color="#DF013A";
-+
+    }
-+
-+
+    gStyle->SetOptFit(0);
-+
+    gr[i]->SetTitle();
-+
+    gr[i]->GetXaxis()->SetTitle("N(Vtx)");
-+
+    gr[i]->GetYaxis()->SetTitle("<MPF>_{data}/<MPF>_{mc}");
-+
+    gr[i]->GetXaxis()->CenterTitle();
-+
+    gr[i]->GetYaxis()->CenterTitle();
-+
+    gr[i]->SetName(name.c_str());
-+
+    gr[i]->SetMarkerSize(0.00001);
-+
+    gr[i]->GetYaxis()->SetRangeUser(0.90,1.10);
-+
+    gr[i]->SetLineColor(TColor::GetColor(color.c_str()));
-+
+    gr[i]->SetMarkerColor(TColor::GetColor(color.c_str()));
-+
+    if(dofit) {
-+
+      gr[i]->Fit("pol1");
-+
+      fit[i] = (TF1*)gr[i]->GetFunction("pol1");
-+
+      fit[i]->SetLineColor(TColor::GetColor(color.c_str()));
-+
+    }
-+
-+
+    stringstream nameplus;
-+
+    nameplus << name << " (slope: " << std::setprecision(3) << fit[i]->GetParameter(1) << ")";
-+
-+
+    leg->AddEntry(gr[i],nameplus.str().c_str(),"l");
-+
+    if(i==0) gr[i]->Draw("AP*");
-+
+    else gr[i]->Draw("P");
-+
+  }
-+
+  leg->Draw();
-+
+  DrawPrelim();
-+
-+
+  CompleteSave(can,"ScenarioComparisonInclusive");
-+
+  delete can;
-+
-+
+}
-+
-+
+void compare_selection(string identifier) {
-+
+  bool recognized_scenario=false;
-+
-+
+  dout << "Running with identifier " << identifier << endl;
-+
+  if(identifier=="Zb30") {
-+
+      LeadingB=TCut ("Zb30_bTagProbCSVBP[0]>0.898");
-+
+      EtaB=TCut("abs(Zb30_pfJetGoodEta[0])<1.3");
-+
+      PhiZcut=TCut("abs(Zb30_pfJetDphiZ[0])>2.7");
-+
+      recognized_scenario=true;
-+
+  }
-+
+  if(identifier=="Zb3010") {
-+
+      LeadingB=TCut ("Zb3010_bTagProbCSVBP[0]>0.898");
-+
+      EtaB=TCut("abs(Zb3010_pfJetGoodEta[0])<1.3");
-+
+      PhiZcut=TCut("abs(Zb3010_pfJetDphiZ[0])>2.7");
-+
+      recognized_scenario=true;
-+
+  }
-+
-+
+  if(identifier=="ZbCHS3010") {
-+
+      LeadingB=TCut ("ZbCHS1510_bTagProbCSVBP[0]>0.898");
-+
+      EtaB=TCut("abs(ZbCHS1510_pfJetGoodEta[0])<1.3");
-+
+      PhiZcut=TCut("abs(ZbCHS1510_pfJetDphiZ[0])>2.7");
-+
+      recognized_scenario=true;
-+
+  }
-+
-+
+  if(identifier=="Zb40") {
-+
+      LeadingB=TCut ("Zb40_bTagProbCSVBP[0]>0.898");
-+
+      EtaB=TCut("abs(Zb40_pfJetGoodEta[0])<1.3");
-+
+      PhiZcut=TCut("abs(Zb40_pfJetDphiZ[0])>2.7");
-+
+      recognized_scenario=true;
-+
+  }
-+
-+
+  if(identifier=="ZbCHS301030") {
-+
+      LeadingB=TCut ("ZbCHS301030_bTagProbCSVBP[0]>0.898");
-+
+      EtaB=TCut("abs(ZbCHS301030_pfJetGoodEta[0])<1.3");
-+
+      PhiZcut=TCut("abs(ZbCHS301030_pfJetDphiZ[0])>2.7");
-+
+      recognized_scenario=true;
-+
+  }
-+
-+
+  if(identifier=="Zb1530") {
-+
+      LeadingB=TCut ("Zb1530_bTagProbCSVBP[0]>0.898");
-+
+      EtaB=TCut("abs(Zb1530_pfJetGoodEta[0])<1.3");
-+
+      PhiZcut=TCut("abs(Zb1530_pfJetDphiZ[0])>2.7");
-+
+      recognized_scenario=true;
-+
+  }
-+
-+
+  if(identifier=="ZbCHS301010") {
-+
+      LeadingB=TCut ("ZbCHS301010_bTagProbCSVBP[0]>0.898");
-+
+      EtaB=TCut("abs(ZbCHS301010_pfJetGoodEta[0])<1.3");
-+
+      PhiZcut=TCut("abs(ZbCHS301010_pfJetDphiZ[0])>2.7");
-+
+      recognized_scenario=true;
-+
+  }
-+
-+
+  if(identifier=="Zb1510") {
-+
+      LeadingB=TCut ("Zb1510_bTagProbCSVBP[0]>0.898");
-+
+      EtaB=TCut("abs(Zb1510_pfJetGoodEta[0])<1.3");
-+
+      PhiZcut=TCut("abs(Zb1510_pfJetDphiZ[0])>2.7");
-+
+      recognized_scenario=true;
-+
+  }
-+
-+
+  if(identifier=="ZbCHS301010") {
-+
+      LeadingB=TCut ("ZbCHS301010_bTagProbCSVBP[0]>0.898");
-+
+      EtaB=TCut("abs(ZbCHS301010_pfJetGoodEta[0])<1.3");
-+
+      PhiZcut=TCut("abs(ZbCHS301010_pfJetDphiZ[0])>2.7");
-+
+      recognized_scenario=true;
-+
+  }
-+
-+
+  if(identifier=="ZbMikko") {
-+
+      LeadingB=TCut ("ZbMikko_bTagProbCSVBP[0]>0.898 && ZbMikko__IsClean");
-+
+      EtaB=TCut("abs(ZbMikko_pfJetGoodEta[0])<1.3 && ZbMikko__IsClean");
-+
+      PhiZcut=TCut("abs(ZbMikko_pfJetDphiZ[0])>2.7 && ZbMikko__IsClean");
-+
+      recognized_scenario=true;
-+
+  }
-+
-+
+  if(identifier=="ZbCHS3010_SecEta3") {
-+
+      LeadingB=TCut ("ZbCHS3010_SecEta3_bTagProbCSVBP[0]>0.898");
-+
+      EtaB=TCut("abs(ZbCHS3010_SecEta3_pfJetGoodEta[0])<1.3");
-+
+      PhiZcut=TCut("abs(ZbCHS3010_SecEta3_pfJetDphiZ[0])>2.7");
-+
+      recognized_scenario=true;
-+
+  }
-+
-+
+  if(identifier=="ZbCHS3010_SecEta5") {
-+
+      LeadingB=TCut ("ZbCHS3010_SecEta5_bTagProbCSVBP[0]>0.898");
-+
+      EtaB=TCut("abs(ZbCHS3010_SecEta5_pfJetGoodEta[0])<1.3");
-+
+      PhiZcut=TCut("abs(ZbCHS3010_SecEta5_pfJetDphiZ[0])>2.7");
-+
+      recognized_scenario=true;
-+
+  }
-+
-+
+  if(identifier=="ZbCHS3010_p5Clean") {
-+
+      LeadingB=TCut ("ZbCHS3010_p5Clean_bTagProbCSVBP[0]>0.898 && ZbCHS3010_p5Clean_IsClean");
-+
+      EtaB=TCut("abs(ZbCHS3010_p5Clean_pfJetGoodEta[0])<1.3 && ZbCHS3010_p5Clean_IsClean");
-+
+      PhiZcut=TCut("abs(ZbCHS3010_p5Clean_pfJetDphiZ[0])>2.7 && ZbCHS3010_p5Clean_IsClean");
-+
+      recognized_scenario=true;
-+
+  }
-+
-+
+  if(identifier=="ZbCHS3010") {
-+
+      LeadingB=TCut ("ZbCHS3010_bTagProbCSVBP[0]>0.898");
-+
+      EtaB=TCut("abs(ZbCHS3010_pfJetGoodEta[0])<1.3");
-+
+      PhiZcut=TCut("abs(ZbCHS3010_pfJetDphiZ[0])>2.7");
-+
+      recognized_scenario=true;
-+
+  }
-+
+  if(identifier=="ZbCHS3010") {
-+
+      LeadingB=TCut ("ZbCHS3010_bTagProbCSVBP[0]>0.898");
-+
+      EtaB=TCut("abs(ZbCHS3010_pfJetGoodEta[0])<1.3");
-+
+      PhiZcut=TCut("abs(ZbCHS3010_pfJetDphiZ[0])>2.7");
-+
+      recognized_scenario=true;
-+
+  }
-+
+  if(identifier=="ZbCHS30") {
-+
+      LeadingB=TCut ("ZbCHS30_bTagProbCSVBP[0]>0.898");
-+
+      EtaB=TCut("abs(ZbCHS30_pfJetGoodEta[0])<1.3");
-+
+      PhiZcut=TCut("abs(ZbCHS30_pfJetDphiZ[0])>2.7");
-+
+      recognized_scenario=true;
-+
+  }
-+
-+
+  assert(recognized_scenario);
-+
-+
+  dout << "Cuts: " << endl;
-+
+  dout << "   " << (const char*) LeadingB << endl;
-+
+  dout << "   " << (const char*) EtaB << endl;
-+
+  dout << "   " << (const char*) PhiZcut << endl;
-+
-+
+  dout << endl << endl;
-+
-+
+//  ScenarioComparison();
-+
+//  ScenarioComparisonInclusive();
-+
-+
+//    DoPUStudy(identifier);
-+
+//    SpecialCutFlow(identifier);
-+
-+
+//    draw_kin_variable("pt",ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&"pt>30&&pt<1000",17000,30,200,"Z p_{T}","Official/"+identifier+"/Zpt__"+identifier,1);
-+
+//    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);
-+
+//    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);
-+
+//    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);
-+
+  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);
-+
+  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);
-+
+  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);
-+
+  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);
-+
+  draw_kin_variable(identifier+"_alpha",ZplusBsel&&LeadingB&&EtaB&&PhiZcut,20,0,2,"#alpha","Official/"+identifier+"/alpha_full__"+identifier,1);
-+
+   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);
-+
-+
+   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");
-+
+   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");
-+
+   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}");
-+
-+
+}
-+
-+
+void GetNumberEventsInsideOutsideAlphaWindow(TCut specialcut, string title) {
-+
-+
+//  TCut cut =  ZplusBsel&&LeadingB&&PhiZcut&&specialcut&&TCut("ZbCHS3010_alpha<0.3");
-+
+  TCut cut = specialcut;
-+
+  TCut in  =  EtaB;
-+
+  TCut out =  TCut("abs(ZbCHS3010_pfJetGoodEta[0])>1.3");
-+
-+
+  TH1F *data_IN  = allsamples.Draw("data_IN",    "mll",1,0,150, "nothing [GeV]", "events", cut&&in, data,luminosity);
-+
+  TH1F *data_OUT = allsamples.Draw("data_OUT",   "mll",1,0,150, "nothing [GeV]", "events", cut&&out,data,luminosity);
-+
+  TH1F *mc_IN    = allsamples.Draw("mc_IN",      "mll",1,0,150, "nothing [GeV]", "events", cut&&in, mc  ,luminosity);
-+
+  TH1F *mc_OUT   = allsamples.Draw("mc_OUT",     "mll",1,0,150, "nothing [GeV]", "events", cut&&out,mc  ,luminosity);
-+
-+
+  dout << title << " : " << endl;
-+
+  dout << "  Data: In " << data_IN->Integral() << " vs out " << data_OUT->Integral() << endl;
-+
+  dout << "  MC  : In " << mc_IN->Integral()   << " vs out " << mc_OUT->Integral()   << endl;
-+
+  dout << "       ratio in  : " << data_IN->Integral()/mc_IN->Integral() << " +/- " << (data_IN->Integral()/mc_IN->Integral()) * sqrt(1.0/mc_IN->Integral()+ 1.0/data_IN->Integral()) << endl;
-+
+  dout << "       ratio out : " << data_OUT->Integral()/mc_OUT->Integral() << " +/- " << (data_OUT->Integral()/mc_OUT->Integral()) * sqrt(1.0/mc_OUT->Integral()+ 1.0/data_OUT->Integral()) << endl;
-+
-+
+  delete data_IN;
-+
+  delete data_OUT;
-+
+  delete mc_IN;
-+
+  delete mc_OUT;
-+
+}
-+
-+
+void TEST_GetNumberEventsInsideOutsideAlphaWindow() {
-+
+  TCanvas *ca = new TCanvas("ca","ca");
-+
+//   GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==0"), "ee");
-+
+//   GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==1"), "mm");
-+
+//   GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2"), "ee/mm");
-+
-+
+  dout << "BASE SELECTION" << endl;
-+
+  GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==0"), "ee");
-+
+  GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2&&id1==1"), "mm");
-+
+  GetNumberEventsInsideOutsideAlphaWindow(TCut("id1==id2"), "ee/mm");
-+
-+
+  dout << "BASE SELECTION: Z+b" << endl;
-+
+  GetNumberEventsInsideOutsideAlphaWindow(ZplusBsel&&TCut("id1==id2&&id1==0"), "ee");
-+
+  GetNumberEventsInsideOutsideAlphaWindow(ZplusBsel&&TCut("id1==id2&&id1==1"), "mm");
-+
+  GetNumberEventsInsideOutsideAlphaWindow(ZplusBsel&&TCut("id1==id2"), "ee/mm");
-+
-+
+  dout << "Leading B" << endl;
-+
+  GetNumberEventsInsideOutsideAlphaWindow(LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==0"), "ee");
-+
+  GetNumberEventsInsideOutsideAlphaWindow(LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==1"), "mm");
-+
+  GetNumberEventsInsideOutsideAlphaWindow(LeadingB&&ZplusBsel&&TCut("id1==id2"), "ee/mm");
-+
-+
+  dout << "PhiZcut" << endl;
-+
+  GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==0"), "ee");
-+
+  GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==1"), "mm");
-+
+  GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2"), "ee/mm");
-+
-+
+  dout << "alpha cut" << endl;
-+
+  GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==0")&&TCut("ZbCHS3010_alpha<0.3"), "ee");
-+
+  GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2&&id1==1")&&TCut("ZbCHS3010_alpha<0.3"), "mm");
-+
+  GetNumberEventsInsideOutsideAlphaWindow(PhiZcut&&LeadingB&&ZplusBsel&&TCut("id1==id2")&&TCut("ZbCHS3010_alpha<0.3"), "ee/mm");
-+
-+
+  delete ca;
-+
-+
-+
+}
-+
-+
-+
+void MetSpectrumDepletionIllustration(int isdata) {
-+
-+
+  TCanvas *can = new TCanvas("can","can");
-+
+  TCut basecut(ZplusBsel&&LeadingB&&EtaB&&PhiZcut);
-+
-+
+  cutWeight=STDWEIGHT*TCut("ZbCHS3010_BTagWgtTUp*(1.0)*(ZbCHS3010_bTagProbCSVBP[0]>0.898)");
-+
+  TH1F *standard  = allsamples.Draw("standard","met[1]",80,0,400, "MET [GeV]", "events", basecut,isdata,luminosity);
-+
-+
+  cutWeight=STDWEIGHT*TCut("ZbCHS3010_BTagWgtTUp*(1.0*softMuon)*(ZbCHS3010_bTagProbCSVBP[0]>0.898)");
-+
+  TH1F *Neutrinos  = allsamples.Draw("Neutrinos","met[1]",80,0,400, "MET [GeV]", "events", basecut,isdata,luminosity);
-+
-+
+  cutWeight=STDWEIGHT*TCut("ZbCHS3010_BTagWgtTUp*(1.0*(1.0-softMuon))*(ZbCHS3010_bTagProbCSVBP[0]>0.898)");
-+
+  TH1F *NoNeutrinos  = allsamples.Draw("NoNeutrinos","met[1]",80,0,400, "MET [GeV]", "events", basecut,isdata,luminosity);
-+
-+
+  standard->Scale(1.0/standard->Integral());
-+
+  Neutrinos->Scale(1.0/Neutrinos->Integral());
-+
+  NoNeutrinos->Scale(1.0/NoNeutrinos->Integral());
-+
-+
+  Neutrinos->SetLineColor(kGreen);
-+
+  NoNeutrinos->SetLineColor(kRed);
-+
-+
+  can->SetLogy(1);
-+
+  TLegend *leg = make_legend();
-+
+  leg->AddEntry(standard,"Standard sample","ep");
-+
+  leg->AddEntry(Neutrinos,"#nu enriched sample","l");
-+
+  leg->AddEntry(NoNeutrinos,"#nu depleted sample","l");
-+
+  leg->SetX1(0.9*leg->GetX1());
-+
+  standard->Draw("e1");
-+
+  Neutrinos->Draw("histo,same");
-+
+  NoNeutrinos->Draw("histo,same");
-+
+  if(isdata==data) leg->SetHeader("Data:");
-+
+  else leg->SetHeader("MC:");
-+
+  leg->Draw();
-+
+  DrawPrelim();
-+
-+
+  if(isdata==data) CompleteSave(can,"NeutrinoMetSpectrum__Data");
-+
+  else CompleteSave(can,"NeutrinoMetSpectrum__MC");
-+
+  delete standard;
-+
+  delete Neutrinos;
-+
+  delete NoNeutrinos;
-+
+  delete can;
-+
+  CleanLegends();
-+
+}
-+
-+
+string PrintVector(vector<float> Systematics) {
-+
+  stringstream fullanswer;
-+
+  fullanswer << " [ ";
-+
+  if(Systematics.size()>0) fullanswer << Systematics[0];
-+
+  for(unsigned int is=1;is<Systematics.size();is++) fullanswer << " , " << Systematics[is];
-+
+  fullanswer << " ] ";
-+
+  return fullanswer.str();
-+
+}
-+
-+
-+
+void do_basic_ZB_analysis(bool DoCompleteAnalysis) {
-+
+  STDWEIGHT=TCut(cutWeight);
-+
+  cutWeight=TCut(STDWEIGHT*TCut("ZbCHS3010_BTagWgtT"));
-+
+  dout << "The lepton requirement is " << (const char*) leptoncut << endl;
-+
+  bool doquick=false;
-+
-+
+   TCanvas *zbcanvas = new TCanvas("zbcanvas","zbcanvas");
-+
-+
-+
+//   MetSpectrumDepletionIllustration(data);
-+
+//   MetSpectrumDepletionIllustration(mc);
-+
-+
+//   compare_selection("ZbCHS3010_p5Clean");
-+
+//   compare_selection("ZbCHS3010");
-+
+//   compare_selection("Zb30");
-+
+//   compare_selection("Zb3010");
-+
+//   compare_selection("Zb1510");
-+
+//   compare_selection("ZbCHS30");
-+
+//   compare_selection("ZbCHS3010");
-+
+//   compare_selection("ZbCHS3010");
-+
+//   compare_selection("ZbCHS301010");
-+
+//   compare_selection("ZbCHS3010_SecEta3");
-+
+//   compare_selection("ZbCHS3010_SecEta5");
-+
+//   compare_selection("Zb1510");
-+
+//   compare_selection("ZbCHS301010");
-+
+//   compare_selection("ZbMikko");
-+
+//   compare_selection("ZbCHS3010");
-+
+//   compare_selection("Zb40");
-+
-+
+  if(!doquick) {
-+
+//    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.)
-+
+//    print_all_b_yields();
-+
+//     draw_mpf_vars();
-+
+//    DrawEvilCutFlow();
-+
-+
+//     draw_Zb_kin_vars();
-+
-+
+   }
-+
-+
+//  GetNumberEventsInsideOutsideAlphaWindow();
-+
-+
+  if(DoCompleteAnalysis) {
-+
+    /* need to run the following scenarios:
-+
+     * - normal
-+
+     * - inclusive
-+
+     * - medium WP
-+
+     * - alpha up
-+
+     * - alpha down
-+
+     */
-+
-+
+    bool fast=true;bool slow=false;//don't change these
-+
-+
+    ZbResultContainer inclusive = new_data_mc_agreement_2d(slow,"ZbCHS3010_alpha","inclusive",TCut("mll>0"),TCut("1.0"),0.3);
-+
+    ZbResultContainer Reference = new_data_mc_agreement_2d(slow,"ZbCHS3010_alpha","reference",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtT"),0.3);
-+
-+
-+
+    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);
-+
+    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);
-+
-+
+    ZbResultContainer effmistagUp = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","effmistagUp",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtTUp"),0.3);
-+
+    ZbResultContainer effmistagDown = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","effmistagDown",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtTDown"),0.3);
-+
+    ZbResultContainer medium = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","medium",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.679"),TCut("ZbCHS3010_BTagWgtM"),0.3);
-+
+    ZbResultContainer loose = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","loose",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.244"),TCut("ZbCHS3010_BTagWgtL"),0.3);
-+
+    ZbResultContainer AlphaUp = new_data_mc_agreement_2d(fast,"ZbCHS3010_alphaUp","AlphaUp",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtT"),0.3);
-+
+    ZbResultContainer AlphaDown = new_data_mc_agreement_2d(fast,"ZbCHS3010_alphaDown","AlphaDown",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtT"),0.3);
-+
-+
+    ZbResultContainer L5inclusive = new_data_mc_agreement_2d(fast,"ZbCHS3010_alphaL5","L5inclusive",TCut("mll>0"),TCut("1.0"),0.3);
-+
+    ZbResultContainer L5Reference = new_data_mc_agreement_2d(fast,"ZbCHS3010_alphaL5","L5reference",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtT"),0.3);
-+
-+
+    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);
-+
+    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);
-+
+    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);
-+
+    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);
-+
+    ZbResultContainer AlphaThresholdVariationDown10inc = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","AlphaThresholdVariationDown10inc",TCut("mll>0"),                              TCut("1.0"),0.3-0.1*0.3);
-+
+    ZbResultContainer AlphaThresholdVariationUp10inc =   new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","AlphaThresholdVariationUp10inc",  TCut("mll>0"),                              TCut("1.0"),0.3+0.1*0.3);
-+
+    ZbResultContainer AlphaThresholdVariationDown30inc = new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","AlphaThresholdVariationDown30inc",TCut("mll>0"),                              TCut("1.0"),0.3-0.3*0.3);
-+
+    ZbResultContainer AlphaThresholdVariationUp30inc =   new_data_mc_agreement_2d(fast,"ZbCHS3010_alpha","AlphaThresholdVariationUp30inc",  TCut("mll>0"),                              TCut("1.0"),0.3+0.3*0.3);
-+
-+
-+
+    ZbResultContainer OnlyOneJet   = new_data_mc_agreement_2d(slow,"ZbCHS3010_alpha","reference",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898&&ZbCHS3010_pfJetGoodNum==1"),TCut("ZbCHS3010_BTagWgtT"),0.3);
-+
+    ZbResultContainer MultipleJets = new_data_mc_agreement_2d(slow,"ZbCHS3010_alpha","reference",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898&&ZbCHS3010_pfJetGoodNum>1"),TCut("ZbCHS3010_BTagWgtT"),0.3);
-+
-+
+    ZbResultContainer OnlyOneJetinc   = new_data_mc_agreement_2d(slow,"ZbCHS3010_alpha","reference",TCut("ZbCHS3010_pfJetGoodNum==1"),TCut("1.0"),0.3);
-+
+    ZbResultContainer MultipleJetsinc = new_data_mc_agreement_2d(slow,"ZbCHS3010_alpha","reference",TCut("ZbCHS3010_pfJetGoodNum>1"),TCut("1.0"),0.3);
-+
-+
+    //and now let's compute some systematics!
-+
+    //1 alpha variation
-+
+    float SysMPF_Alpha_down = Reference.MPF_Result_FaceValue.getValue()-AlphaDown.MPF_Result_FaceValue.getValue();
-+
+    float SysMPF_Alpha_up = Reference.MPF_Result_FaceValue.getValue()-AlphaUp.MPF_Result_FaceValue.getValue();
-+
-+
+    vector<float> SysMPF_Alpha_down_Bins;
-+
+    vector<float> SysMPF_Alpha_up_Bins;
-+
+    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());
-+
+    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());
-+
-+
-+
+    //2 btagging efficiency/mistag correction
-+
+    float SysMPF_EFF_down = Reference.MPF_Result_FaceValue.getValue()-effmistagDown.MPF_Result_FaceValue.getValue();
-+
+    float SysMPF_EFF_up = Reference.MPF_Result_FaceValue.getValue()-effmistagUp.MPF_Result_FaceValue.getValue();
-+
-+
+    vector<float> SysMPF_EFF_down_Bins;
-+
+    vector<float> SysMPF_EFF_up_Bins;
-+
+    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());
-+
+    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());
-+
-+
+    //3 Presence of soft muons (->neutrino question)
-+
+    float SysMPF_Neutrino_up = Reference.MPF_Result_FaceValue.getValue()-NeutrinoQ.MPF_Result_FaceValue.getValue();
-+
+    float SysMPF_Neutrino_down = Reference.MPF_Result_FaceValue.getValue()-ANeutrinoQ.MPF_Result_FaceValue.getValue();
-+
-+
+    vector<float> SysMPF_Neutrino_up_Bins;
-+
+    vector<float> SysMPF_Neutrino_down_Bins;
-+
+    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());
-+
+    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());
-+
-+
-+
+    //4 purity
-+
+    zbcanvas->cd();
-+
+    gStyle->SetOptFit(0);
-+
+    write_info(__FUNCTION__,"Don't know the exact purity at the moment, still needs to be determined!");
-+
+    TH1F *purityhisto = new TH1F("purityhisto","purityhisto",1,0,1);
-+
+    purityhisto->GetXaxis()->SetTitle("Purity");
-+
+    purityhisto->GetYaxis()->SetTitle("C_{abs}");
-+
+    purityhisto->GetXaxis()->CenterTitle();
-+
+    purityhisto->GetYaxis()->CenterTitle();
-+
-+
+    TGraphErrors *gSysMPF_Purity = new TGraphErrors();
-+
-+
+    gSysMPF_Purity->SetPoint(0,PurityInclusive,inclusive.MPF_Result_FaceValue.getValue());
-+
+    gSysMPF_Purity->SetPointError(0,0.0,inclusive.MPF_Result_FaceValue.getError());
-+
-+
+    gSysMPF_Purity->SetPoint(1,PurityLoose,loose.MPF_Result_FaceValue.getValue());//x value?
-+
+    gSysMPF_Purity->SetPointError(1,0.0,loose.MPF_Result_FaceValue.getError());
-+
+    gSysMPF_Purity->SetPoint(2,PurityMedium,medium.MPF_Result_FaceValue.getValue());//x value?
-+
+    gSysMPF_Purity->SetPointError(2,0.0,medium.MPF_Result_FaceValue.getError());
-+
+    gSysMPF_Purity->SetPoint(3,PurityTight,Reference.MPF_Result_FaceValue.getValue());//x value?
-+
+    gSysMPF_Purity->SetPointError(3,0.0,Reference.MPF_Result_FaceValue.getError());
-+
-+
-+
+    float min,max;
-+
+    FindMinMax(gSysMPF_Purity,min,max);
-+
+    purityhisto->SetStats(0);
-+
+    purityhisto->GetYaxis()->SetRangeUser(min,max);
-+
-+
+    gSysMPF_Purity->Fit("pol1");
-+
+    TPolyLine *purity_fit_uncert = GetFitUncertaintyShape(gSysMPF_Purity, "pol1", min, max,0.0,1.0);
-+
+    TF1* pufit = (TF1*)gSysMPF_Purity->GetFunction("pol1");
-+
-+
+    TF1* pufit_c = new TF1("pufit_c","[0]+[1]*x");
-+
+    pufit_c->SetParameters(pufit->GetParameters());
-+
+    pufit_c->SetLineColor(TColor::GetColor("#0101DF"));
-+
+    purity_fit_uncert->SetFillColor(TColor::GetColor("#5353FC"));
-+
+    pufit->SetLineColor(TColor::GetColor("#5353FC"));
-+
-+
+    purityhisto->Draw();
-+
+    purity_fit_uncert->Draw("f");
-+
+    gSysMPF_Purity->Draw("P*");
-+
+    pufit_c->Draw("same");
-+
+    DrawPrelim();
-+
-+
+    float ExtrapolatedResult=pufit->GetParameter(0)+1.0*pufit->GetParameter(1);
-+
+    float dExtrapolatedResult=sqrt(pufit->GetParError(0)*pufit->GetParError(0)+1.0*1.0*pufit->GetParError(1)*pufit->GetParError(1));
-+
-+
+    stringstream summary;
-+
+    summary << "#splitline{Reference: " << std::setprecision(4) << Reference.MPF_Result_FaceValue.getValue() << "+/-" << std::setprecision(4) << Reference.MPF_Result_FaceValue.getError() << "}{";
-+
+    summary << "Extrapolation: " << std::setprecision(4) << ExtrapolatedResult << " +/- " << std::setprecision(4) << dExtrapolatedResult << "}";
-+
-+
+    dout << "Have found an extrapolated result at 1.0 of " << ExtrapolatedResult << " ;  will use the difference between this result and the face value as ";
-+
+    dout << "purity-related systematic uncertainty" << endl;
-+
-+
+    TText *infobox = write_title(summary.str());
-+
+    infobox->SetX(0.75);
-+
+    infobox->SetTextSize(0.03);
-+
+    infobox->SetY(0.75);
-+
+    infobox->Draw();
-+
-+
+    CompleteSave(zbcanvas,"Systematics/Purity");
-+
-+
+    dout << "Bin-by-bin systematics : " << endl;
-+
+    dout << "    Oversmearing : " << endl;
-+
+    dout << "      down: " << PrintVector(SysMPF_Alpha_down_Bins) << endl;
-+
+    dout << "      up  : " << PrintVector(SysMPF_Alpha_up_Bins) << endl;
-+
+    dout << "    eff/mistag variation: " << endl;
-+
+    dout << "      down: " << PrintVector(SysMPF_EFF_down_Bins) << endl;
-+
+    dout << "      up  : " << PrintVector(SysMPF_EFF_up_Bins) << endl;
-+
+    dout << "    neutrinos: " << endl;
-+
+    dout << "      down: " << PrintVector(SysMPF_Neutrino_down_Bins) << endl;
-+
+    dout << "      up  : " << PrintVector(SysMPF_Neutrino_up_Bins) << endl;
-+
-+
-+
+    float SysMPF_Purity = abs(ExtrapolatedResult-Reference.MPF_Result_FaceValue.getValue());
-+
-+
+    float sys_alpha    = abs(SysMPF_Alpha_down)>abs(SysMPF_Alpha_up)?abs(SysMPF_Alpha_down):abs(SysMPF_Alpha_up);
-+
+    float sys_eff      = abs(SysMPF_EFF_down)>abs(SysMPF_EFF_up)?abs(SysMPF_EFF_down):abs(SysMPF_EFF_up);
-+
+    float sys_purity   = abs(SysMPF_Purity);
-+
+    float sys_neutrino = abs(SysMPF_Neutrino_down)>abs(SysMPF_Neutrino_up)?abs(SysMPF_Neutrino_down):abs(SysMPF_Neutrino_up);
-+
+    float sys          = sqrt(sys_alpha*sys_alpha+sys_eff*sys_eff+sys_purity*sys_purity+sys_neutrino*sys_neutrino);
-+
-+
+    dout << "MPF METHOD:" << endl;
-+
+    dout << "   Systematics :    down    up     (selected)" << endl;
-+
+    dout << "      Oversmearing         : " << SysMPF_Alpha_down << "   ,  " << SysMPF_Alpha_up << "   ( " << sys_alpha << ") " << endl;
-+
+    dout << "      eff/mistag variation : " << SysMPF_EFF_down << "   ,  " << SysMPF_EFF_up << "   ( " << sys_eff << ") " << endl;
-+
+    dout << "      purity               : " << SysMPF_Purity << "   ( " << sys_purity << " )" << endl;
-+
+    dout << "      neutrino             : " << SysMPF_Neutrino_down << "   ,  " << SysMPF_Neutrino_up << "   ( " << sys_neutrino << ") " << endl;
-+
+    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;
-+
+    dout << endl << endl;
-+
+    dout << "   FINAL RESULTS : " << endl;
-+
+    dout << "     C_{abs}^{b}   = " << Reference.MPF_Result_FaceValue << " (stat) +/- " << sys << " (sys) " << endl;
-+
+    dout << "     C_{abs}^{inc} = " << inclusive.MPF_Result_FaceValue << " (stat) " << endl;
-+
-+
+    float sysfinal = Reference.MPF_Result_FaceValue.getValue()/inclusive.MPF_Result_FaceValue.getValue();
-+
+    sysfinal *= sqrt(sys*sys/(Reference.MPF_Result_FaceValue.getValue()*Reference.MPF_Result_FaceValue.getValue()));
-+
+    //yes, this could be simplified but it would look like a bug.
-+
-<
+void do_basic_ZB_analysis() {
-<
+  TCanvas *zbcanvas = new TCanvas("zbcanvas","zbcanvas");
-<
+  write_warning(__FUNCTION__,"TESTING NEW DATA MC AGREEMENT");
-<
-<
+  /*
-<
+  write_warning(__FUNCTION__,"Not doing b yields (already implemented)");
-<
+//  print_all_b_yields();
-<
+  write_warning(__FUNCTION__,"Not drawing mpf variables (already implemented)");
-<
+//  draw_mpf_vars();
-<
+  write_warning(__FUNCTION__,"Not drawing mpf variables (already implemented)");
-<
+  draw_Zb_kin_vars();
-<
-<
+  data_mc_agreement_2d();
-<
+  */
-<
+  new_data_mc_agreement_2d();
->
->
->
+    dout << "Note: Varying the alpha cut (->face value) we get the following results: " << endl;
->
+    dout << "-30%: " << AlphaThresholdVariationDown30.MPF_Result_FaceValue/AlphaThresholdVariationDown30inc.MPF_Result_FaceValue << " (stat)" << endl;
->
+    dout << "-10%: " << AlphaThresholdVariationDown10.MPF_Result_FaceValue/AlphaThresholdVariationDown10inc.MPF_Result_FaceValue << " (stat)" << endl;
->
+    dout << "  0%: " << Reference.MPF_Result_FaceValue/inclusive.MPF_Result_FaceValue << " (stat) +/- " << endl;
->
+    dout << "+10%: " << AlphaThresholdVariationUp10.MPF_Result_FaceValue/AlphaThresholdVariationUp10inc.MPF_Result_FaceValue << " (stat)" << endl;
->
+    dout << "+30%: " << AlphaThresholdVariationUp30.MPF_Result_FaceValue/AlphaThresholdVariationUp30inc.MPF_Result_FaceValue << " (stat)" << endl;
->
->
+    dout << "          XCheck for Rabs, varying the alpha cut (->face value) we get the following results: " << endl;
->
+    dout << "             -30%: " << AlphaThresholdVariationDown30.Rabs_Result_FaceValue/AlphaThresholdVariationDown30inc.Rabs_Result_FaceValue << " (stat)" << endl;
->
+    dout << "             -10%: " << AlphaThresholdVariationDown10.Rabs_Result_FaceValue/AlphaThresholdVariationDown10inc.Rabs_Result_FaceValue << " (stat)" << endl;
->
+    dout << "               0%: " << Reference.Rabs_Result_FaceValue/inclusive.Rabs_Result_FaceValue << " (stat) +/- " << endl;
->
+    dout << "             +10%: " << AlphaThresholdVariationUp10.Rabs_Result_FaceValue/AlphaThresholdVariationUp10inc.Rabs_Result_FaceValue << " (stat)" << endl;
->
+    dout << "             +30%: " << AlphaThresholdVariationUp30.Rabs_Result_FaceValue/AlphaThresholdVariationUp30inc.Rabs_Result_FaceValue << " (stat)" << endl;
->
->
->
+    dout << "  Also checked what happens when you only use events with one jet and only with more than one jets: " << endl;
->
+    dout << "      Only 1      : " << OnlyOneJet.MPF_Result_FaceValue/OnlyOneJetinc.MPF_Result_FaceValue << " (stat) " << endl;
->
+    dout << "      More than 1 : " << MultipleJets.MPF_Result_FaceValue/MultipleJetsinc.MPF_Result_FaceValue << " (stat) " << endl;
->
->
+    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());
->
+    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());
->
+    float sys_alphavar = sys_alphavar_up>sys_alphavar_down?sys_alphavar_up:sys_alphavar_down;
->
->
+    dout << " -> Additional systematic (from 30%) : " << sys_alphavar << endl;
->
+    sysfinal=sqrt(sysfinal*sysfinal+sys_alphavar*sys_alphavar);
->
+    dout << "Final systematic is : " << sysfinal << endl;
->
->
+    dout << "   ******************************************** <  FINAL RESULT > ******************************************** " << endl;
->
+    dout << "     C_{corr}      = " << Reference.MPF_Result_FaceValue/inclusive.MPF_Result_FaceValue << " (stat) +/- " << sysfinal << " (sys) " << endl;
->
+    dout << "   ******************************************** < /FINAL RESULT > ******************************************** " << endl;
->
->
+    dout << "Note that if L5 corrections are applied, we get " << L5Reference.MPF_Result_FaceValue/L5inclusive.MPF_Result_FaceValue << " (stat) " << endl;
->
->
+    LeadingB=TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898");
->
->
+  } else {
->
+    ZbResultContainer Reference = new_data_mc_agreement_2d(false,"ZbCHS3010_alpha","reference",TCut("ZbCHS3010_bTagProbCSVBP[0]>0.898"),TCut("ZbCHS3010_BTagWgtT"),0.3);
->
+  }
->
+  delete zbcanvas;
-+
-+
-+
-+
-+
+//*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*  DELETE EVERYTHING BELOW THIS LINE /*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*//
-+
-+
-+
+ const char* concatenate(string bla1, string bla2) {
-+
+  stringstream bla;
-+
+  bla << bla1 << bla2;
-+
+  return bla.str().c_str();
-+
+}
-+
-+
+void SpecialCutFlow(string identifier) {
-+
+  stringstream MegaCut;
-+
-+
-+
-+
+  MegaCut << "   1*(" << (const char*) (ZplusBsel&&TCut(concatenate(identifier,"_pfJetGoodNumBtag>0"))) << ")";
-+
+  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB) << ")";
-+
+  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB) << ")";
-+
+  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut) << ")";
-+
+  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>30&&pt<1000").c_str())) << ")";
-+
+  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>30&&pt<1000&&"+identifier+"_alpha<0.3").c_str())) << ")";
-+
+  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>30&&pt<1000&&"+identifier+"_alpha<0.2").c_str())) << ")";
-+
+  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>30&&pt<1000&&"+identifier+"_alpha<0.15").c_str())) << ")";
-+
+  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>30&&pt<1000&&"+identifier+"_alpha<0.1").c_str())) << ")";
-+
+  MegaCut << " + 1*(" << (const char*) (ZplusBsel&&LeadingB&&EtaB&&PhiZcut&&TCut(((string)"pt>30&&pt<1000&&"+identifier+"_alpha<0.05").c_str())) << ")";
-+
+  MegaCut << " ";
-+
-+
+  TCanvas *can = new TCanvas("can","can");
-+
+  can->SetLogy(1);
-+
+  TCut basecut(ZplusBsel);
-+
-+
-+
-+
+  TLegend *leg = allsamples.allbglegend();
-+
-+
-+
+  TH1F *data  = allsamples.Draw("data",    MegaCut.str(),11,-0.5,10.5, "", "events", basecut,data,luminosity);
-+
+  THStack mcm = allsamples.DrawStack("mc", MegaCut.str(),11,-0.5,10.5, "", "events", basecut,mc,luminosity);
-+
-+
+  float runningsum=0;
-+
+  for(int i=data->GetNbinsX();i>0;i--) {
-+
+    runningsum+=data->GetBinContent(i);
-+
+    data->SetBinContent(i,runningsum);
-+
+  }
-+
-+
+  float minimum=data->GetMaximum();
-+
-+
+  TH1F* h;
-+
+  TIter nextOF(mcm.GetHists());
-+
-+
+  float mcsum=0;
-+
+  float zb=0;
-+
+  while ( h = (TH1F*)nextOF() ) {
-+
+    float runningsum=0;
-+
+    minimum=data->GetMaximum();//done deliberately at each step so get the minimum of the last (!) contribution
-+
+    for(int i=h->GetNbinsX();i>0;i--) {
-+
+      runningsum+=h->GetBinContent(i);
-+
+      h->SetBinContent(i,runningsum);
-+
+      if(runningsum<minimum&&runningsum>0) minimum=runningsum;
-+
+    }
-+
+    if(Contains(h->GetName(),"Z_b_")) {
-+
+      zb=h->GetBinContent(h->GetNbinsX());
-+
+    }
-+
+    mcsum+=h->GetBinContent(h->GetNbinsX());
-+
+  }
-+
-+
-+
-+
+  data->SetMinimum(0.05*minimum);
-+
+  can->cd(1)->SetBottomMargin(0.25);
-+
+  data->GetXaxis()->SetBinLabel(1,"Pre-selection");
-+
+  data->GetXaxis()->SetBinLabel(2,"Contains b jet");
-+
+  data->GetXaxis()->SetBinLabel(3,"Leading jet is b-jet");
-+
+  data->GetXaxis()->SetBinLabel(4,"|#eta_{b}|<1.3 ");
-+
+  data->GetXaxis()->SetBinLabel(5,"|#delta#phi(b,Z)|>2.7");
-+
+  data->GetXaxis()->SetBinLabel(6,"30<p_{T}^{Z}<1000 GeV");
-+
+  data->GetXaxis()->SetBinLabel(7,"#alpha < 0.3");
-+
+  data->GetXaxis()->SetBinLabel(8,"#alpha < 0.2");
-+
+  data->GetXaxis()->SetBinLabel(9,"#alpha < 0.15");
-+
+  data->GetXaxis()->SetBinLabel(10,"#alpha < 0.1");
-+
+  data->GetXaxis()->SetBinLabel(11,"#alpha < 0.05");
-+
+  data->GetXaxis()->LabelsOption("v");
-+
-+
+  stringstream purityinfo;
-+
+  purityinfo << "Purity: " << std::setprecision(3) << 100*zb/mcsum << " %";
-+
+  stringstream neventsinfo;
-+
+  neventsinfo << "Nevents: " << data->GetBinContent(data->GetNbinsX());
-+
+  TH1F *crap = new TH1F("crap","",1,0,1);
-+
+  crap->SetLineColor(kWhite);
-+
+  leg->AddEntry(crap,purityinfo.str().c_str(),"l");
-+
+  leg->AddEntry(crap,neventsinfo.str().c_str(),"l");
-+
-+
+  data->Draw("e1");
-+
+  mcm.Draw("histo,same");
-+
+  data->Draw("same,e1,axis");
-+
+  data->Draw("same,e1");
-+
+//  data->Draw("same,TEXT");
-+
-+
+  leg->Draw();
-+
-+
+  DrawPrelim();
-+
-+
+  CompleteSave(can,"CutFlow___"+identifier);
-+
-+
+  delete crap;
-+
+  delete data;
-+
+  CleanLegends();
-+
+  DeleteStack(mcm);
-+
-+
-+
+}

Diff Legend

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

Comparing UserCode/cbrown/Development/Plotting/Modules/ZbTools.C (file contents): Revision 1.3 by buchmann, Thu Jul 12 12:53:41 2012 UTC vs. Revision 1.28 by buchmann, Fri Apr 26 06:52:09 2013 UTC

Diff Legend

Comparing UserCode/cbrown/Development/Plotting/Modules/ZbTools.C (file contents):
Revision 1.3 by buchmann, Thu Jul 12 12:53:41 2012 UTC vs.
Revision 1.28 by buchmann, Fri Apr 26 06:52:09 2013 UTC