Plotting/Modules/Plotting_Functions.C

#include <iostream>
#include <vector>
#include <sys/stat.h>

#include <TCut.h>
#include <TROOT.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TColor.h>
#include <TPaveText.h>
#include <TRandom.h>
#include <TH1.h>
#include <TH2.h>
#include <TF1.h>
#include <TSQLResult.h>
#include <TProfile.h>

//#include "TTbar_stuff.C"
using namespace std;

using namespace PlottingSetup;

void find_peaks(float &MCPeak,float &MCPeakError, float &DataPeak, float &DataPeakError, float &MCSigma, float &DataSigma, stringstream &result)
{
  TH1F *rawJZBeemmMC      = allsamples.Draw("rawJZBeemmMC",jzbvariablemc,100,-50,50, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc, luminosity);
  TH1F *rawJZBeemmData    = allsamples.Draw("rawJZBeemmData",jzbvariabledata,100, -50,50, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data, luminosity);
  TH1F *rawJZBemMC      = allsamples.Draw("rawJZBemMC",jzbvariablemc,100,-50,50, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc, luminosity);
  TH1F *rawJZBemData    = allsamples.Draw("rawJZBemData",jzbvariabledata,100, -50,50, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data, luminosity);
  TH1F *rawttbarjzbeemmMC;

  if(method==doKM) {
    //we only need this histo for the KM fitting...
    rawttbarjzbeemmMC = allsamples.Draw("rawttbarjzbeemmMC",jzbvariablemc,100, -50,50, "JZB [GeV]", "events",cutmass&&cutOSSF&&cutnJets,mc,luminosity,allsamples.FindSample("TTJet"));
    MCPeak=find_peak(rawJZBeemmMC, rawttbarjzbeemmMC, -40, 40, mc, MCPeakError,MCSigma,method);
    DataPeak=find_peak(rawJZBeemmData, rawJZBeemmData, -40, 40, data, DataPeakError,DataSigma,method);
  }
  else {
    TH1F *reducedMC = (TH1F*)rawJZBeemmMC->Clone();
    TH1F *reducedData = (TH1F*)rawJZBeemmData->Clone();
    reducedMC->Add(rawJZBemMC,-1);
    reducedData->Add(rawJZBemData,-1);
    //this is Kostas' way of doing it - we subtract em to get rid of some of the ttbar contribution (in reality, of flavor-symmetric contribution)
    MCPeak=find_peak(reducedMC, rawttbarjzbeemmMC, -40, 40, mc, MCPeakError,MCSigma,method);
    DataPeak=find_peak(reducedData, rawJZBeemmData, -40, 40, data, DataPeakError,DataSigma,method);
    
  }
  
    
//  MCPeak=find_peak(rawJZBeemmMC, rawttbarjzbeemmMC, -40, 40, mc, MCPeakError,MCSigma,method);
//  DataPeak=find_peak(rawJZBeemmData, rawJZBeemmData, -40, 40, data, DataPeakError,DataSigma,method);
  cout << "We have found the peak in Data at " << DataPeak << " +/- " << DataPeakError << " with sigma=" << DataSigma << " +/- ?? (not impl.)" << endl;
  result << "We have found the peak in Data at " << DataPeak << " +/- " << DataPeakError << " with sigma=" << DataSigma << " +/- ?? (not impl.)" << endl;
  cout << "We have found the peak in MC at " << MCPeak << " +/- " << MCPeakError << " with sigma=" << MCSigma << " +/- ?? (not impl.)" << endl;
  result << "We have found the peak in MC at " << MCPeak << " +/- " << MCPeakError << " with sigma=" << MCSigma << " +/- ?? (not impl.)" << endl;
}

void make_special_mll_plot(int nbins, float min, float max, bool logscale,string xlabel) {
  
  TCanvas *ckin = new TCanvas("ckin","Kinematic Plots (in the making)",600,600);

  TH1F *datahistoOSSF = allsamples.Draw("datahistoOSSF","mll",nbins,min,max, xlabel, "events",cutOSSF&&cutnJets&&basiccut,data,luminosity);
  THStack mcstackOSSF = allsamples.DrawStack("mcstackOSSF","mll",nbins,min,max, xlabel, "events",cutOSSF&&cutnJets&&basiccut,mc,luminosity);
  TH1F *datahistoOSOF = allsamples.Draw("datahistoOSSF","mll",nbins,min,max, xlabel, "events",cutOSOF&&cutnJets&&basiccut,data,luminosity);
  
  if(logscale) ckin->SetLogy(1);
  datahistoOSSF->SetMarkerSize(DataMarkerSize);
  datahistoOSSF->GetXaxis()->SetTitle(xlabel.c_str());
  datahistoOSSF->GetXaxis()->CenterTitle();
  datahistoOSSF->GetYaxis()->SetTitle("events");
  datahistoOSSF->GetYaxis()->CenterTitle();
  
  datahistoOSOF->SetMarkerSize(DataMarkerSize);
  datahistoOSSF->SetMarkerSize(DataMarkerSize);
  datahistoOSSF->Draw();
  mcstackOSSF.Draw("same");
  datahistoOSSF->Draw("same");
  datahistoOSOF->SetMarkerColor(TColor::GetColor("#FF4000"));
  datahistoOSOF->SetLineColor(TColor::GetColor("#FF4000"));
  datahistoOSOF->SetMarkerStyle(21);
  datahistoOSOF->Draw("same");
  TLegend *kinleg = allsamples.allbglegend("",datahistoOSSF);
  kinleg->AddEntry(datahistoOSOF,"OSOF (data)","p");
  kinleg->Draw();
  CompleteSave(ckin,"kin/mll_ossf_osof_distribution");
  
  delete datahistoOSOF;
  delete datahistoOSSF;
}
  

void make_plot(string variable, int nbins, float min, float max, bool logscale,string xlabel, string filename,bool isPF=true) {
//  TCut basiccut("(pfJetGoodNum>=2&&pfJetGoodID[0])&&(pfJetGoodNum>=2&&pfJetGoodID[1])&&(passed_triggers||!is_data)");
  TCut ibasiccut=basiccut;
  
  if(isPF) ibasiccut=basiccut&&"pfjzb[0]>-998";
  //Step 1: Adapt the variable (if we're dealing with PF we need to adapt the variable!)
  if(isPF) {
    if(variable=="mll") variable="pfmll[0]";
    if(variable=="jetpt[0]") variable="pfJetGoodPt[0]";
    if(variable=="jeteta[0]") variable="pfJetGoodEta[0]";
    if(variable=="pt") variable="pfpt[0]";
    if(variable=="pt1") variable="pfpt1[0]";
    if(variable=="eta1") variable="pfeta1[0]";
    if(variable=="jzb[1]") variable="pfjzb[0]";
    //if(variable=="pfJetGoodNum") variable="pfJetGoodNum"; // pointless.
  }

  //Step 2: Refine the cut
  TCanvas *ckin = new TCanvas("ckin","Kinematic Plots (in the making)",600,600);
  ckin->SetLogy(logscale);
  TCut cut;
  cut=cutmass&&cutOSSF&&cutnJets&&ibasiccut;
  if(filename=="nJets") cut=cutmass&&cutOSSF&&ibasiccut;
  if(filename=="nJets_nocuts_except_mll_ossf") cut=cutmass&&cutOSSF;
  if(filename=="mll") cut=cutOSSF&&cutnJets&&ibasiccut;
  if(filename=="mll_ee") cut=cutOSSF&&cutnJets&&ibasiccut&&"id1==0";
  if(filename=="mll_mm") cut=cutOSSF&&cutnJets&&ibasiccut&&"id1==1";
  if(filename=="mll_inclusive"||filename=="mll_inclusive_highrange") cut=cutOSSF;
  if(filename=="mll_inclusive_ee") cut=cutOSSF&&"id1==0";
  if(filename=="mll_inclusive_mm") cut=cutOSSF&&"id1==1";
  if(filename=="pfJetGoodEta_0") cut=cutOSSF&&cutmass&&ibasiccut;

  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(variable=="pfJetGoodPt[0]") datahisto->SetMaximum(10*datahisto->GetMaximum());
  if(variable=="pt") datahisto->SetMaximum(10*datahisto->GetMaximum());
  if(filename=="mll_inclusive") datahisto->SetMinimum(1);
  datahisto->Draw("e1");
  mcstack.Draw("same");
  datahisto->Draw("same,e1");
  TLegend *kinleg = allsamples.allbglegend();
  kinleg->Draw();
  TText* write_cut = write_cut_on_canvas(decipher_cut(cut,basicqualitycut));
  write_cut->Draw();
  TText* write_variable = write_text(0.99,0.01,variable);
  write_variable->SetTextAlign(31);
  write_variable->SetTextSize(0.02);
//  write_variable->Draw();
  if(isPF) CompleteSave(ckin,"kin/"+filename+"__PF");
  else CompleteSave(ckin,"kin/"+filename);
  datahisto->Delete();
}

void do_kinematic_plots(bool doPF=false)
{
  bool dolog=true;
  bool nolog=false;
  make_plot("mll",50,50,150,dolog,"m_{ll} [GeV]","mll",doPF);
  make_plot("mll",50,50,150,dolog,"m_{ll} [GeV]","mll_ee",doPF);
  make_plot("mll",50,50,150,dolog,"m_{ll} [GeV]","mll_mm",doPF);
  make_plot("mll",100,50,150,dolog,"m_{ll} [GeV]","mll_inclusive",doPF);
  make_plot("mll",100,50,150,dolog,"m_{ll} [GeV]","mll_inclusive_ee",doPF);
  make_plot("mll",100,50,150,dolog,"m_{ll} [GeV]","mll_inclusive_mm",doPF);
  make_plot("mll",300,50,350,dolog,"m_{ll} [GeV]","mll_inclusive_highrange",doPF);
  make_plot("jetpt[0]",40,0,200,dolog,"leading jet p_{T} [GeV]","pfJetGoodPt_0",doPF);
  make_plot("jeteta[0]",40,-5,5,nolog,"leading jet #eta","pfJetGoodEta_0",doPF);
  make_plot("pt",50,0,200,dolog,"Z p_{T} [GeV]","Zpt",doPF);
  make_plot("pt1",50,0,100,nolog,"p_{T} [GeV]","pt",doPF);
  make_plot("eta1",40,-5,5,nolog,"#eta_{l}","eta",doPF);
  make_plot("jzb[1]",100,-150,150,dolog,"JZB [GeV]","jzb_ossf",doPF);
  make_plot("pfJetGoodNum",8,0.5,8.5,dolog,"nJets","nJets",doPF);
  make_plot("pfJetGoodNum",8,0.5,8.5,dolog,"nJets","nJets_nocuts_except_mll_ossf",doPF);
  make_special_mll_plot(30,50,150,dolog,"m_{ll} [GeV]");
}

void do_kinematic_PF_plots()
{
  do_kinematic_plots(true);
}

void signal_bg_comparison()
{
  TCanvas *can = new TCanvas("can","Signal Background Comparison Canvas");
  int sbg_nbins=100;
  float sbg_min=-200;
  float sbg_max=400;
  can->SetLogy(1);
  TH1F *JZBplotZJETs = allsamples.Draw("JZBplotZJETs",jzbvariablemc,sbg_nbins,sbg_min,sbg_max, "JZB [GeV]", "events",cutmass&&cutOSSF&&cutnJets,mc,luminosity,allsamples.FindSample("DYJetsToLL"));
  TH1F *JZBplotLM4 = allsamples.Draw("JZBplotLM4",jzbvariablemc,sbg_nbins,sbg_min,sbg_max, "JZB [GeV]", "events",cutmass&&cutOSSF&&cutnJets,mc,luminosity,allsamples.FindSample("LM4"));
  TH1F *JZBplotTtbar = allsamples.Draw("JZBplotTtbar",jzbvariablemc,sbg_nbins,sbg_min,sbg_max, "JZB [GeV]", "events",cutmass&&cutOSSF&&cutnJets,mc,luminosity,allsamples.FindSample("TTJets"));
  
  JZBplotTtbar->SetLineColor(allsamples.GetColor("TTJet"));
  JZBplotTtbar->SetLineColor(allsamples.GetColor("TTJet"));
  JZBplotTtbar->SetLineWidth(2);
  JZBplotZJETs->SetFillColor(allsamples.GetColor("DY"));
  JZBplotZJETs->SetLineColor(kBlack);
  JZBplotLM4->SetLineWidth(2);
  JZBplotLM4->SetLineStyle(2);
  JZBplotZJETs->SetMinimum(0.05);
  JZBplotLM4->Scale(20.0);
  JZBplotZJETs->Draw("histo");
  JZBplotLM4->Draw("histo,same");
  JZBplotTtbar->Draw("histo,same");
  TLegend *signal_bg_comparison_leg = make_legend("MC:");
  signal_bg_comparison_leg->AddEntry(JZBplotZJETs,"Z+Jets","f");
  signal_bg_comparison_leg->AddEntry(JZBplotTtbar,"t#bar{t}","f");
  signal_bg_comparison_leg->AddEntry(JZBplotLM4,"20xLM4","f");
  signal_bg_comparison_leg->Draw();
  CompleteSave(can,"jzb_bg_vs_signal_distribution");
  JZBplotLM4->Scale(1.0/20);
  JZBplotZJETs->SetMinimum(0.004);
  JZBplotZJETs->Draw("histo");
  JZBplotLM4->Draw("histo,same");
  JZBplotTtbar->Draw("histo,same");
  TLegend *signal_bg_comparison_leg2 = make_legend("MC:");
  signal_bg_comparison_leg2->AddEntry(JZBplotZJETs,"Z+Jets","f");
  signal_bg_comparison_leg2->AddEntry(JZBplotTtbar,"t#bar{t}","f");
  signal_bg_comparison_leg2->AddEntry(JZBplotLM4,"LM4","f");
  signal_bg_comparison_leg2->Draw();
  CompleteSave(can,"jzb_bg_vs_signal_distribution_unscaled_LM4");
}

//TF1 *BpredFunc,*BpredFuncP,*BpredFuncN;

vector<TF1*> do_cb_fit_to_plot(TH1F *histo, float Sigma, float doingfitacrosstheboard=false) {
  
//  BpredFunc = new TF1("BpredFunc",InvCrystalBall,0,histo->GetBinLowEdge(histo->GetNbinsX())+histo->GetBinWidth(histo->GetNbinsX()),5);
  TF1 *BpredFunc = new TF1("BpredFunc",InvCrystalBall,0,1000,5);
  BpredFunc->SetParameter(0,histo->GetBinContent(1));
  if(doingfitacrosstheboard) BpredFunc->SetParameter(0,histo->GetMaximum());
  BpredFunc->SetParameter(1,0.);
  if(method==1) BpredFunc->SetParameter(2,10*Sigma);//KM
  else BpredFunc->SetParameter(2,Sigma);//Gaussian based methods
  if(method==-99) BpredFunc->SetParameter(2,2.0*Sigma);//Kostas
  BpredFunc->SetParameter(3,1.8);
  BpredFunc->SetParameter(4,2.5);
  histo->Fit(BpredFunc,"QN0");
  BpredFunc->SetLineColor(kBlue);
  BpredFunc->SetLineWidth(1);
  
  TF1 *BpredFuncP = new TF1("BpredFuncP",InvCrystalBallP,-1000,histo->GetBinLowEdge(histo->GetNbinsX())+histo->GetBinWidth(histo->GetNbinsX()),5);
  TF1 *BpredFuncN = new TF1("BpredFuncN",InvCrystalBallN,-1000,histo->GetBinLowEdge(histo->GetNbinsX())+histo->GetBinWidth(histo->GetNbinsX()),5);
  
  BpredFuncP->SetParameters(BpredFunc->GetParameters());
  BpredFuncP->SetLineColor(kBlue);
  BpredFuncP->SetLineWidth(1);
  BpredFuncP->SetLineStyle(2);
  
  BpredFuncN->SetParameters(BpredFunc->GetParameters());
  BpredFuncN->SetLineColor(kBlue);
  BpredFuncN->SetLineWidth(1);
  BpredFuncN->SetLineStyle(2);
  
  vector<TF1*> functions;
  functions.push_back(BpredFuncN);
  functions.push_back(BpredFunc);
  functions.push_back(BpredFuncP);
  return functions;
}
void do_prediction_plot(string jzb, TCanvas *globalcanvas, float sigma, float high, bool is_data, bool do_fit=true)
{
  int nbins=100;
  if(is_data) nbins=50;
  float low=0;
  float hi=500;
  
  TH1F *RcorrJZBeemm   = allsamples.Draw("RcorrJZBeemm",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,is_data, luminosity);
  TH1F *LcorrJZBeemm   = allsamples.Draw("LcorrJZBeemm",("-"+jzb).c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,is_data, luminosity);
  TH1F *RcorrJZBem     = allsamples.Draw("RcorrJZBem",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,is_data, luminosity);
  TH1F *LcorrJZBem     = allsamples.Draw("LcorrJZBem",("-"+jzb).c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,is_data, luminosity);
  
  TH1F *RcorrJZBSBem   = allsamples.Draw("RcorrJZBSBem",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets,is_data, luminosity);
  TH1F *LcorrJZBSBem   = allsamples.Draw("LcorrJZBSBem",("-"+jzb).c_str(),nbins,low,hi, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets,is_data, luminosity);

  TH1F *RcorrJZBSBeemm = allsamples.Draw("RcorrJZBSBeemm",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets,is_data, luminosity);
  TH1F *LcorrJZBSBeemm = allsamples.Draw("LcorrJZBSBeemm",("-"+jzb).c_str(),nbins,low,hi, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets,is_data, luminosity);
  
  TH1F *lm4RcorrJZBeemm   = allsamples.Draw("lm4RcorrJZBeemm",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,is_data, luminosity,allsamples.FindSample("LM4"));
  
  TH1F *Bpred = (TH1F*)LcorrJZBeemm->Clone("Bpred");
  Bpred->Add(RcorrJZBem,1.0/3);
  Bpred->Add(LcorrJZBem,-1.0/3);
  Bpred->Add(RcorrJZBSBem,1.0/3);
  Bpred->Add(LcorrJZBSBem,-1.0/3);
  Bpred->Add(RcorrJZBSBeemm,1.0/3);
  Bpred->Add(LcorrJZBSBeemm,-1.0/3);
  globalcanvas->cd();
  globalcanvas->SetLogy(1);
  Bpred->SetLineColor(kRed);
  Bpred->SetStats(0);
  Bpred->GetXaxis()->SetRangeUser(0,high);
  Bpred->Draw("hist");
  RcorrJZBeemm->SetMarkerStyle(20);
  RcorrJZBeemm->Draw("e1x0,same");
  lm4RcorrJZBeemm->SetLineColor(TColor::GetColor("#088A08"));
  lm4RcorrJZBeemm->Draw("histo,same");
  RcorrJZBeemm->SetMarkerSize(DataMarkerSize);
  
  TLegend *legBpred = make_legend("",0.6,0.5);
  if(is_data)
  {
    vector<TF1*> functions = do_cb_fit_to_plot(Bpred,sigma);
    functions[0]->Draw("same");functions[1]->Draw("same");functions[2]->Draw("same");
    legBpred->AddEntry(RcorrJZBeemm,"observed (data)","p");
    legBpred->AddEntry(Bpred,"predicted (data)","l");
    legBpred->AddEntry(functions[0],"predicted fit (data)","l");
    legBpred->AddEntry(functions[1],"stat. uncert.","l");
    legBpred->AddEntry(lm4RcorrJZBeemm,"LM4","l");
    legBpred->Draw();
    CompleteSave(globalcanvas,"Bpred_Data");
  }
  else {
    legBpred->AddEntry(RcorrJZBeemm,"MC true B","l");
    legBpred->AddEntry(Bpred,"data-driven B","l");
    legBpred->AddEntry(lm4RcorrJZBeemm,"LM4","l");
    legBpred->SetY1(0.65);
    legBpred->Draw();
    
    CompleteSave(globalcanvas,"Bpred_MC");
  }
  
  cout << "Splitting up the different estimates ... " << endl;
  
  TH1F *Bpredem = (TH1F*)LcorrJZBeemm->Clone();
  Bpredem->Add(RcorrJZBem);
  Bpredem->Add(LcorrJZBem,-1);
  TH1F *BpredSBem = (TH1F*)LcorrJZBeemm->Clone();
  BpredSBem->Add(RcorrJZBSBem);
  Bpred->Add(LcorrJZBSBem,-1);
  TH1F *BpredSBeemm = (TH1F*)LcorrJZBeemm->Clone();  
  BpredSBeemm->Add(RcorrJZBSBeemm);
  BpredSBeemm->Add(LcorrJZBSBeemm,-1.0);
  globalcanvas->cd();
  globalcanvas->SetLogy(1);
  Bpredem->SetLineColor(kRed);
  Bpredem->SetStats(0);
  Bpredem->GetXaxis()->SetRangeUser(0,high);
  Bpredem->Draw("hist");
  BpredSBem->SetLineColor(TColor::GetColor("#0B6138"));
  BpredSBem->Draw("hist,same");
  BpredSBeemm->SetLineColor(kBlue);
  BpredSBeemm->Draw("hist,same");
  RcorrJZBeemm->SetMarkerStyle(20);
  RcorrJZBeemm->Draw("e1x0,same");
  RcorrJZBeemm->SetMarkerSize(DataMarkerSize);
  
  TLegend *legBpredc = make_legend("",0.6,0.5);
  if(is_data)
  {
    legBpredc->AddEntry(RcorrJZBeemm,"observed (data)","p");
    legBpredc->AddEntry(Bpredem,"predicted (e#mu method)","l");
    legBpredc->AddEntry(BpredSBem,"predicted (OF sidebands)","l");
    legBpredc->AddEntry(BpredSBeemm,"predicted (SF sidebands)","l");
    legBpredc->Draw();
    CompleteSave(globalcanvas,"Bpred_Data_comparison");
  }
  else {
    legBpredc->AddEntry(RcorrJZBeemm,"observed (MC)","p");
    legBpredc->AddEntry(Bpredem,"predicted (e#mu method)","l");
    legBpredc->AddEntry(BpredSBem,"predicted (OF sidebands)","l");
    legBpredc->AddEntry(BpredSBeemm,"predicted (SF sidebands)","l");
    legBpredc->Draw();
    legBpredc->Draw();
    CompleteSave(globalcanvas,"Bpred_MC_comparison");
  }
  delete RcorrJZBeemm;
  delete LcorrJZBeemm;
  delete RcorrJZBem;
  delete LcorrJZBem;
  delete RcorrJZBSBem;
  delete LcorrJZBSBem;
  delete RcorrJZBSBeemm;
  delete LcorrJZBSBeemm;
  delete lm4RcorrJZBeemm;
}

void do_prediction_plots(string mcjzb, string datajzb, float DataSigma, float MCSigma) {
  TCanvas *globalcanvas = new TCanvas("globalcanvas","Prediction Canvas");
  do_prediction_plot(datajzb,globalcanvas,DataSigma,350,data);
  do_prediction_plot(mcjzb,globalcanvas,MCSigma,300,mc);
}

void do_ratio_plot(int is_data,vector<float> binning, string jzb, TCanvas *can, float high=-9999) {
  bool do_data=0;
  bool dosignal=0;
  if(is_data==1) do_data=1;
  if(is_data==2) dosignal=1;
  TH1F *RcorrJZBeemm   = allsamples.Draw("RcorrJZBeemm",jzb.c_str(),binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,do_data, luminosity,dosignal);
  TH1F *LcorrJZBeemm   = allsamples.Draw("LcorrJZBeemm",("-"+jzb).c_str(),binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,do_data, luminosity,dosignal);
  TH1F *RcorrJZBem     = allsamples.Draw("RcorrJZBem",jzb.c_str(),binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,do_data, luminosity,dosignal);
  TH1F *LcorrJZBem     = allsamples.Draw("LcorrJZBem",("-"+jzb).c_str(),binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,do_data, luminosity,dosignal);
  
  TH1F *RcorrJZBSBem   = allsamples.Draw("RcorrJZBSbem",jzb.c_str(),binning, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets,do_data, luminosity,dosignal);
  TH1F *LcorrJZBSBem   = allsamples.Draw("LcorrJZBSbem",("-"+jzb).c_str(),binning, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets,do_data, luminosity,dosignal);
  TH1F *RcorrJZBSBeemm = allsamples.Draw("RcorrJZBSbeemm",jzb.c_str(),binning, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets,do_data, luminosity,dosignal);
  TH1F *LcorrJZBSbeemm = allsamples.Draw("LcorrJZBSbeemm",("-"+jzb).c_str(),binning, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets,do_data, luminosity,dosignal);
  
  TH1F *Bpred = (TH1F*)LcorrJZBeemm->Clone("Bpred");
  Bpred->Add(RcorrJZBem,1.0/3);
  Bpred->Add(LcorrJZBem,-1.0/3);
  Bpred->Add(RcorrJZBSBem,1.0/3);
  Bpred->Add(LcorrJZBSBem,-1.0/3);
  Bpred->Add(RcorrJZBSBeemm,1.0/3);
  Bpred->Add(LcorrJZBSbeemm,-1.0/3);
  can->cd();
  can->SetLogy(0);
  Bpred->SetLineColor(kRed);
  Bpred->SetStats(0);
  if(high>0) Bpred->GetXaxis()->SetRangeUser(0,high);
  TH1F *JZBratioforfitting=(TH1F*)RcorrJZBeemm->Clone("JZBratioforfitting");
  JZBratioforfitting->Divide(Bpred);
  TGraphAsymmErrors *JZBratio = histRatio(RcorrJZBeemm,Bpred,is_data,binning);
  
  
  JZBratio->SetTitle("");
  JZBratio->GetYaxis()->SetRangeUser(0.0,9.0);
//  if(is_data==1) JZBratio->GetXaxis()->SetRangeUser(0,350);
  
  TF1 *pol0 = new TF1("pol0","[0]",0,1000);
  TF1 *pol0d = new TF1("pol0","[0]",0,1000);
//  straightline_fit->SetParameter(0,1);
  JZBratioforfitting->Fit(pol0,"Q0R","",0,30);
  pol0d->SetParameter(0,pol0->GetParameter(0));
  
  JZBratio->GetYaxis()->SetTitle("Observed/Predicted");
  JZBratio->GetXaxis()->SetTitle("JZB [GeV]");
  JZBratio->GetYaxis()->SetNdivisions(519);
  JZBratio->GetYaxis()->SetRangeUser(0.0,9.0);
  JZBratio->GetYaxis()->CenterTitle();
  JZBratio->GetXaxis()->CenterTitle();
  JZBratio->SetMarkerSize(DataMarkerSize);
  JZBratio->Draw("AP");
  /////----------------------------
  TPaveText *writeresult = new TPaveText(0.2,0.78,0.49,0.91,"blNDC");
  writeresult->SetFillStyle(4000);
  writeresult->SetFillColor(kWhite);
  writeresult->SetTextFont(42);
  ostringstream jzb_agreement_data_text;
  jzb_agreement_data_text<< setprecision(2) << "mean =" << pol0->GetParameter(0) << " #pm " << setprecision(1) <<  pol0->GetParError(0);
  if(is_data==1) fitresultconstdata=pol0->GetParameter(0);// data
  if(is_data==0) fitresultconstmc=pol0->GetParameter(0); // monte carlo, no signal
/*  if(is_data) writeresult->AddText("Data closure test");
  else writeresult->AddText("MC closure test");
  */
  writeresult->AddText(jzb_agreement_data_text.str().c_str());
  writeresult->Draw("same");
  pol0d->Draw("same");
  TF1 *topline =  new TF1("","1.2",0,1000);
  TF1 *bottomline =  new TF1("","0.6",0,1000);
  topline->SetLineColor(kBlue);
  topline->SetLineStyle(2);
  bottomline->SetLineColor(kBlue);
  bottomline->SetLineStyle(2);
  topline->Draw("same");
  bottomline->Draw("same");
  TF1 *oneline = new TF1("","1.0",0,1000);
  oneline->SetLineColor(kBlue);
  oneline->SetLineStyle(1); 
  oneline->Draw("same"); 
  TLegend *phony_leg = make_legend("ratio");//this line is just to have the default CMS Preliminary (...) on the canvas as well.
  TLegend *leg = new TLegend(0.60,0.5,0.94,0.89);
  leg->SetTextFont(42);
  if(is_data==1) leg->SetHeader("Ratio (data)");
  else leg->SetHeader("Ratio (MC)");
  
  leg->SetFillStyle(4000);
  leg->SetFillColor(kWhite);
  leg->SetTextFont(42);
//  leg->AddEntry(topline,"+20\% sys envelope","l");
  leg->AddEntry(JZBratio,"obs/pred","p");
  leg->AddEntry(oneline,"ratio = 1","l");
  leg->AddEntry(pol0d,"fit in [0,30] GeV","l");
  leg->AddEntry(bottomline,"sys. envelope","l");
  leg->Draw("same");
  if(is_data==1) CompleteSave(can, "jzb_ratio_data");
  if(is_data==0) CompleteSave(can, "jzb_ratio_mc");
  if(is_data==2) CompleteSave(can, "jzb_ratio_mc_BandS");//special case, MC with signal!
  
  delete RcorrJZBeemm;
  delete LcorrJZBeemm;
  delete RcorrJZBem;
  delete LcorrJZBem;
  
  delete RcorrJZBSBem;
  delete LcorrJZBSBem;
  delete RcorrJZBSBeemm;
  delete LcorrJZBSbeemm;
}

void do_ratio_plots(string mcjzb,string datajzb) {
  TCanvas *globalc = new TCanvas("globalc","Ratio Plot Canvas");
  globalc->SetLogy(0);
  vector<float> ratio_binning; 
  ratio_binning.push_back(0);
  ratio_binning.push_back(5);
  ratio_binning.push_back(10);
  ratio_binning.push_back(20);
  ratio_binning.push_back(50);
  ratio_binning.push_back(100);
  //ratio_binning.push_back(350);
  ratio_binning.push_back(200);
  ratio_binning.push_back(500);
  
  do_ratio_plot(mc,ratio_binning,mcjzb,globalc);
  do_ratio_plot(data,ratio_binning,datajzb,globalc);
  do_ratio_plot(mcwithsignal,ratio_binning,mcjzb,globalc);
}

string give_jzb_expression(float peak, int type) {
  stringstream val;
  if(type==data) {
    if(peak<0) val << jzbvariabledata << "+" << TMath::Abs(peak);
    if(peak>0) val << jzbvariabledata << "-" << TMath::Abs(peak);
    if(peak==0) val << jzbvariabledata;
  }
  if(type==mc) {
    if(peak<0) val << jzbvariablemc << "+" << TMath::Abs(peak);
    if(peak>0) val << jzbvariablemc << "-" << TMath::Abs(peak);
    if(peak==0) val << jzbvariablemc;
  }
  return val.str();
}

string centralupdownname(int num) {
  if (num==0) return "central";
  if (num==1) return "down";
  if (num==2) return "up";
  return "centralupdownnameERROR";
}

string dataormc(int isdata) {
  if(isdata) return "Data";
  else return "MC";
}

int histocounter=0;
string give_histo_number(string basename="", int isdata=-1, int centralcounter=-1) {
  histocounter++;
  stringstream histocounternum;
  if(isdata==-1) {
  if(basename=="") histocounternum << histocounter;
  else histocounternum << basename << "_" << histocounter;
  }
  else {
    histocounternum << basename << "_" << dataormc(isdata) << "_" << centralupdownname(centralcounter);
  }
  return histocounternum.str();
}

int central=0;
int down=1;
int up=2;

void lepton_comparison_plots() {
  TCanvas *can = new TCanvas("can","Lepton Comparison Canvas");
  can->SetLogy(1);
  TH1F *eemc = allsamples.Draw("eemc","mll",50,50,150, "mll [GeV]", "events", cutOSSF&&cutnJets&&"(id1==0)",mc, luminosity,allsamples.FindSample("DYJetsToLL"));
  TH1F *mmmc = allsamples.Draw("mmmc","mll",50,50,150, "mll [GeV]", "events", cutOSSF&&cutnJets&&"(id1==1)",mc, luminosity,allsamples.FindSample("DYJetsToLL"));
  eemc->SetLineColor(kBlue);
  mmmc->SetLineColor(kRed);
  eemc->SetMinimum(0.1);
  eemc->SetMaximum(10*eemc->GetMaximum());
  eemc->Draw("histo");
  mmmc->Draw("histo,same");
  TLegend *leg = make_legend();
  leg->AddEntry(eemc,"ZJets->ee (MC)","l");
  leg->AddEntry(mmmc,"ZJets->#mu#mu (MC)","l");
  leg->Draw("same");
  CompleteSave(can, "mll_effratio_mc");
  
  TH1F *eed = allsamples.Draw("eed","mll",50,50,150, "mll [GeV]", "events", cutOSSF&&cutnJets&&"(id1==0)",data, luminosity);
  TH1F *mmd = allsamples.Draw("mmd","mll",50,50,150, "mll [GeV]", "events", cutOSSF&&cutnJets&&"(id1==1)",data, luminosity);
  eed->SetLineColor(kBlue);
  mmd->SetLineColor(kRed);
  eed->SetMinimum(0.1);
  eed->SetMaximum(10*eed->GetMaximum());
  eed->Draw("histo");
  mmd->Draw("histo,same");
  TLegend *leg2 = make_legend();
  leg2->AddEntry(eed,"ZJets->ee (data)","l");
  leg2->AddEntry(mmd,"ZJets->#mu#mu (data)","l");
  leg2->Draw();
  CompleteSave(can, "mll_effratio_data");
  
  TH1F *jeed = allsamples.Draw("jeed",jzbvariabledata,    90,-100,350, "JZB [GeV]", "events", cutOSSF&&cutnJets&&"(id1==0)",data, luminosity);
  TH1F *jmmd = allsamples.Draw("jmmd",jzbvariabledata,    90,-100,350, "JZB [GeV]", "events", cutOSSF&&cutnJets&&"(id1==1)",data, luminosity);
  TH1F *jeemmd = allsamples.Draw("jeemmd",jzbvariabledata,90,-100,350, "JZB [GeV]", "events", cutOSSF&&cutnJets,data, luminosity);
  cout << "ee : " << jeed->GetMean() << "+/-" << jeed->GetMeanError() << endl;
  cout << "ee : " << jmmd->GetMean() << "+/-" << jmmd->GetMeanError() << endl;
  cout << "eemd : " << jeemmd->GetMean() << "+/-" << jeemmd->GetMeanError() << endl;
  jeemmd->SetLineColor(kBlack);
  jeemmd->SetMarkerStyle(25);
  jeed->SetLineColor(kBlue);
  jmmd->SetLineColor(kRed);
  jeed->SetMinimum(0.1);
  jeed->SetMaximum(10*eed->GetMaximum());
  jeemmd->DrawNormalized();
  jeed->DrawNormalized("histo,same");
  jmmd->DrawNormalized("histo,same");
  jeemmd->DrawNormalized("same");
  TLegend *jleg2 = make_legend("Data");
  jleg2->AddEntry(jeemmd,"ee and #mu#mu","p");
  jleg2->AddEntry(jeed,"ee","l");
  jleg2->AddEntry(jmmd,"#mu#mu","l");
  jleg2->Draw();
  CompleteSave(can,"jzb_effratio_data");
  
  TH1F *zjeed = allsamples.Draw("zjeed",jzbvariabledata,    90,-100,350, "JZB [GeV]", "events", cutOSSF&&cutnJets&&"(id1==0)",mc,  luminosity,allsamples.FindSample("DYJetsToLL_TuneZ2"));
  TH1F *zjmmd = allsamples.Draw("zjmmd",jzbvariabledata,    90,-100,350, "JZB [GeV]", "events", cutOSSF&&cutnJets&&"(id1==1)",mc,  luminosity,allsamples.FindSample("DYJetsToLL_TuneZ2"));
  TH1F *zjeemmd = allsamples.Draw("zjeemmd",jzbvariabledata,90,-100,350, "JZB [GeV]", "events", cutOSSF&&cutnJets,mc,  luminosity,allsamples.FindSample("DYJetsToLL_TuneZ2"));
  cout << "Z+Jets ee : " << zjeed->GetMean() << "+/-" << zjeed->GetMeanError() << endl;
  cout << "Z+Jets ee : " << zjmmd->GetMean() << "+/-" << zjmmd->GetMeanError() << endl;
  cout << "Z+Jets eemd : " << zjeemmd->GetMean() << "+/-" << zjeemmd->GetMeanError() << endl;
  zjeemmd->SetLineColor(kBlack);
  zjeemmd->SetMarkerStyle(25);
  zjeed->SetLineColor(kBlue);
  zjmmd->SetLineColor(kRed);
  zjeed->SetMinimum(0.1);
  zjeed->SetMaximum(10*eed->GetMaximum());
  zjeemmd->DrawNormalized();
  zjeed->DrawNormalized("histo,same");
  zjmmd->DrawNormalized("histo,same");
  zjeemmd->DrawNormalized("same");
  TLegend *zjleg2 = make_legend("Data");
  zjleg2->AddEntry(jeemmd,"ee and #mu#mu","p");
  zjleg2->AddEntry(jeed,"ee","l");
  zjleg2->AddEntry(jmmd,"#mu#mu","l");
  zjleg2->Draw();
  CompleteSave(can,"jzb_effratio_ZJets");
  
  TH1F *ld = allsamples.Draw("ld","mll",50,50,150, "mll [GeV]", "events", cutOSSF&&cutnJets,data, luminosity);
  ld->DrawNormalized("e1");
  eed->DrawNormalized("histo,same");
  mmd->DrawNormalized("histo,same");
  TLegend *leg3 = make_legend();
  leg3->AddEntry(ld,"ZJets->ll (data)","p");
  leg3->AddEntry(eed,"ZJets->ee (data)","l");
  leg3->AddEntry(mmd,"ZJets->#mu#mu (data)","l");
  leg3->Draw();
  CompleteSave(can,"mll_effratio_data__all_compared");
  /*
  TH1F *jzbld  = allsamples.Draw("jzbld",jzbvariable,75,-150,150, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data, luminosity);
  TH1F *jzbemd = allsamples.Draw("jzbemd",jzbvariable,75,-150,150, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data, luminosity);
  */
  TH1F *jzbld  = allsamples.Draw("jzbld",jzbvariabledata,92,-110,350, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data, luminosity);
  TH1F *jzbemd = allsamples.Draw("jzbemd",jzbvariabledata,92,-110,350, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data, luminosity);
  jzbld->SetMarkerColor(kBlack);
  jzbld->SetMarkerStyle(26);
  jzbemd->SetMarkerStyle(25);
  jzbemd->SetMarkerColor(kRed);
  jzbemd->SetLineColor(kRed);
  jzbld->SetMinimum(0.35);
  jzbld->Draw("e1");
  jzbemd->Draw("e1,same");
  TLegend *leg4 = make_legend("Data");
  leg4->AddEntry(jzbld,"ee or #mu#mu","p");
  leg4->AddEntry(jzbemd,"e#mu","p");
  leg4->Draw();
  CompleteSave(can,"jzb_eemumu_emu_data");
  
  TH1F *ttbarjzbld  = allsamples.Draw("ttbarjzbld",jzbvariablemc,110,-150,400, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data, luminosity,allsamples.FindSample("TTJet"));
  TH1F *ttbarjzbemd = allsamples.Draw("ttbarjzbemd",jzbvariablemc,110,-150,400, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data, luminosity,allsamples.FindSample("TTJet"));
  ttbarjzbld->SetLineColor(allsamples.GetColor("TTJet"));
  ttbarjzbemd->SetLineColor(allsamples.GetColor("TTJet"));
  ttbarjzbld->SetLineWidth(2);
  ttbarjzbemd->SetLineWidth(2);
  ttbarjzbld->Draw("histo");
  ttbarjzbemd->SetLineStyle(2);
  ttbarjzbemd->Draw("histo,same");
  TLegend *leg5 = make_legend();
  leg5->AddEntry(ttbarjzbld,"t#bar{t}->(ee or #mu#mu)","l");
  leg5->AddEntry(ttbarjzbemd,"t#bar{t}->e#mu","l");
  leg5->Draw();
  CompleteSave(can,"ttbar_emu_mc");
  
  
}

bool is_OF(TCut cut) {
  string scut = (const char*) cut;
  if((int)scut.find("id1!=id2")>-1) return true;
  if((int)scut.find("id1==id2")>-1) return false;
  return false;
}

void draw_pure_jzb_histo(TCut cut,string variable,string savename, TCanvas *can,float min,float max,int nbins) {
  can->cd();
  can->SetLogy(1);
  string xlabel="JZB [GeV]";

  TH1F *datahisto = allsamples.Draw("datahisto",variable,nbins,min,max, xlabel, "events",cut,data,luminosity);
  THStack mcstack = allsamples.DrawStack("mcstack",variable,nbins,min,max, xlabel, "events",cut,mc,luminosity);
  
  datahisto->SetMinimum(0.1);
  //if(savename=="jzb_OSOF") datahisto->SetMaximum(10);
  datahisto->SetMarkerSize(DataMarkerSize);
  datahisto->Draw("e1");
  mcstack.Draw("same");
  datahisto->Draw("same,e1");
  
  TLegend *leg;
  if(is_OF(cut)) leg = allsamples.allbglegend("Opposite Flavor",datahisto);
  else leg = allsamples.allbglegend("Same Flavor",datahisto);
  leg->Draw();
  string write_cut = decipher_cut(cut,"");
  TText *writeline1 = write_cut_on_canvas(write_cut.c_str());
  writeline1->SetTextSize(0.035);
  writeline1->Draw();
  CompleteSave(can, ("jzb/"+savename));
  
  datahisto->Delete();
  mcstack.Delete();
}

Double_t GausR(Double_t *x, Double_t *par) {
  return gRandom->Gaus(x[0],par[0]);
}
  
void jzb_plots(string mcjzb, string datajzb) {
  TCanvas *can = new TCanvas("can","JZB Plots Canvas");
  float max=350;
  float min=-120;
  int nbins=(max-min)/5.0; // we want 5 GeV/bin
  
  stringstream ss;
//  ss << "GausRandom(" << datajzb << ",0)";
  draw_pure_jzb_histo(cutOSSF&&cutnJets&&cutmass,datajzb,"jzb_OSSF",can,min,max,nbins);
  draw_pure_jzb_histo(cutOSOF&&cutnJets&&cutmass,datajzb,"jzb_OSOF",can,min,max,nbins);
  draw_pure_jzb_histo(cutOSSF&&cutnJets&&sidebandcut,datajzb,"jzb_OSSF_SB",can,min,max,nbins);
  draw_pure_jzb_histo(cutOSOF&&cutnJets&&sidebandcut,datajzb,"jzb_OSOF_SB",can,min,max,nbins);
  draw_pure_jzb_histo(cutOSSF&&cutnJets&&cutmass,datajzb,"jzb_OSSF_rebin",can,min,max,nbins/4.0);
  draw_pure_jzb_histo(cutOSOF&&cutnJets&&cutmass,datajzb,"jzb_OSOF_rebin",can,min,max,nbins/4.0);
  draw_pure_jzb_histo(cutOSSF&&cutnJets&&sidebandcut,datajzb,"jzb_OSSF_SB_rebin",can,min,max,nbins/4.0);
  draw_pure_jzb_histo(cutOSOF&&cutnJets&&sidebandcut,datajzb,"jzb_OSOF_SB_rebin",can,min,max,nbins/4.0);
}

void calculate_all_yields(string mcdrawcommand,vector<float> jzb_cuts) {
  cout << "Calculating background yields in MC:" << endl;
  jzb_cuts.push_back(14000);
  TH1F *allbgs = allsamples.Draw("allbgs",jzbvariablemc,jzb_cuts, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc,luminosity);
  float cumulative=0;
  for(int i=allbgs->GetNbinsX();i>=1;i--) {
    cumulative+=allbgs->GetBinContent(i);
    cout << "Above " << allbgs->GetBinLowEdge(i) << " GeV/c : " << cumulative << endl;
  }
}

void rediscover_the_top(string mcjzb, string datajzb) {
  cout << "Hi! today we are going to (try to) rediscover the top!" << endl;
  TCanvas *c3 = new TCanvas("c3","c3");
  c3->SetLogy(1);
  vector<float> binning;
  //binning=allsamples.get_optimal_binsize(mcjzb,cutmass&&cutOSSF&&cutnJets,20,50,800);
  /*
  binning.push_back(50);
  binning.push_back(100);
  binning.push_back(150);
  binning.push_back(200);
  binning.push_back(500);
  
  
  TH1F *dataprediction = allsamples.Draw("dataprediction",    "-"+datajzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,  luminosity);
  TH1F *puresignal     = allsamples.Draw("puresignal",    datajzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data,  luminosity);
//  TH1F *puresignal     = allsamples.Draw("puresignal",        mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc,  luminosity,allsamples.FindSample("TTJets"));
  TH1F *observed       = allsamples.Draw("observed",          datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
  /*
  ofstream myfile;
  TH1F *ratio = (TH1F*)observed->Clone();
  ratio->Divide(dataprediction);
  ratio->GetYaxis()->SetTitle("Ratio obs/pred");
  ratio->Draw();
  c3->SaveAs("testratio.png");
  myfile.open ("ShapeFit_log.txt");
  establish_upper_limits(observed,dataprediction,puresignal,"LM4",myfile);
  myfile.close();
  */
  
  
  int nbins=100;
  float low=0;
  float hi=500;
  TCanvas *c4 = new TCanvas("c4","c4",900,900);
  c4->Divide(2,2);
  c4->cd(1);
  c4->cd(1)->SetLogy(1);
  TH1F *datapredictiont = allsamples.Draw("datapredictiont",    "-"+datajzb, nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,  luminosity);
  TH1F *datapredictiono = allsamples.Draw("datapredictiono",    "-"+datajzb, nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data,  luminosity);
  datapredictiont->Add(datapredictiono,-1);
  cout << "Second way of doing this !!!! Analytical shape to the left :-D" << endl;
  vector<TF1*> functions = do_cb_fit_to_plot(datapredictiont,10);
  datapredictiont->SetMarkerColor(kRed);
  datapredictiont->SetLineColor(kRed);
  datapredictiont->Draw();
  functions[1]->Draw("same");
  TText *title1 = write_title("Top Background Prediction (JZB<0, with osof subtr)");
  title1->Draw();
  
  c4->cd(2);
  c4->cd(2)->SetLogy(1);
  TH1F *observedt = allsamples.Draw("observedt", datajzb,  nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,  luminosity);
  observedt->Draw();
  datapredictiont->Draw("histo,same");
  functions[1]->Draw("same");
  TText *title2 = write_title("Observed and predicted background");
  title2->Draw();
  
  c4->cd(3);
  c4->cd(3)->SetLogy(1);
//  TH1F *ratio = (TH1F*)observedt->Clone();

  TH1F *analytical_background_prediction= new TH1F("analytical_background_prediction","",nbins,low,hi);
  for(int i=0;i<=nbins;i++) {
    analytical_background_prediction->SetBinContent(i+1,functions[1]->Eval(((hi-low)/((float)nbins))*(i+0.5)));
    analytical_background_prediction->SetBinError(i+1,TMath::Sqrt(functions[1]->Eval(((hi-low)/((float)nbins))*(i+0.5))));
  }
  analytical_background_prediction->GetYaxis()->SetTitle("JZB [GeV]");
  analytical_background_prediction->GetYaxis()->CenterTitle();
  TH1F *analyticaldrawonly = (TH1F*)analytical_background_prediction->Clone();
  analytical_background_prediction->SetFillColor(TColor::GetColor("#3399FF"));
  analytical_background_prediction->SetMarkerSize(0);
  analytical_background_prediction->Draw("e5");
  analyticaldrawonly->Draw("histo,same");
  functions[1]->Draw("same");
  TText *title3 = write_title("Analytical bg pred histo");
  title3->Draw();
  
  c4->cd(4);
//  c4->cd(4)->SetLogy(1);
  vector<float> ratio_binning;
  ratio_binning.push_back(0);
  ratio_binning.push_back(5);
  ratio_binning.push_back(10);
  ratio_binning.push_back(20);
  ratio_binning.push_back(50);
//  ratio_binning.push_back(60);
/*
  ratio_binning.push_back(51);
  ratio_binning.push_back(52);
  ratio_binning.push_back(53);
  ratio_binning.push_back(54);
  ratio_binning.push_back(55);
  ratio_binning.push_back(56);
  ratio_binning.push_back(57);
  ratio_binning.push_back(58);
  ratio_binning.push_back(59);
  ratio_binning.push_back(60);
//  ratio_binning.push_back(70);*/
//  ratio_binning.push_back(80);
//  ratio_binning.push_back(90);
ratio_binning.push_back(80);
//  ratio_binning.push_back(110);
  ratio_binning.push_back(500);
  
  TH1F *observedtb = allsamples.Draw("observedtb", datajzb,  ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,  luminosity);
  TH1F *datapredictiontb = allsamples.Draw("datapredictiontb",    "-"+datajzb, ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,  luminosity);
  TH1F *datapredictiontbo = allsamples.Draw("datapredictiontbo",    "-"+datajzb, ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data,  luminosity);
  datapredictiontb->Add(datapredictiontbo,-1);
  TH1F *analytical_background_predictionb = allsamples.Draw("analytical_background_predictionb",datajzb, ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets&&"mll<2",data,  luminosity);
  for(int i=0;i<=ratio_binning.size();i++) {
    analytical_background_predictionb->SetBinContent(i+1,functions[1]->Eval(analytical_background_predictionb->GetBinCenter(i)));
    analytical_background_predictionb->SetBinError(i+1,TMath::Sqrt(functions[1]->Eval(analytical_background_predictionb->GetBinCenter(i))));
  }
  
  TH1F *ratio = (TH1F*) observedtb->Clone();
  ratio->Divide(datapredictiontb);
  
  for (int i=0;i<=ratio_binning.size();i++) {
    cout << observedtb->GetBinLowEdge(i+1) << ";"<<observedtb->GetBinContent(i+1) << ";" << datapredictiontb->GetBinContent(i+1) << " --> " << ratio->GetBinContent(i+1) << "+/-" << ratio->GetBinError(i+1) << endl;
  }
  
//  ratio->Divide(datapredictiontb);
//  ratio->Divide(analytical_background_predictionb);
//  TGraphAsymmErrors *JZBratio= histRatio(observedtb,analytical_background_predictionb,data,ratio_binning);
//  ratio->Divide(analytical_background_prediction);
//  ratio->Divide(datapredictiont);
//  ratio->GetYaxis()->SetTitle("obs/pred");
//  JZBratio->Draw("AP");
  ratio->GetYaxis()->SetRangeUser(0,10);
  ratio->Draw();
//analytical_background_predictionb->Draw();
//  JZBratio->SetTitle("");
  TText *title4 = write_title("Ratio of observed to predicted");
  title4->Draw();
  
//  CompleteSave(c4,"test/ttbar_discovery_dataprediction___analytical_function");
  CompleteSave(c4,"test/ttbar_discovery_dataprediction__analytical__new_binning_one_huge_bin_from_80");

  
}

void calculate_upper_limits(string mcjzb, string datajzb) {
  write_warning("calculate_upper_limits","Upper limit calculation temporarily deactivated");
//  write_warning("calculate_upper_limits","Calculation of SUSY upper limits has been temporarily suspended in favor of top discovery");
//  rediscover_the_top(mcjzb,datajzb);
/*  
  TCanvas *c3 = new TCanvas("c3","c3");
  c3->SetLogy(1);
  vector<float> binning;
  //binning=allsamples.get_optimal_binsize(mcjzb,cutmass&&cutOSSF&&cutnJets,20,50,800);
  binning.push_back(50);
  binning.push_back(100);
  binning.push_back(150);
  binning.push_back(200);
  binning.push_back(500);
  TH1F *datapredictiona = allsamples.Draw("datapredictiona",    "-"+datajzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc,  luminosity);
  TH1F *datapredictionb = allsamples.Draw("datapredictionb",    "-"+datajzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc,  luminosity);
  TH1F *datapredictionc = allsamples.Draw("datapredictionc",    datajzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc,  luminosity);
  TH1F *dataprediction = (TH1F*)datapredictiona->Clone();
  dataprediction->Add(datapredictionb,-1);
  dataprediction->Add(datapredictionc);
  TH1F *puresignal     = allsamples.Draw("puresignal",        mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
  TH1F *signalpred     = allsamples.Draw("signalpred",    "-"+mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
  TH1F *signalpredlo   = allsamples.Draw("signalpredlo",  "-"+mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
  TH1F *signalpredro   = allsamples.Draw("signalpredro",      mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
  TH1F *puredata       = allsamples.Draw("puredata",          datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
  signalpred->Add(signalpredlo,-1);
  signalpred->Add(signalpredro);
  puresignal->Add(signalpred,-1);//subtracting signal contamination
  ofstream myfile;
  myfile.open ("ShapeFit_log.txt");
  establish_upper_limits(puredata,dataprediction,puresignal,"LM4",myfile);
  myfile.close();
*/ 
}

void susy_scan_axis_labeling(TH2F *histo) {
  histo->GetXaxis()->SetTitle("#Chi_{2}^{0}-LSP");
  histo->GetXaxis()->CenterTitle();
  histo->GetYaxis()->SetTitle("m_{#tilde{q}}");
  histo->GetYaxis()->CenterTitle();
}

void scan_susy_space(string mcjzb, string datajzb) {
  TCanvas *c3 = new TCanvas("c3","c3");
  vector<float> binning;
  binning=allsamples.get_optimal_binsize(mcjzb,cutmass&&cutOSSF&&cutnJets,20,50,800);
  /*
  binning.push_back(50);
  binning.push_back(75);
  binning.push_back(100);
  binning.push_back(150);
  binning.push_back(200);
  binning.push_back(500);
  */
  float arrbinning[binning.size()];
  for(int i=0;i<binning.size();i++) arrbinning[i]=binning[i];
  TH1F *puredata   = allsamples.Draw("puredata",  datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
  puredata->SetMarkerSize(DataMarkerSize);
  TH1F *allbgs   = allsamples.Draw("allbgs",  "-"+datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
  TH1F *allbgsb   = allsamples.Draw("allbgsb",  "-"+datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data,luminosity);
  TH1F *allbgsc   = allsamples.Draw("allbgsc",  datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data,luminosity);
  allbgs->Add(allbgsb,-1);
  allbgs->Add(allbgsc);
  int ndata=puredata->Integral();
  ofstream myfile;
  myfile.open ("susyscan_log.txt");
  TFile *susyscanfile = new TFile("/scratch/fronga/SMS/T5z_SqSqToQZQZ_38xFall10.root");
  TTree *suevents = (TTree*)susyscanfile->Get("events");
  TH2F *exclusionmap = new TH2F("exclusionmap","",20,0,500,20,0,1000);
  TH2F *exclusionmap1s = new TH2F("exclusionmap1s","",20,0,500,20,0,1000);
  TH2F *exclusionmap2s = new TH2F("exclusionmap2s","",20,0,500,20,0,1000);
  TH2F *exclusionmap3s = new TH2F("exclusionmap3s","",20,0,500,20,0,1000);
  
  susy_scan_axis_labeling(exclusionmap);
  susy_scan_axis_labeling(exclusionmap1s);
  susy_scan_axis_labeling(exclusionmap2s);
  susy_scan_axis_labeling(exclusionmap3s);
  
  Int_t MyPalette[100];
  Double_t r[]    = {0., 0.0, 1.0, 1.0, 1.0};
  Double_t g[]    = {0., 0.0, 0.0, 1.0, 1.0};
  Double_t b[]    = {0., 1.0, 0.0, 0.0, 1.0};
  Double_t stop[] = {0., .25, .50, .75, 1.0};
  Int_t FI = TColor::CreateGradientColorTable(5, stop, r, g, b, 100);
  for (int i=0;i<100;i++) MyPalette[i] = FI+i;
   
  gStyle->SetPalette(100, MyPalette);
  
  for(int m23=50;m23<500;m23+=25) {
    for (int m0=(2*(m23-50)+150);m0<=1000;m0+=50)
    {
      c3->cd();
      stringstream drawcondition;
      drawcondition << "pfJetGoodNum>=3&&(TMath::Abs(masses[0]-"<<m0<<")<10&&TMath::Abs(masses[2]-masses[3]-"<<m23<<")<10)&&mll>5&&id1==id2";
      TH1F *puresignal = new TH1F("puresignal","puresignal",binning.size()-1,arrbinning);
      TH1F *puresignall= new TH1F("puresignall","puresignal",binning.size()-1,arrbinning);
      stringstream drawvar,drawvar2;
      drawvar<<mcjzb<<">>puresignal";
      drawvar2<<"-"<<mcjzb<<">>puresignall";
      suevents->Draw(drawvar.str().c_str(),drawcondition.str().c_str());
      suevents->Draw(drawvar2.str().c_str(),drawcondition.str().c_str());
      if(puresignal->Integral()<60) {
        delete puresignal;
        continue;
      }
      puresignal->Add(puresignall,-1);//we need to correct for the signal contamination - we effectively only see (JZB>0)-(JZB<0) !!
      puresignal->Scale(ndata/(20*puresignal->Integral()));//normalizing it to 5% of the data
      stringstream saveas;
      saveas<<"Model_Scan/m0_"<<m0<<"__m23_"<<m23;
      cout << "PLEASE KEEP IN MIND THAT SIGNAL CONTAMINATION IS NOT REALLY TAKEN CARE OF YET DUE TO LOW STATISTICS! SHOULD BE SOMETHING LIKE THIS : "<< endl;
//        TH1F *signalpredlo   = allsamples.Draw("signalpredlo",  "-"+mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
//        TH1F *signalpredro   = allsamples.Draw("signalpredro",      mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
//        TH1F *puredata       = allsamples.Draw("puredata",          datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
//        signalpred->Add(signalpredlo,-1);
//        signalpred->Add(signalpredro);
//        puresignal->Add(signalpred,-1);//subtracting signal contamination
//---------------------
//      cout << "(m0,m23)=("<<m0<<","<<m23<<") contains " << puresignal->Integral() << endl;
//    TH1F *puresignal = allsamples.Draw("puresignal",mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
      vector<float> results=establish_upper_limits(puredata,allbgs,puresignal,saveas.str(),myfile);  
      if(results.size()==0) {
        delete puresignal;
        continue;
      }
      exclusionmap->Fill(m23,m0,results[0]);
      exclusionmap1s->Fill(m23,m0,results[1]);
      exclusionmap2s->Fill(m23,m0,results[2]);
      exclusionmap3s->Fill(m23,m0,results[3]);
      delete puresignal;
      cout << "(m0,m23)=("<<m0<<","<<m23<<") : 3 sigma at " << results[3] << endl;
    }
  }//end of model scan for loop
  
  cout << "Exclusion Map contains" << exclusionmap->Integral() << " (integral) and entries: " << exclusionmap->GetEntries() << endl;
  c3->cd();
  exclusionmap->Draw("CONTZ");
  CompleteSave(c3,"Model_Scan/CONT/Model_Scan_Mean_values");
  exclusionmap->Draw("COLZ");
  CompleteSave(c3,"Model_Scan/COL/Model_Scan_Mean_values");
  
  exclusionmap1s->Draw("CONTZ");
  CompleteSave(c3,"Model_Scan/CONT/Model_Scan_1sigma_values");
  exclusionmap1s->Draw("COLZ");
  CompleteSave(c3,"Model_Scan/COL/Model_Scan_1sigma_values");
  
  exclusionmap2s->Draw("CONTZ");
  CompleteSave(c3,"Model_Scan/CONT/Model_Scan_2sigma_values");
  exclusionmap2s->Draw("COLZ");
  CompleteSave(c3,"Model_Scan/COL/Model_Scan_2sigma_values");
  
  exclusionmap3s->Draw("CONTZ");
  CompleteSave(c3,"Model_Scan/CONT/Model_Scan_3sigma_values");
  exclusionmap3s->Draw("COLZ");
  CompleteSave(c3,"Model_Scan/COL/Model_Scan_3sigma_values");
  
  TFile *exclusion_limits = new TFile("exclusion_limits.root","RECREATE");
  exclusionmap->Write();
  exclusionmap1s->Write();
  exclusionmap2s->Write();
  exclusionmap3s->Write();
  exclusion_limits->Close();
  susyscanfile->Close();
  
  myfile.close();
}

void draw_ttbar_and_zjets_shape_for_one_configuration(string mcjzb, string datajzb, int leptontype=-1, int scenario=0,bool floating=false) {
  //Step 1: Establishing cuts
  stringstream jetcutstring;
  string writescenario="";

  if(scenario==0) jetcutstring << "(pfJetGoodNum>=3)&&"<<(const char*) basicqualitycut;
  if(scenario==1) jetcutstring << "(pfJetPt[0]>50&&pfJetPt[1]>50)&&"<<(const char*)basicqualitycut;
  TCut jetcut(jetcutstring.str().c_str());
  string leptoncut="mll>0";
  if(leptontype==0||leptontype==1) {
    if(leptontype==0) {
      leptoncut="id1==0";
      writescenario="__ee";
    }
    else {
      leptoncut="id1==1";
      writescenario="__ee";
    }
  }
  TCut lepcut(leptoncut.c_str());
  
  TCanvas *c5 = new TCanvas("c5","c5",1500,500);
  TCanvas *c6 = new TCanvas("c6","c6");
  c5->Divide(3,1);

  //STEP 2: Extract Zjets shape in data
  c5->cd(1);
  c5->cd(1)->SetLogy(1);
  TCut massat40("mll>40");
  TH1F *ossfleft   = allsamples.Draw("ossfleft",  "-"+datajzb,40,0,200, "JZB [GeV]", "events", massat40&&cutOSSF&&jetcut&&lepcut,data,luminosity);
  TH1F *osofleft   = allsamples.Draw("osofleft",  "-"+datajzb,40,0,200, "JZB [GeV]", "events", massat40&&cutOSOF&&jetcut&&lepcut,data,luminosity);
  ossfleft->SetLineColor(kRed);
  ossfleft->SetMarkerColor(kRed);
  ossfleft->Add(osofleft,-1); 
  vector<TF1*> functions = do_cb_fit_to_plot(ossfleft,10);
  ossfleft->SetMarkerSize(DataMarkerSize);
  ossfleft->Draw();
  functions[0]->Draw("same");functions[1]->Draw("same");functions[2]->Draw("same");
  TF1 *zjetsfunc = (TF1*) functions[1]->Clone();
  TF1 *zjetsfuncN = (TF1*) functions[0]->Clone();
  TF1 *zjetsfuncP = (TF1*) functions[2]->Clone();
  zjetsfunc->Draw("same");zjetsfuncN->Draw("same");zjetsfuncP->Draw("same");
  TLegend *leg1 = new TLegend(0.6,0.6,0.89,0.80);
  leg1->SetFillColor(kWhite);
  leg1->SetLineColor(kWhite);
  leg1->AddEntry(ossfleft,"OSSF (sub),JZB<peak","p");
  leg1->AddEntry(zjetsfunc,"OSSF fit ('zjets')","l");
  leg1->Draw("same");
  TText *titleleft = write_title("Extracting Z+Jets shape");
  titleleft->Draw();
  
  //Step 3: Extract ttbar shape (in data or MC?)
  c5->cd(2);
  c5->cd(2)->SetLogy(1);
  TH1F *osof;
  TText *titlecenter;
  bool frommc=false;
  if(frommc) {
    osof = allsamples.Draw("osof",datajzb,40,-200,200, "JZB [GeV]", "events", massat40&&cutOSSF&&jetcut&&lepcut,mc,luminosity,allsamples.FindSample("TTJets"));
    titlecenter = write_title("Extracting ttbar shape (from ossf MC)");
  }
  else {
    osof = allsamples.Draw("osof",datajzb,40,-200,200, "JZB [GeV]", "events", massat40&&cutOSOF&&jetcut&&lepcut,data,luminosity);
    titlecenter = write_title("Extracting ttbar shape (from osof data)");
  }
  osof->SetMarkerSize(DataMarkerSize);
  osof->Draw();
  vector<TF1*> ttbarfunctions = do_cb_fit_to_plot(osof,35,true);
  ttbarfunctions[0]->SetLineColor(kRed); ttbarfunctions[0]->SetLineStyle(2);    ttbarfunctions[0]->Draw("same");
  ttbarfunctions[1]->SetLineColor(kRed);                                        ttbarfunctions[1]->Draw("same");
  ttbarfunctions[2]->SetLineColor(kRed); ttbarfunctions[2]->SetLineStyle(2);    ttbarfunctions[2]->Draw("same");

  TLegend *leg2 = new TLegend(0.15,0.8,0.4,0.89);
  leg2->SetFillColor(kWhite);
  leg2->SetLineColor(kWhite);
  if(frommc) {
    leg2->AddEntry(osof,"t#bar{t} OSSF, MC","p");
    leg2->AddEntry(ttbarfunctions[1],"Fit to t#bar{t} OSSF,MC","l");
  } else {
    leg2->AddEntry(osof,"OSOF","p");
    leg2->AddEntry(ttbarfunctions[1],"Fit to OSOF","l");
  }
  leg2->Draw("same");
  titlecenter->Draw();

//--------------------------------------------------------------------------------------------------------------------------------
  //STEP 4: Present it!
  // actually: if we wanna let it float we need to do that first :-) 
  c5->cd(3);
  c5->cd(3)->SetLogy(1);
  TH1F *observed   = allsamples.Draw("observed",datajzb,100,0,500, "JZB [GeV]", "events", massat40&&cutOSSF&&jetcut&&lepcut,data,luminosity);
  observed->SetMarkerSize(DataMarkerSize);

  TF1 *logparc = new TF1("logparc",InvCrystalBall,0,1000,5);    logparc->SetLineColor(kRed);
  TF1 *logparcn = new TF1("logparcn",InvCrystalBallN,0,1000,5); logparcn->SetLineColor(kRed);   logparcn->SetLineStyle(2);
  TF1 *logparcp = new TF1("logparcp",InvCrystalBallP,0,1000,5); logparcp->SetLineColor(kRed);   logparcp->SetLineStyle(2);
  
  TF1 *zjetsc = new TF1("zjetsc",InvCrystalBall,0,1000,5);      zjetsc->SetLineColor(kBlue);    
  TF1 *zjetscn = new TF1("zjetscn",InvCrystalBallN,0,1000,5);   zjetscn->SetLineColor(kBlue);   zjetscn->SetLineStyle(2);
  TF1 *zjetscp = new TF1("zjetscp",InvCrystalBallP,0,1000,5);   zjetscp->SetLineColor(kBlue);   zjetscp->SetLineStyle(2);
  
  TF1 *ZplusJetsplusTTbar = new TF1("ZplusJetsplusTTbar", DoubleInvCrystalBall,0,1000,10);      ZplusJetsplusTTbar->SetLineColor(kBlue);
  TF1 *ZplusJetsplusTTbarP= new TF1("ZplusJetsplusTTbarP",DoubleInvCrystalBallP,0,1000,10);     ZplusJetsplusTTbarP->SetLineColor(kBlue);       ZplusJetsplusTTbarP->SetLineStyle(2);
  TF1 *ZplusJetsplusTTbarN= new TF1("ZplusJetsplusTTbarN",DoubleInvCrystalBallN,0,1000,10);     ZplusJetsplusTTbarN->SetLineColor(kBlue);       ZplusJetsplusTTbarN->SetLineStyle(2);
  
  zjetsc->SetParameters(zjetsfunc->GetParameters());
  zjetscp->SetParameters(zjetsfunc->GetParameters());
  zjetscn->SetParameters(zjetsfunc->GetParameters());
  
  TH1F *observeda   = allsamples.Draw("observeda",datajzb,53,80,350, "JZB [GeV]", "events", massat40&&cutOSSF&&jetcut&&lepcut,data,luminosity);
  //blublu
  logparc->SetParameters(ttbarfunctions[1]->GetParameters());
  logparcn->SetParameters(ttbarfunctions[1]->GetParameters());
  logparcp->SetParameters(ttbarfunctions[1]->GetParameters());
  if(floating) {
    cout << "TTbar contribution assumed (before fitting) : " << logparc->GetParameter(0) << endl;
    logparc->SetParameters(ttbarfunctions[1]->GetParameters());
    for(int i=0;i<10;i++) {
      if(i<5) ZplusJetsplusTTbar->FixParameter(i,zjetsfunc->GetParameter(i));
      if(i>=5) {
        if (i>5) ZplusJetsplusTTbar->FixParameter(i,logparc->GetParameter(i-5));
        if (i==5) ZplusJetsplusTTbar->SetParameter(i,logparc->GetParameter(i-5));
      }
    }//end of setting parameters
    observeda->Draw("same");
    ZplusJetsplusTTbar->Draw("same");
    observeda->Fit(ZplusJetsplusTTbar);
    cout << "--> Quality of Z+Jets / TTbar fit : chi2/ndf = " << ZplusJetsplusTTbar->GetChisquare() << "/" << ZplusJetsplusTTbar->GetNDF() << endl;
    ZplusJetsplusTTbar->Draw("same");
    ZplusJetsplusTTbarP->SetParameters(ZplusJetsplusTTbar->GetParameters());
    ZplusJetsplusTTbarN->SetParameters(ZplusJetsplusTTbar->GetParameters());
    cout << "TTbar contribution found (after fitting) : " << ZplusJetsplusTTbar->GetParameter(5) << endl;
    float factor = ZplusJetsplusTTbar->GetParameter(5) / logparc->GetParameter(0);
    cout << "FACTOR: " << factor << endl; 
    logparc->SetParameter(0,factor*ttbarfunctions[1]->GetParameter(0));
    logparcn->SetParameter(0,factor*ttbarfunctions[1]->GetParameter(0));
    logparcp->SetParameter(0,factor*ttbarfunctions[1]->GetParameter(0));
  }
  
  c5->cd(3);
  c5->cd(3)->SetLogy(1);
  observed->Draw();
  zjetsc->Draw("same");zjetscn->Draw("same");zjetscp->Draw("same");
  logparc->Draw("same");
  logparcn->Draw("same");
  logparcp->Draw("same");

  TLegend *leg3 = new TLegend(0.6,0.6,0.89,0.80);
  leg3->SetFillColor(kWhite);
  leg3->SetLineColor(kWhite);
  leg3->AddEntry(observed,"OSSF,JZB>peak","p");
  leg3->AddEntry(ttbarfunctions[1],"OSOF fit ('ttbar')","l");
  leg3->AddEntry(zjetsfunc,"OSSF,JZB<0 fit ('zjets')","l");
  leg3->Draw("same");
  TText *titleright = write_title("Summary of shapes and observed shape");
  titleright->Draw();
  
  c6->cd()->SetLogy(1);
  observed->Draw();  
  zjetsc->Draw("same");zjetscn->Draw("same");zjetscp->Draw("same");
  logparc->Draw("same");
  logparcn->Draw("same");
  logparcp->Draw("same");
  leg3->Draw("same");
  titleright->Draw();
  
  if(scenario==0) {
    CompleteSave(c5,"Shapes2/Making_of___3jetsabove30"+writescenario);
    CompleteSave(c5->cd(1),"Shapes2/Making_of___3jetsabove30"+writescenario+"__cd1");
    CompleteSave(c5->cd(2),"Shapes2/Making_of___3jetsabove30"+writescenario+"__cd2");
    CompleteSave(c5->cd(3),"Shapes2/Making_of___3jetsabove30"+writescenario+"__cd3");
    CompleteSave(c6,"Shapes2/Background_Shapes___3jetsabove30"+writescenario);
  } else {
    CompleteSave(c5,"Shapes2/Making_of___2jetsabove50"+writescenario);
    CompleteSave(c5->cd(1),"Shapes2/Making_of___2jetsabove50"+writescenario+"__cd1");
    CompleteSave(c5->cd(2),"Shapes2/Making_of___2jetsabove50"+writescenario+"__cd2");
    CompleteSave(c5->cd(3),"Shapes2/Making_of___2jetsabove50"+writescenario+"__cd3");
    CompleteSave(c6,"Shapes2/Background_Shapes___2jetsabove50"+writescenario);
  }
  cout << "Statistics about our fits: " << endl; 
  cout << "Z+Jets shape: Chi2/ndf = " << zjetsfunc->GetChisquare() << "/" << ossfleft->GetNbinsX() << endl;
  cout << "ttbar shape: Chi2/ndf = " << ttbarfunctions[1]->GetChisquare() << "/" << osof->GetNbinsX() << endl;

  c6->cd();
  TLegend *additionallegend = new TLegend(0.6,0.6,0.89,0.89);
  additionallegend->SetFillColor(kWhite);
  additionallegend->SetLineColor(kWhite);
  additionallegend->AddEntry(observed,"Data","p");
  additionallegend->AddEntry(ZplusJetsplusTTbar,"Fitted Z+jets & TTbar","l");
  additionallegend->AddEntry(zjetsc,"Z+jets","l");
  additionallegend->AddEntry(logparc,"TTbar","l");
  observed->Draw();
  ZplusJetsplusTTbar->SetLineColor(kGreen);
  ZplusJetsplusTTbarP->SetLineColor(kGreen);
  ZplusJetsplusTTbarN->SetLineColor(kGreen);
  ZplusJetsplusTTbarP->SetLineStyle(2);
  ZplusJetsplusTTbarN->SetLineStyle(2);
  TF1 *ZplusJetsplusTTbar2 = new TF1("ZplusJetsplusTTbar2",DoubleInvCrystalBall,0,1000,10);
  ZplusJetsplusTTbar2->SetParameters(ZplusJetsplusTTbar->GetParameters());
  ZplusJetsplusTTbar2->SetLineColor(kGreen);
  ZplusJetsplusTTbarP->SetFillColor(TColor::GetColor("#81F781"));
  ZplusJetsplusTTbarN->SetFillColor(kWhite);
  ZplusJetsplusTTbarP->Draw("fcsame");
  ZplusJetsplusTTbarN->Draw("fcsame");
  TH1F *hZplusJetsplusTTbar = (TH1F*)ZplusJetsplusTTbar2->GetHistogram();
  TH1F *hZplusJetsplusTTbarN = (TH1F*)ZplusJetsplusTTbarN->GetHistogram();
  TH1F *hZplusJetsplusTTbarP = (TH1F*)ZplusJetsplusTTbarP->GetHistogram();
  hZplusJetsplusTTbar->SetMarkerSize(0);
  hZplusJetsplusTTbarP->SetMarkerSize(0);
  hZplusJetsplusTTbarN->SetMarkerSize(0);
  for (int i=1;i<=hZplusJetsplusTTbar->GetNbinsX();i++) {
    float newerror=hZplusJetsplusTTbarP->GetBinContent(i)-hZplusJetsplusTTbar->GetBinContent(i);
    hZplusJetsplusTTbar->SetBinError(i,newerror);
    if(hZplusJetsplusTTbar->GetBinContent(i)<0.05) hZplusJetsplusTTbar->SetBinContent(i,0); //avoiding a displaying probolem
  }
  hZplusJetsplusTTbarP->SetFillColor(kGreen);
  hZplusJetsplusTTbarN->SetFillColor(kWhite);
  hZplusJetsplusTTbarN->Draw("same");
  
  ZplusJetsplusTTbar2->SetMarkerSize(0);
  ZplusJetsplusTTbar2->Draw("same");

  zjetsc->Draw("same");zjetscn->Draw("same");zjetscp->Draw("same");
  logparc->Draw("same");
  logparcn->Draw("same");
  logparcp->Draw("same");
  additionallegend->Draw("same");
    if(scenario==0) {
    CompleteSave(c6,"Shapes2/Background_Shapes___3jetsabove30__allfits__"+writescenario);
  } else {
    CompleteSave(c6,"Shapes2/Background_Shapes___2jetsabove50__allfits__"+writescenario);
  }
//--------------------------------------------------------------------------------------------------------------------------------
}


void draw_ttbar_and_zjets_shape(string mcjzb, string datajzb) {
  int all_leptons=-1;
  int electrons_only=0;
  int mu_only=1;
  int twojetswith50gev=1;
  int threejetswith30gev=0;
/*  
  draw_ttbar_and_zjets_shape_for_one_configuration(mcjzb,datajzb,all_leptons,twojetswith50gev);
  draw_ttbar_and_zjets_shape_for_one_configuration(mcjzb,datajzb,all_leptons,threejetswith30gev);
  
  draw_ttbar_and_zjets_shape_for_one_configuration(mcjzb,datajzb,electrons_only,twojetswith50gev);
  draw_ttbar_and_zjets_shape_for_one_configuration(mcjzb,datajzb,electrons_only,threejetswith30gev);

  draw_ttbar_and_zjets_shape_for_one_configuration(mcjzb,datajzb,mu_only,twojetswith50gev);
  draw_ttbar_and_zjets_shape_for_one_configuration(mcjzb,datajzb,mu_only,threejetswith30gev);
  */

  draw_ttbar_and_zjets_shape_for_one_configuration(mcjzb,datajzb,all_leptons,threejetswith30gev,true);
}

void find_correction_factors(string &jzbvardata,string &jzbvarmc) {
  //first: colorful plots
  TCanvas *cancorr = new TCanvas("cancorr","Canvas for Response Correction");
  cancorr->SetLogz();
  gStyle->SetPalette(1);
  TCut zptforresponsepresentation("pt<600&&TMath::Abs(91.2-mll)<20&&id1==id2&&(ch1*ch2<0)&&((sumJetPt[1]/pt)<5.0)");
  TH2F *niceresponseplotd = new TH2F("niceresponseplotd","",400,0,600,400,0,5);
  (allsamples.collection)[allsamples.FindSample("AllData")[0]].events->Draw("sumJetPt[1]/pt:pt>>niceresponseplotd",zptforresponsepresentation);
  niceresponseplotd->SetStats(0);
  niceresponseplotd->GetXaxis()->SetTitle("Z pt [GeV]");
  niceresponseplotd->GetYaxis()->SetTitle("Response");
  niceresponseplotd->GetXaxis()->CenterTitle();
  niceresponseplotd->GetYaxis()->CenterTitle();
  niceresponseplotd->Draw("COLZ");
  TProfile * pro = (TProfile*)niceresponseplotd->ProfileX();
  pro->Rebin(4);
  pro->SetMarkerSize(DataMarkerSize);
  pro->Draw("same,e1");
  DrawPrelim();
  CompleteSave(cancorr,"ResponseCorrection/Response_Correction_Illustration_Data");
  
  (allsamples.collection)[allsamples.FindSample("DY")[0]].events->Draw("sumJetPt[1]/pt:pt>>niceresponseplotd",zptforresponsepresentation);
  niceresponseplotd->SetStats(0);
  niceresponseplotd->GetXaxis()->SetTitle("Z pt [GeV]");
  niceresponseplotd->GetYaxis()->SetTitle("Response");
  niceresponseplotd->Draw("COLZ");
  (allsamples.collection)[allsamples.FindSample("DY")[0]].events->Draw("sumJetPt[1]/pt:pt",zptforresponsepresentation,"PROF,same");
  TProfile * prof = (TProfile*)niceresponseplotd->ProfileX();
//  prof->Rebin(4);
  prof->SetMarkerSize(DataMarkerSize);
  prof->Draw("same,e1");
  DrawPrelim();
  CompleteSave(cancorr,"ResponseCorrection/Response_Correction_Illustration_MC");
  
  //Step 2: Getting the result
  TCut zptcutforresponse("pt>30&&pt<300&&TMath::Abs(91.2-mll)<20&&id1==id2&&(ch1*ch2<0)");
  TH1F *dresponse   = allsamples.Draw("dresponse","sumJetPt[1]/pt",100,0,5, "JZB [GeV]", "events", zptcutforresponse,data,luminosity);
  float datacorrection=dresponse->GetMean();
  stringstream jzbvardatas;
  if(datacorrection>1) jzbvardatas<<"(jzb[1]-"<<datacorrection-1<<"*pt)";
  if(datacorrection<1) jzbvardatas<<"(jzb[1]+"<<1-datacorrection<<"*pt)";
  jzbvardata=jzbvardatas.str();
  TH1F *mcresponse  = allsamples.Draw("mcresponse","sumJetPt[1]/pt",100,0,5, "JZB [GeV]", "events", zptcutforresponse,mc,luminosity,allsamples.FindSample("DYJetsToLL_TuneZ2_M"));
  float mccorrection=mcresponse->GetMean();
  stringstream jzbvarmcs;
  if(mccorrection>1) jzbvarmcs<<"(jzb[1]-"<<mccorrection-1<<"*pt)";
  if(mccorrection<1) jzbvarmcs<<"(jzb[1]+"<<1-mccorrection<<"*pt)";
  jzbvarmc=jzbvarmcs.str();
  cout << "JZB Z pt correction summary : " << endl;
  cout << "  Data: The response is " << dresponse->GetMean() << "  --> jzb variable is now : " << jzbvardata << endl;
  cout << "  MC  : The response is " << mcresponse->GetMean() << "  --> jzb variable is now : " << jzbvarmc << endl;
}

void pick_up_events(string cut) {
  cout << "Picking up events with cut " << cut << endl;
  allsamples.PickUpEvents(cut);
}

void test() {
  TCanvas *testcanv = new TCanvas("testcanv","testcanv");
  testcanv->cd();
  switch_overunderflow(true);
  TH1F *ptdistr   = allsamples.Draw("ptdistr","pt1",100,30,200, "p_{T} [GeV]", "events", cutOSSF,data,luminosity);
  switch_overunderflow(false);
  ptdistr->Draw();
  testcanv->SaveAs("test.png");
  cout << "HELLO there!" << endl;
}
Revision:	1.18
Committed:	Fri Jul 15 06:52:40 2011 UTC (13 years, 10 months ago) by buchmann
Content type:	text/plain
Branch:	MAIN
Changes since 1.17:	+11 -0 lines
Log Message:	Added more plots (mll split up for ee/mm); fixed one potential memory leak