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 <TLegendEntry.h>

using namespace std;

using namespace PlottingSetup;

void todo() {
  dout << "My to do list: " << endl;
  dout << "  - ExperimentalModule::Poisson_ratio_plot : Get the second part to work!" << endl;
  dout << "  - compare_onpeak_offpeak_signal_distributions: Implement this function ... " << endl;
  dout << "Info : The trigger requirement is currently set to " << (const char*) passtrig << endl;
  dout << "Info : The pt requirement is currently set to " << (const char*) passtrig << endl;
  dout << "Info : The mll requirement is currently set to " << (const char*) cutmass << endl;
  dout << "Info : The lepton requirement is currently set to " << (const char*) leptoncut << endl;
  dout << "Info : The weight applied to all MC is " << (const char*) cutWeight << endl;
}  


void find_one_peak_combination(TCut specialcut, float &MCPeak,float &MCPeakError, float &DataPeak, float &DataPeakError, float &MCSigma, float &MCSigmaError, float &DataSigma, float& DataSigmaError, stringstream &result, bool doPUreweighting = true, string saveas="")
{
  // Temporarily switch off PU reweighting, if asked
  TCut weightbackup=cutWeight;
  if ( !doPUreweighting ) cutWeight="1.0";
  
  int nbins=100;
  if(PlottingSetup::DoBTag) nbins=25;
  
  TCanvas *tempcan = new TCanvas("tempcan","Temporary canvas for peak finding preparations");
  TH1F *rawJZBeemmMC      = allsamples.Draw("rawJZBeemmMC",jzbvariablemc,nbins,-50,50, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets&&specialcut,mc, luminosity);
  TH1F *rawJZBeemmData    = allsamples.Draw("rawJZBeemmData",jzbvariabledata,nbins, -50,50, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets&&specialcut,data, luminosity);
  TH1F *rawJZBemMC      = allsamples.Draw("rawJZBemMC",jzbvariablemc,nbins,-50,50, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets&&specialcut,mc, luminosity);
  TH1F *rawJZBemData    = allsamples.Draw("rawJZBemData",jzbvariabledata,nbins, -50,50, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets&&specialcut,data, luminosity);
  TH1F *rawttbarjzbeemmMC;

  if(method==doKM) {
    //we only need this histo for the KM fitting...
    rawttbarjzbeemmMC = allsamples.Draw("rawttbarjzbeemmMC",jzbvariablemc,nbins, -50,50, "JZB [GeV]", "events",cutmass&&cutOSSF&&cutnJets&&specialcut,mc,luminosity,allsamples.FindSample("TTJet"));
    MCPeak=find_peak(rawJZBeemmMC, rawttbarjzbeemmMC, -40, 40, mc, MCPeakError,MCSigma,MCSigmaError,method,saveas);
    DataPeak=find_peak(rawJZBeemmData, rawJZBeemmData, -40, 40, data, DataPeakError,DataSigma,DataSigmaError,method,saveas);
    delete rawttbarjzbeemmMC;
  }
  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,MCSigmaError,method,saveas);
    DataPeak=find_peak(reducedData, rawJZBeemmData, -40, 40, data, DataPeakError,DataSigma,DataSigmaError,method,saveas);
    delete reducedMC;
    delete reducedData;
  }

  // Revert to original PU reweighting
  if ( !doPUreweighting ) cutWeight = weightbackup;
    
//  MCPeak=find_peak(rawJZBeemmMC, rawttbarjzbeemmMC, -40, 40, mc, MCPeakError,MCSigma,method);
//  DataPeak=find_peak(rawJZBeemmData, rawJZBeemmData, -40, 40, data, DataPeakError,DataSigma,method);
  dout   << "   We have found the peak in Data at " << DataPeak << " +/- " << DataPeakError << " with sigma=" << DataSigma << " +/- " << DataSigmaError << endl;
  result << "   We have found the peak in Data at " << DataPeak << " +/- " << DataPeakError << " with sigma=" << DataSigma << " +/- " << DataSigmaError << endl;
  dout   << "   We have found the peak in MC at " << MCPeak << " +/- " << MCPeakError << " with sigma=" << MCSigma << " +/- " << MCSigmaError << endl;
  result << "   We have found the peak in MC at " << MCPeak << " +/- " << MCPeakError << " with sigma=" << MCSigma << " +/- " << MCSigmaError << endl;
  delete rawJZBeemmData;
  delete rawJZBeemmMC;
  delete rawJZBemData;
  delete rawJZBemMC;
  delete tempcan;
}

void find_peaks(float &MCPeak,float &MCPeakError, float &DataPeak, float &DataPeakError, stringstream &result, bool doPUreweighting, stringstream &datajzb, stringstream &mcjzb)
{
  
  bool DoInvidualeemmPeaks=false;
  
  float mcpeak, datapeak;
  float mcpeakerr, datapeakerr;
  
  float mceepeak,mcmmpeak;
  float mceepeakerr,mcmmpeakerr;
  
  float datammpeak,dataeepeak;
  float datammpeakerr,dataeepeakerr;
  
  float mcSigma,mcSigmaError, dataSigma, dataSigmaError;
  
  dout << "Finding global peak : " << endl;
  find_one_peak_combination(TCut(""),mcpeak,mcpeakerr, datapeak,datapeakerr,mcSigma,mcSigmaError, dataSigma,dataSigmaError,result,doPUreweighting,"");
  
  if(DoInvidualeemmPeaks) {
    dout << "Finding peak for electrons : " << endl;
    find_one_peak_combination(TCut("id1==0"),mceepeak,mceepeakerr,dataeepeak,dataeepeakerr,mcSigma,mcSigmaError,dataSigma,dataSigmaError,result,doPUreweighting,"_ele");
    dout << "Finding peak for muons : " << endl;
    find_one_peak_combination(TCut("id1==1"),mcmmpeak,mcmmpeakerr,datammpeak,datammpeakerr,mcSigma,mcSigmaError,dataSigma,dataSigmaError,result,doPUreweighting,"_mu");
    
    datajzb << "(" << jzbvariabledata;
    mcjzb << "(" << jzbvariablemc;
    
    if(dataeepeak>0) datajzb << "- (id1==id2)*(id1==0)*" << TMath::Abs(dataeepeak) << " ";
    else datajzb << "+ (id1==id2)*(id1==0)*" << TMath::Abs(dataeepeak) << " ";
    
    if(datammpeak>0) datajzb << "- (id1==id2)*(id1==1)*" << TMath::Abs(datammpeak) << " ";
    else datajzb << "+ (id1==id2)*(id1==1)*" << TMath::Abs(datammpeak) << " ";
    
    if(datapeak>0) datajzb << "- (id1!=id2)*" << TMath::Abs(datapeak) << " ";
    else datajzb << "+ (id1!=id2)*" << TMath::Abs(datapeak) << " ";
    
    datajzb << ")";
    
    if(mceepeak>0) mcjzb << "- (id1==id2)*(id1==0)*" << TMath::Abs(mceepeak) << " ";
    else mcjzb << "+ (id1==id2)*(id1==0)*" << TMath::Abs(mceepeak) << " ";
    
    if(mcmmpeak>0) mcjzb << "- (id1==id2)*(id1==1)*" << TMath::Abs(mcmmpeak) << " ";
    else mcjzb << "+ (id1==id2)*(id1==1)*" << TMath::Abs(mcmmpeak) << " ";
    
    if(mcpeak>0) mcjzb << "- (id1!=id2)*" << TMath::Abs(mcpeak) << " ";
    else mcjzb << "+ (id1!=id2)*" << TMath::Abs(mcpeak) << " ";
    
    mcjzb << ")";
  } else {
    datajzb << "(" << jzbvariabledata;
      mcjzb << "(" << jzbvariablemc;
    
    if(datapeak>0) datajzb << "- " << TMath::Abs(datapeak) << " ";
    else datajzb << "+ " << TMath::Abs(datapeak) << " ";
    
    datajzb << ")";
    
    if(mcpeak>0) mcjzb << "- " << TMath::Abs(mcpeak) << " ";
    else mcjzb << "+ " << TMath::Abs(mcpeak) << " ";
    
    mcjzb << ")";
  }
  
  
}

void make_special_obs_pred_mll_plot(string datajzb, string mcjzb, float jzbthreshold, float binWidth = 5.0) {
  float min=70.0;
  float max=115.0;
  if(!PlottingSetup::RestrictToMassPeak) {
    min=20;
    max=300;
  }
  int nbins=int((max-min)/binWidth);
  
  TCanvas *ckin = new TCanvas("ckin","Kinematic Plots (in the making)",600,600);
  
  stringstream largerzeroS;
  largerzeroS << "(" << datajzb << ">" << jzbthreshold << ")";
  TCut largerzeroD(largerzeroS.str().c_str());
  
  stringstream smallerzeroS;
  smallerzeroS << "(" << datajzb << "<-" << jzbthreshold << ")";
  TCut smallerzeroD(smallerzeroS.str().c_str());


  stringstream largerzeroMS;
  largerzeroMS << "(" << mcjzb << ">" << jzbthreshold << ")";
  TCut largerzeroM(largerzeroMS.str().c_str());
  
  TH1F *RcorrJZBeemm     = allsamples.Draw("RcorrJZBeemm",       "mll",nbins,min,max, "m_{ll} [GeV]", "events", cutmass&&cutOSSF&&cutnJets&&largerzeroD,data,luminosity);
  THStack mcRcorrJZBeemm = allsamples.DrawStack("mcRcorrJZBeemm","mll",nbins,min,max, "m_{ll} [GeV}", "events", cutmass&&cutOSSF&&cutnJets&&largerzeroM,mc,luminosity);
  TH1F *LcorrJZBeemm     = allsamples.Draw("LcorrJZBeemm",       "mll",nbins,min,max, "m_{ll} [GeV]", "events", cutmass&&cutOSSF&&cutnJets&&smallerzeroD,data,luminosity);
  TH1F *RcorrJZBem       = allsamples.Draw("RcorrJZBem",         "mll",nbins,min,max, "m_{ll} [GeV]", "events", cutmass&&cutOSOF&&cutnJets&&largerzeroD,data,luminosity);
  TH1F *LcorrJZBem       = allsamples.Draw("LcorrJZBem",         "mll",nbins,min,max, "m_{ll} [GeV]", "events", cutmass&&cutOSOF&&cutnJets&&smallerzeroD,data,luminosity);

  TH1F *RcorrJZBSBem;
  TH1F *LcorrJZBSBem;
  TH1F *RcorrJZBSBeemm;
  TH1F *LcorrJZBSBeemm;
  
//  TH1F *RcorrJZBeemmNoS;

  if(PlottingSetup::RestrictToMassPeak) {
    RcorrJZBSBem   = allsamples.Draw("RcorrJZBSBem",  "mll",nbins,min,max, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets&&largerzeroD,data, luminosity);
    LcorrJZBSBem   = allsamples.Draw("LcorrJZBSBem",  "mll",nbins,min,max, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets&&smallerzeroD,data, luminosity);
    RcorrJZBSBeemm = allsamples.Draw("RcorrJZBSBeemm","mll",nbins,min,max, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets&&largerzeroD,data, luminosity);
    LcorrJZBSBeemm = allsamples.Draw("LcorrJZBSBeemm","mll",nbins,min,max, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets&&smallerzeroD,data, luminosity);
  }
  
  // Separate predictions
  TH1F* SFN = (TH1F*)LcorrJZBeemm->Clone("SFN");
  TH1F* OFP = (TH1F*)RcorrJZBem->Clone("OFP");
  TH1F* OFN = (TH1F*)LcorrJZBem->Clone("OFN");
  if(PlottingSetup::RestrictToMassPeak) {
    OFP->Scale(1.0/3.0);
    OFP->Add(RcorrJZBSBem,1.0/3.);
    OFP->Add(RcorrJZBSBeemm,1.0/3.);
    OFN->Scale(1.0/3.0);
    OFN->Add(LcorrJZBSBem,1.0/3.);
    OFN->Add(LcorrJZBSBeemm,1.0/3.);
  } 

  TH1F* Bpred = (TH1F*)SFN->Clone("Bpred");
  Bpred->Add(OFP);
  Bpred->Add(OFN,-1);
  Bpred->SetLineColor(kRed);
  
  RcorrJZBeemm->SetTitleOffset(1.3,"y");
  RcorrJZBeemm->Draw();
  mcRcorrJZBeemm.Draw("same");
  Bpred->Draw("histo,same");
  RcorrJZBeemm->Draw("same");
  
  TLegend *leg = allsamples.allbglegend();
  leg->SetX1(0.58);
  leg->AddEntry(RcorrJZBeemm,"observed (data)","l");
  leg->AddEntry(Bpred,"predicted (data)","l");
  leg->Draw("same");
  
  stringstream saveas;
  saveas << "kin/Mll_After_Cut/Cut_At" << jzbthreshold;
  CompleteSave(ckin,saveas.str());
  
  // Draw all predictions overlayed
  unsigned int w = gStyle->GetHistLineWidth()+1; // Make line a bit wider, since we dash it
  SFN->SetLineColor(kGreen+2);
  SFN->SetLineStyle(2);
  SFN->SetLineWidth(w);
  OFP->SetLineColor(kBlue+2);
  OFP->SetLineStyle(2);
  OFP->SetLineWidth(w);
  OFN->SetLineColor(kMagenta+2);
  OFN->SetLineStyle(3);
  OFN->SetLineWidth(w);
  
  RcorrJZBeemm->Draw();
  SFN->Draw("histo,same");
  OFP->Draw("histo,same");
  OFN->Draw("histo,same");
  Bpred->Draw("histo,same");
  RcorrJZBeemm->Draw("same");
  
  TLegend *leg2 = make_legend("",0.52,0.7);
  leg2->AddEntry(RcorrJZBeemm,"observed (data)","lp");
  leg2->AddEntry(Bpred,"predicted (data)","l");
  leg2->AddEntry(SFN,  "  SF JZB<0","l");
  leg2->AddEntry(OFN,  "  OF JZB<0","l");
  leg2->AddEntry(OFP,  "  OF JZB>0","l");
  leg2->Draw("same");
  
  saveas.str("");
  saveas << "kin/Mll_After_Cut/Cut_At" << jzbthreshold << "_nomc";
  CompleteSave(ckin,saveas.str());
  
  delete RcorrJZBeemm;
  delete LcorrJZBeemm;
  delete RcorrJZBem;
  delete LcorrJZBem;
  if(PlottingSetup::RestrictToMassPeak) {
    delete RcorrJZBSBeemm;
    delete LcorrJZBSBeemm;
    delete RcorrJZBSBem;
    delete LcorrJZBSBem;
  }
  delete Bpred;
  delete ckin;
}

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("datahistoOSOF","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("#FE642E"));
  datahistoOSOF->SetLineColor(kRed);
  datahistoOSOF->SetMarkerStyle(21);
  datahistoOSOF->Draw("same");

  // Try to re-arrange legend...
  TLegend *bgleg = allsamples.allbglegend("",datahistoOSSF);
  TLegend *kinleg = make_legend();
  kinleg->AddEntry(datahistoOSSF,"SF (data)","p");
  kinleg->AddEntry(datahistoOSOF,"OF (data)","p");
  TIter next(bgleg->GetListOfPrimitives());
  TObject* obj;
  // Copy the nice bkgd legend skipping the "data"
  while ( (obj = next()) ) 
    if ( strcmp(((TLegendEntry*)obj)->GetObject()->GetName(),"datahistoOSSF") )
      kinleg->GetListOfPrimitives()->Add(obj);

  kinleg->Draw();
  CompleteSave(ckin,"kin/mll_ossf_osof_distribution");
  
  delete datahistoOSOF;
  delete datahistoOSSF;
  delete ckin;
}
  

void draw_ratio_plot(TH1* hdata, THStack& hmc, float ymin=0.5, float ymax=1.5) {

  // Make a histogram from stack
  TIter next(hmc.GetHists());
  TObject* obj;
  TH1* hratio = NULL;
  while ( (obj = next()) ) {
    if ( !hratio ) {
      hratio = (TH1*)obj->Clone();
      hratio->SetName("hratio");
    } else hratio->Add( (TH1*)obj );
  }
  hratio->Divide(hdata);
  hratio->SetMaximum(ymax);
  hratio->SetMinimum(ymin);
  hratio->SetMarkerStyle(2);
  hratio->SetLineWidth(1);
  hratio->GetYaxis()->SetLabelSize(0.08);
  hratio->GetXaxis()->SetLabelSize(0.0);

  TPad* rpad = new TPad("rpad","",0.15,0.73,0.4,0.88);
  rpad->SetTopMargin(0.0);
  rpad->SetBottomMargin(0.0);
  rpad->SetRightMargin(0.0);
  rpad->Draw();
  rpad->cd();
//  hratio->GetXaxis()->SetNdivisions(0);
  hratio->GetYaxis()->SetNdivisions(502,false);
  hratio->Draw("e1x0");

  TF1* oneline = new TF1("","1.0",0,1000);
  oneline->SetLineColor(kBlue);
  oneline->SetLineStyle(1); 
  oneline->SetLineWidth(1);
  oneline->Draw("same"); 
}

void make_OFSF_plot(string variable, string addcut, string legendTitle, int nbins, float min, float max, bool logscale,
                   string xlabel, string filename, bool isPF=true, bool plotratio=true, bool loadlastminmax=false, float legendPosition=0.3) {

  TCut ibasiccut=basiccut;
  bool draw_separation_lines=false;

  if(isPF) ibasiccut=basiccut&&"pfjzb[0]>-998";

  if(addcut != "") ibasiccut = ibasiccut && addcut.c_str();

  TCut cutSF;
  TCut cutOF;

  cutOF=cutOSOF&&cutnJets&&ibasiccut;
  cutSF=cutOSSF&&cutnJets&&ibasiccut;


  TCanvas *ckin = new TCanvas("ckin","Kinematic Plots (in the making)",600,600);
  ckin->SetLogy(logscale);
  TH1F *datahistoSF = allsamples.Draw("datahistoSF",variable,nbins,min,max, xlabel, "events",cutSF,data,luminosity);
  TH1F *datahistoOF = allsamples.Draw("datahistoOF",variable,nbins,min,max, xlabel, "events",cutOF,data,luminosity);
  string signal("LM3");
  TH1F* signalhisto = new TH1F("signalhisto",signal.c_str(),nbins,min,max);
  int idx = signalsamples.FindSample(signal)[0];
  (signalsamples.collection)[idx].events->Project("signalhisto",variable.c_str(),cutSF);
  signalhisto->Scale((signalsamples.collection)[idx].weight*luminosity);
  signalhisto->SetLineColor((signalsamples.collection)[idx].samplecolor);
  datahistoSF->SetMarkerSize(DataMarkerSize);
  datahistoOF->SetLineColor(kRed);

  if ( !logscale ) {
    datahistoSF->SetMinimum(0); // Defaults
    datahistoOF->SetMinimum(0); // Defaults
  } else {
    datahistoSF->SetMinimum(0.5);
    datahistoOF->SetMinimum(0.5);
  }
  if (logscale) {
    datahistoSF->SetMaximum(5.3*datahistoSF->GetMaximum());
    datahistoOF->SetMaximum(5.3*datahistoSF->GetMaximum());
  } else {
    datahistoSF->SetMaximum(2.0*datahistoSF->GetMaximum());
    datahistoOF->SetMaximum(2.0*datahistoOF->GetMaximum());
  }

  datahistoSF->GetXaxis()->SetTitle(xlabel.c_str());
  datahistoOF->GetXaxis()->SetTitle(xlabel.c_str());
  signalhisto->GetXaxis()->SetTitle(xlabel.c_str());
  datahistoSF->GetYaxis()->SetTitle("N. Events");
  datahistoOF->GetYaxis()->SetTitle("N. Events");
  signalhisto->GetYaxis()->SetTitle("N. Events");
  datahistoSF->GetXaxis()->CenterTitle();
  datahistoOF->GetXaxis()->CenterTitle();
  signalhisto->GetXaxis()->CenterTitle();
  datahistoSF->GetYaxis()->CenterTitle();
  datahistoOF->GetYaxis()->CenterTitle();
  signalhisto->GetYaxis()->CenterTitle();

  TLegend *mleg = new TLegend(0.2+legendPosition, 0.72, 0.55+legendPosition, 0.87, legendTitle.c_str());
  mleg->SetFillColor(kWhite);
  mleg->SetTextFont(42);
  mleg->SetTextSize(0.048);
  mleg->SetLineWidth(0);
  mleg->SetBorderSize(0);
  mleg->AddEntry(datahistoSF, "S. Flavor", "PL");
  mleg->AddEntry(datahistoOF, "O. Flavor", "L");
  mleg->AddEntry(signalhisto, "LM3", "L");

  datahistoSF->Draw("E1");
  datahistoOF->Draw("HIST,SAMES");
  signalhisto->Draw("HIST,SAMES");
  mleg->Draw();
  DrawPrelim();
  CompleteSave(ckin, "./SFOF/" + filename);

  datahistoSF->Delete();
  datahistoOF->Delete();
  signalhisto->Delete();
  delete mleg;
  delete ckin;

}


float lastrange_min=0;
float lastrange_max=0;

void make_kin_plot(string variable, string addcut, int nbins, float min, float max, bool logscale,
                   string xlabel, string filename, bool isPF=true, bool plotratio=true, bool loadlastminmax=false ) {
//  TCut basiccut("(pfJetGoodNum>=2&&pfJetGoodID[0])&&(pfJetGoodNum>=2&&pfJetGoodID[1])&&(passed_triggers||!is_data)");
  TCut ibasiccut=basiccut;
  bool draw_separation_lines=false;
  
  if(isPF) ibasiccut=basiccut&&"pfjzb[0]>-998";

  if(addcut != "") ibasiccut = ibasiccut && addcut.c_str();

  //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=="pt2") variable="pfpt2[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
  TCut cut;
  cut=cutmass&&cutOSSF&&cutnJets&&ibasiccut;
  if(filename=="nJets") cut=cutmass&&cutOSSF&&ibasiccut;
  if(filename=="nJets_osof") cut=cutmass&&cutOSOF&&ibasiccut;
  if(filename=="nJets_nocuts_except_mll_ossf") cut=cutmass&&cutOSSF;
  if(filename=="mll") {
        cut=cutOSSF&&cutnJets&&ibasiccut;
        draw_separation_lines=true;
  }
  if(filename=="mll_ee") cut=cutOSSF&&cutnJets&&ibasiccut&&"id1==0";
  if(filename=="mll_osof") {
        cut=cutOSOF&&cutnJets&&ibasiccut;
        draw_separation_lines=true;
  }
  if(filename=="mll_mm") cut=cutOSSF&&cutnJets&&ibasiccut&&"id1==1";
  if(Contains(filename,"aboveJZB")) cut=cutOSSF&&cutnJets&&ibasiccut;
  if(Contains(filename,"mll_ee_above")) cut=cut&&"id1==0";
  if(Contains(filename,"mll_mm_above")) cut=cut&&"id1==1";
  if(Contains(filename,"mll_osof_aboveJZB")) cut=cutOSOF&&cutnJets&&ibasiccut;
  if(filename=="mll_inclusive"||filename=="mll_inclusive_highrange") cut=cutOSSF;
  if(filename=="mll_inclusive_osof") cut=cutOSOF;
  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&&cutnJets;
  if(filename=="pfJetGoodPt_0") cut=cutOSSF&&cutmass&&ibasiccut&&cutnJets;

  TCanvas *ckin = new TCanvas("ckin","Kinematic Plots (in the making)",600,600);
  ckin->SetLogy(logscale);
  TH1F *datahisto = allsamples.Draw("datahisto",variable,nbins,min,max, xlabel, "events",cut,data,luminosity);
  datahisto->SetMarkerSize(DataMarkerSize);
  if ( !logscale ) datahisto->SetMinimum(0); // Defaults
  else datahisto->SetMinimum(0.5);
  // Custom max.
  if(variable=="TMath::Abs(pfJetDphiMet[0])") datahisto->SetMaximum(1.5*datahisto->GetMaximum());
  if(variable=="pfJetGoodPt[0]") datahisto->SetMaximum(10*datahisto->GetMaximum());
  if(variable=="pt") datahisto->SetMaximum(10*datahisto->GetMaximum());
  if(filename=="mll_inclusive"||filename=="mll_inclusive_mm"||filename=="mll_inclusive_ee") datahisto->SetMinimum(1);
  if(filename=="mll_osof") datahisto->SetMaximum(10*datahisto->GetMaximum());
  if(filename=="mll_osof") datahisto->SetMinimum(9);
  if (logscale) datahisto->SetMaximum(5.3*datahisto->GetMaximum());
  else datahisto->SetMaximum(1.3*datahisto->GetMaximum());

  cout << "******** Cut used : " << (const char*) cut << endl;
  if(loadlastminmax) {
        datahisto->SetMinimum(lastrange_min);
        datahisto->SetMaximum(lastrange_max);
        if(logscale) {
                datahisto->SetMinimum(pow(10,lastrange_min));
                datahisto->SetMaximum(pow(10,lastrange_max));
        }               
  }

  // Draw signal by hand (for some reason I don't manage to use the sample class: it adds it to the stack!)
  string signal("LM3");
  TH1F* signalhisto = new TH1F("signalhisto",signal.c_str(),nbins,min,max);
  int idx = signalsamples.FindSample(signal)[0];
  (signalsamples.collection)[idx].events->Project("signalhisto",variable.c_str(),cut);
  signalhisto->Scale((signalsamples.collection)[idx].weight*luminosity);
  signalhisto->SetLineColor(kOrange);

  THStack mcstack   = allsamples.DrawStack("mcstack",  variable,nbins,min,max,xlabel,"events",cut,mc,luminosity);
  datahisto->Draw("e1");
  ckin->Update();
  mcstack.Draw("same");
  
  datahisto->Draw("same,e1");
  TLegend *kinleg = allsamples.allbglegend();
  kinleg->Draw();
  if(filename=="mll_osof") {
        kinleg->SetHeader("Opposite flavor");
        kinleg->SetX1(0.58);
  }
  if(filename=="mll") {
        kinleg->SetHeader("Same flavor");
        kinleg->SetX1(0.58);
  }  
  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);

  TLine *lowerboundary;
  TLine *upperboundary;

  if(RestrictToMassPeak&&draw_separation_lines) {
        Color_t linecolor=kBlue;
        float linemin=pow(10,ckin->GetUymin());
        if(filename=="mll_osof") linemin=pow(10,lastrange_min);
        lowerboundary = new TLine(71,linemin,71,datahisto->GetMaximum());
        upperboundary = new TLine(111,linemin,111,datahisto->GetMaximum());
        lowerboundary->SetLineColor(linecolor);
        lowerboundary->SetLineStyle(2);
        upperboundary->SetLineColor(linecolor);
        upperboundary->SetLineStyle(2);
  }

  lastrange_min=ckin->GetUymin();
  lastrange_max=ckin->GetUymax();


  if ( plotratio ) {
        TPad *kinpad = new TPad("kinpad","kinpad",0,0,1,1);
        kinpad->cd();
        kinpad->SetLogy(logscale);
        datahisto->Draw("e1");
        mcstack.Draw("same");
        signalhisto->Draw("same");
        datahisto->Draw("same,e1");
        datahisto->Draw("same,axis");
        if(RestrictToMassPeak&&draw_separation_lines) {
          lowerboundary->Draw("same");
          upperboundary->Draw("same");
        }
          
        kinleg->AddEntry("signalhisto",signal.c_str(),"l");
        kinleg->Draw();
        write_cut->Draw();
        DrawPrelim();
        string saveas="kin/"+filename;
        if(isPF) saveas="kin/"+filename+"__PF";
        save_with_ratio(datahisto,mcstack,kinpad->cd(),saveas);
//      if(isPF) CompleteSave(with_ratio,"kin/"+filename+"__PF_withratio");
//      else CompleteSave(with_ratio,"kin/"+filename+"_withratio");
//      delete with_ratio;
  } else {
        if(isPF) CompleteSave(ckin,"kin/"+filename+"__PF");
        else CompleteSave(ckin,"kin/"+filename);
  }
  datahisto->Delete();
  delete ckin;
}

void make_JES_plot(TCut cut, string name) {

  int nbins=10;
  float min=-0.5;
  float max = 9.5;
  bool logscale=true;
  string xlabel="nJets";

  TCanvas *ckin = new TCanvas("ckin","Kinematic Plots (in the making)",600,600);
  ckin->SetLogy(logscale);
  TH1F *datahisto = allsamples.Draw("datahisto","pfJetGoodNum40",nbins,min,max, xlabel, "events",cut,data,luminosity);
  datahisto->SetMarkerSize(DataMarkerSize);
  THStack mcstack = allsamples.DrawStack("mcstack","pfJetGoodNum40",nbins,min,max, xlabel, "events",cut,mc,luminosity);
  TH1F *JESup = allsamples.Draw("JESup","pfJetGoodNum40p1sigma",nbins,min,max, xlabel, "events",cut,mc,luminosity);
  TH1F *JESdn = allsamples.Draw("JESdn","pfJetGoodNum40n1sigma",nbins,min,max, xlabel, "events",cut,mc,luminosity);

  datahisto->SetMinimum(1);
  datahisto->SetMaximum(5.3*datahisto->GetMaximum()); // in line with kinematic plots style

  float xs[nbins],ys[nbins],exs[nbins],eys[nbins];
  for(int i=1;i<JESup->GetNbinsX();i++) {
    float up=JESup->GetBinContent(i);
    float dn=JESdn->GetBinContent(i);
    xs[i]=JESup->GetBinCenter(i);
    ys[i]=0.5*(up+dn);
    exs[i]=0.5*JESup->GetBinWidth(i);
    eys[i]=0.5*TMath::Abs(up-dn);
  }

  TGraphAsymmErrors *JESunc = new TGraphAsymmErrors(nbins, xs,ys,exs,exs,eys,eys);
  JESunc->SetFillColor(TColor::GetColor("#00ADE1"));
  JESunc->SetFillStyle(3002);
  datahisto->Draw("e1");
  mcstack.Draw("same");
  JESunc->Draw("2");
  datahisto->Draw("same,e1");
  TLegend *kinleg = allsamples.allbglegend();
  kinleg->AddEntry(JESunc,"JES uncertainty","f");
  kinleg->Draw();
  CompleteSave(ckin,"Systematics/JES"+name);
  datahisto->Delete();
  delete ckin;

}

void do_kinematic_plots(string mcjzb, string datajzb, bool doPF=false)
{
  bool dolog=true;
  bool nolog=false;
  if(doPF) write_warning(__FUNCTION__,"Please use caution when trying to produce PF plots; not all versions of the JZB trees have these variables!");
  float mll_low=50;
  float mll_hi=160;
  if(!PlottingSetup::RestrictToMassPeak) {
        mll_low=20;
        mll_hi=300;
  }

  /*
  make_OFSF_plot("pfJetGoodNum40", "mll>20&&mll<40&&met[4]>100", "20 GeV < mll < 40 GeV", 5, 3, 8, true, "N. jets", "njetsLow", false, false, false, 0.2);
  make_OFSF_plot("pfJetGoodNum40", "mll>40&&mll<70&&met[4]>100", "40 GeV < mll < 70 GeV", 5, 3, 8, true, "N. jets", "njetsMed", false, false, false, 0.2);
  make_OFSF_plot("pfJetGoodPt[0]/pfJetGoodPt[1]", "mll>20&&mll<40&&met[4]>100", "20 GeV < mll < 40 GeV", 20, 1, 10, true, "pt_{j}^{1}/pt_{j}^{2}", "jpt1pt2Low", false, false, false, 0.2);
  make_OFSF_plot("pfJetGoodPt[0]/pfJetGoodPt[1]", "mll>40&&mll<70&&met[4]>100", "40 GeV < mll < 70 GeV", 20, 1, 10, true, "pt_{j}^{1}/pt_{j}^{2}", "jpt1pt2Med", false, false, false, 0.2);
  make_OFSF_plot("TMath::Abs(pfJetDphiMet[0])", "mll>20&&mll<40&&met[4]>100", "20 GeV < mll < 40 GeV", 16, 0, 3.2, false, "|#Delta#phi(jet1,MET)|", "dphijetmetLow", false, false, false, 0.0);
  make_OFSF_plot("TMath::Abs(pfJetDphiMet[0])", "mll>40&&mll<70&&met[4]>100", "40 GeV < mll < 70 GeV", 16, 0, 3.2, false, "|#Delta#phi(jet1,MET)|", "dphijetmetMed", false, false, false, 0.0);
  make_OFSF_plot("TMath::Abs(dphi)", "mll>20&&mll<40&&met[4]>100", "20 GeV < mll < 40 GeV", 16, 0, 3.2, false, "|#Delta#phi(l1,l2)|", "dphiLow", false, false, false, 0.2);
  make_OFSF_plot("TMath::Abs(dphi)", "mll>40&&mll<70&&met[4]>100", "40 GeV < mll < 70 GeV", 16, 0, 3.2, false, "|#Delta#phi(l1,l2)|", "dphiMed", false, false, false, 0.2);
  make_OFSF_plot("TMath::Abs(dphiMet1)", "mll>20&&mll<40&&met[4]>100", "20 GeV < mll < 40 GeV", 16, 0, 3.2, false, "|#Delta#phi(l1,MET)|", "dphiMet1Low", false, false, false, 0.2);
  make_OFSF_plot("TMath::Abs(dphiMet1)", "mll>40&&mll<70&&met[4]>100", "40 GeV < mll < 70 GeV", 16, 0, 3.2, false, "|#Delta#phi(l1,MET)|", "dphiMet1Med", false, false, false, 0.2);
  make_OFSF_plot("TMath::Abs(dphiMet2)", "mll>20&&mll<40&&met[4]>100", "20 GeV < mll < 40 GeV", 16, 0, 3.2, false, "|#Delta#phi(l2,MET)|", "dphiMet2Low", false, false, false, 0.2);
  make_OFSF_plot("TMath::Abs(dphiMet2)", "mll>40&&mll<70&&met[4]>100", "40 GeV < mll < 70 GeV", 16, 0, 3.2, false, "|#Delta#phi(l2,MET)|", "dphiMet2Med", false, false, false, 0.2);
  make_OFSF_plot("TMath::Min(TMath::Abs(dphiMet1), TMath::Abs(dphiMet2))", "mll>20&&mll<40&&met[4]>100", "20 GeV < mll < 40 GeV", 16, 0, 3.2, false, "Min(|#Delta#phi(l,MET)|)", "dphilcLow", false, false, false, 0.2);
  make_OFSF_plot("TMath::Min(TMath::Abs(dphiMet1), TMath::Abs(dphiMet2))", "mll>40&&mll<70&&met[4]>100", "40 GeV < mll < 70 GeV", 16, 0, 3.2, false, "Min(|#Delta#phi(l,MET)|)", "dphilcMed", false, false, false, 0.2);
  make_OFSF_plot("TMath::Min(TMath::Abs(pfJetDphiMet[0]), TMath::Min(TMath::Abs(pfJetDphiMet[1]), TMath::Abs(pfJetDphiMet[2])))", "mll>20&&mll<40&&met[4]>100", "20 GeV < mll < 40 GeV", 16, 0, 3.2, false, "Min(|#Delta#phi(jet,MET)|)", "dphijcLow", false, false, false, 0.2);
  make_OFSF_plot("TMath::Min(TMath::Abs(pfJetDphiMet[0]), TMath::Min(TMath::Abs(pfJetDphiMet[1]), TMath::Abs(pfJetDphiMet[2])))", "mll>20&&mll<40&&met[4]>100", "40 GeV < mll < 70 GeV", 16, 0, 3.2, false, "Min(|#Delta#phi(jet,MET)|)", "dphijcMed", false, false, false, 0.2);
  make_OFSF_plot("TMath::Min((TMath::Pi()-TMath::Abs(dphiMet1)), (TMath::Pi() - TMath::Abs(dphiMet2)))", "mll>20&&mll<40&&met[4]>100", "20 GeV < mll < 40 GeV", 16, 0, 3.2, false, "Min(#Pi - |#Delta#phi(l,MET)|)", "dphilcoLow", false, false, false, 0.2);
  make_OFSF_plot("TMath::Min((TMath::Pi()-TMath::Abs(dphiMet1)), (TMath::Pi() - TMath::Abs(dphiMet2)))", "mll>20&&mll<40&&met[4]>100", "40 GeV < mll < 70 GeV", 16, 0, 3.2, false, "Min(#Pi - |#Delta#phi(l,MET)|)", "dphilcoMed", false, false, false, 0.2);
  make_OFSF_plot("TMath::Min((TMath::Pi()-TMath::Abs(pfJetDphiMet[0])), TMath::Min( (TMath::Pi()-TMath::Abs(pfJetDphiMet[1])), (TMath::Pi() - TMath::Abs(pfJetDphiMet[2]))))", "mll>20&&mll<40&&met[4]>100", "20 GeV < mll < 40 GeV", 16, 0, 3.2, false, "Min(#Pi - |#Delta#phi(jet,MET)|)", "dphijcoLow", false, false, false, 0.2);
  make_OFSF_plot("TMath::Min((TMath::Pi()-TMath::Abs(pfJetDphiMet[0])), TMath::Min( (TMath::Pi()-TMath::Abs(pfJetDphiMet[1])), (TMath::Pi() - TMath::Abs(pfJetDphiMet[2]))))", "mll>20&&mll<40&&met[4]>100", "40 GeV < mll < 70 GeV", 16, 0, 3.2, false, "Min(#Pi - |#Delta#phi(jet,MET)|)", "dphijcoMed", false, false, false, 0.2);


  make_kin_plot("met[4]","",70,0,350,dolog,"MET [GeV]","met",doPF,true);

  make_kin_plot("mll","",(int)((mll_hi-mll_low)/5),mll_low,mll_hi,dolog,"m_{ll} [GeV]","mll",doPF,true);
  make_kin_plot("mll","",(int)((mll_hi-mll_low)/5),mll_low,mll_hi,dolog,"m_{ll} [GeV]","mll_osof",doPF,true,true);
  make_kin_plot("mll","",(int)((mll_hi-mll_low)/5),mll_low,mll_hi,dolog,"m_{ll} [GeV]","mll_ee",doPF,true);
  make_kin_plot("mll","",(int)((mll_hi-mll_low)/5),mll_low,mll_hi,dolog,"m_{ll} [GeV]","mll_mm",doPF,true);
  make_kin_plot("mll","",(int)((mll_hi-mll_low))/5,mll_low,mll_hi,dolog,"m_{ll} [GeV]","mll_inclusive",doPF,true);
  make_kin_plot("mll","",(int)((mll_hi-mll_low))/5,mll_low,mll_hi,dolog,"m_{ll} [GeV]","mll_inclusive_ee",doPF,true);
  make_kin_plot("mll","",(int)((mll_hi-mll_low))/5,mll_low,mll_hi,dolog,"m_{ll} [GeV]","mll_inclusive_mm",doPF,true);
  make_kin_plot("mll","",(int)((mll_hi-mll_low))/5,mll_low,mll_hi,dolog,"m_{ll} [GeV]","mll_inclusive_osof",doPF,true);
  make_kin_plot("mll","",(int)((350-mll_low))/5,mll_low,350,dolog,"m_{ll} [GeV]","mll_inclusive_highrange",doPF);
  if(!doPF) make_special_mll_plot((int)((mll_hi-mll_low)/5),mll_low,mll_hi,dolog,"m_{ll} [GeV]");

  make_kin_plot("pfJetGoodPt[0]/pfJetGoodPt[1]","",45,1,10,dolog,"pt_{j}^{1}/pt_{j}^{2}","j1j2ratio",doPF,true);
  make_kin_plot("TMath::Abs(pfJetDphiMet[0])","",32,0,3.2,nolog,"|#Delta#phi(jet1,MET)|","dphiJ1MET",doPF,true);

  make_kin_plot("pfJetGoodNum40","",9,-0.5,8.5,dolog,"nJets","nJets",doPF);
  make_kin_plot("pfJetGoodNum40","",9,-0.5,8.5,dolog,"nJets","nJets_osof",doPF);
  make_kin_plot("pfJetGoodNum40","",9,-0.5,8.5,dolog,"nJets","nJets_nocuts_except_mll_ossf",doPF);

  make_kin_plot("numVtx","",(int)(30.5-(-0.5)),-0.5,30.5,nolog,"N(Vtx)","numVtx",doPF);
//  make_kin_plot("jetpt[0]","",40,0,200,dolog,"leading jet p_{T} [GeV]","pfJetGoodPt_0",doPF);
//  make_kin_plot("jeteta[0]","",40,-5,5,nolog,"leading jet #eta","pfJetGoodEta_0",doPF);
  make_kin_plot("pt","",50,0,500,dolog,"Z p_{T} [GeV]","Zpt",doPF);
  make_kin_plot("pt1","",50,0,200,nolog,"p_{T} [GeV]","pt1",doPF);
  make_kin_plot("pt2","",50,0,200,nolog,"p_{T} [GeV]","pt2",doPF);
  make_kin_plot("eta1","",40,-3,3,nolog,"#eta_{l}","eta",doPF);
  make_kin_plot("jzb[1]","",100,-150,200,dolog,"JZB [GeV]","jzb_ossf",doPF);
  stringstream jzbcut;
  jzbcut << "((is_data&&("<<datajzb<<")>100)||(!is_data&&("<<mcjzb<<")>100))";
  make_kin_plot("mll",jzbcut.str(),(int)((mll_hi-mll_low)/5),mll_low,mll_hi,nolog,"m_{ll} [GeV]","mll_aboveJZB100",doPF,true);
  make_kin_plot("mll",jzbcut.str(),(int)((mll_hi-mll_low)/5),mll_low,mll_hi,nolog,"m_{ll} [GeV]","mll_osof_aboveJZB100",doPF,true);
  make_kin_plot("mll",jzbcut.str(),(int)((mll_hi-mll_low)/5),mll_low,mll_hi,nolog,"m_{ll} [GeV]","mll_ee_aboveJZB100",doPF,true);
  make_kin_plot("mll",jzbcut.str(),(int)((mll_hi-mll_low)/5),mll_low,mll_hi,nolog,"m_{ll} [GeV]","mll_mm_aboveJZB100",doPF,true);
  stringstream jzbcut2;
  jzbcut2 << "((is_data&&("<<datajzb<<")>150)||(!is_data&&("<<mcjzb<<")>150))";
  make_kin_plot("mll",jzbcut2.str(),(int)((mll_hi-mll_low)/5),mll_low,mll_hi,nolog,"m_{ll} [GeV]","mll_aboveJZB150",doPF,true);
  make_kin_plot("mll",jzbcut2.str(),(int)((mll_hi-mll_low)/5),mll_low,mll_hi,nolog,"m_{ll} [GeV]","mll_osof_aboveJZB150",doPF,true);
  make_kin_plot("mll",jzbcut2.str(),(int)((mll_hi-mll_low)/5),mll_low,mll_hi,nolog,"m_{ll} [GeV]","mll_ee_aboveJZB150",doPF,true);
  make_kin_plot("mll",jzbcut2.str(),(int)((mll_hi-mll_low)/5),mll_low,mll_hi,nolog,"m_{ll} [GeV]","mll_mm_aboveJZB150",doPF,true);
  stringstream jzbcut3;
  jzbcut3 << "((is_data&&("<<datajzb<<")>50)||(!is_data&&("<<mcjzb<<")>50))";
  make_kin_plot("mll",jzbcut3.str(),(int)((mll_hi-mll_low)/5),mll_low,mll_hi,nolog,"m_{ll} [GeV]","mll_aboveJZB50",doPF,true);   
  make_kin_plot("mll",jzbcut3.str(),(int)((mll_hi-mll_low)/5),mll_low,mll_hi,nolog,"m_{ll} [GeV]","mll_osof_aboveJZB50",doPF,true,true);
  make_kin_plot("mll",jzbcut3.str(),(int)((mll_hi-mll_low)/5),mll_low,mll_hi,nolog,"m_{ll} [GeV]","mll_ee_aboveJZB50",doPF,true);   
  make_kin_plot("mll",jzbcut3.str(),(int)((mll_hi-mll_low)/5),mll_low,mll_hi,nolog,"m_{ll} [GeV]","mll_mm_aboveJZB50",doPF,true);   

  make_kin_plot("mll","met[4]>100",(int)((mll_hi-mll_low)/5),mll_low,mll_hi,nolog,"m_{ll} [GeV] (MET>100GeV)","mll_met100",doPF,true);

  make_special_obs_pred_mll_plot(datajzb,mcjzb,0);
  make_special_obs_pred_mll_plot(datajzb,mcjzb,50);
  make_special_obs_pred_mll_plot(datajzb,mcjzb,80);
  make_special_obs_pred_mll_plot(datajzb,mcjzb,100);
//   make_special_obs_pred_mll_plot(datajzb,mcjzb,150);
//   make_special_obs_pred_mll_plot(datajzb,mcjzb,200);
//   make_special_obs_pred_mll_plot(datajzb,mcjzb,250);

  make_JES_plot(cutmass&&cutOSSF&&basiccut,"_ossf");
  make_JES_plot(cutmass&&cutOSOF&&basiccut,"_osof");
  */
  
  string regioncut[4] =    {"mll>0&&","mll>20&&mll<70&&","mll>70&&mll<110&&","mll>110&&"};
  string regionname[4] = {"global","lowmass","Zregion","highmass"};
  bool StoreRatioPlots=false;
  for(int iregion=0;iregion<4;iregion++) {
    cout << regionname[iregion] << " : " << regioncut[iregion] << endl;

    //Z properties
    make_kin_plot("pt","met[4]>100&&id1==0&&"+regioncut[iregion]+(const char*)TCut(cutmass&&cutOSSF&&cutnJets),40,0,200,dolog,"p_{T}^{ee} [GeV] (MET>100GeV)","Z/ZPT_ee_met100_"+regionname[iregion],doPF,StoreRatioPlots);
    make_kin_plot("pt","met[4]>100&&id1==1&&"+regioncut[iregion]+(const char*)TCut(cutmass&&cutOSSF&&cutnJets),40,0,200,dolog,"p_{T}^{#mu#mu} [GeV] (MET>100GeV)","Z/ZPT_mm_met100_"+regionname[iregion],doPF,StoreRatioPlots);
    
    make_kin_plot("phi","met[4]>100&&id1==0&&"+regioncut[iregion]+(const char*)TCut(cutmass&&cutOSSF&&cutnJets),32,-3.2,3.2,nolog,"#phi^{ee} (MET>100GeV)","Z/ZPhi_ee_met100_"+regionname[iregion],doPF,StoreRatioPlots);
    make_kin_plot("phi","met[4]>100&&id1==1&&"+regioncut[iregion]+(const char*)TCut(cutmass&&cutOSSF&&cutnJets),32,-3.2,3.2,nolog,"#phi^{#mu#mu} (MET>100GeV)","Z/ZPhi_mm_met100_"+regionname[iregion],doPF,StoreRatioPlots);
    
    make_kin_plot("eta","met[4]>100&&id1==0&&"+regioncut[iregion]+(const char*)TCut(cutmass&&cutOSSF&&cutnJets),30,-3.0,3.0,nolog,"#eta^{ee} (MET>100GeV)","Z/ZEta_ee_met100_"+regionname[iregion],doPF,StoreRatioPlots);
    make_kin_plot("eta","met[4]>100&&id1==1&&"+regioncut[iregion]+(const char*)TCut(cutmass&&cutOSSF&&cutnJets),30,-3.0,3.0,nolog,"#eta^{#mu#mu} (MET>100GeV)","Z/ZEta_mm_met100_"+regionname[iregion],doPF,StoreRatioPlots);
    
    // lepton properties
    make_kin_plot("pt1","met[4]>100&&id1==0&&"+regioncut[iregion]+(const char*)TCut(cutmass&&cutnJets),40,0,200,dolog,"p_{T}^{e1} [GeV] (MET>100GeV)","Lepton/LeadingElectronPt_met100_"+regionname[iregion],doPF,StoreRatioPlots);
    make_kin_plot("pt1","met[4]>100&&id1==1&&"+regioncut[iregion]+(const char*)TCut(cutmass&&cutnJets),40,0,200,dolog,"p_{T}^{#mu1} [GeV] (MET>100GeV)","Lepton/LeadingMuonPt_met100_"+regionname[iregion],doPF,StoreRatioPlots);
    make_kin_plot("pt2","met[4]>100&&id2==0&&"+regioncut[iregion]+(const char*)TCut(cutmass&&cutnJets),40,0,200,dolog,"p_{T}^{e2} [GeV] (MET>100GeV)","Lepton/SubLeadingElectronPt_met100_"+regionname[iregion],doPF,StoreRatioPlots);
    make_kin_plot("pt2","met[4]>100&&id2==1&&"+regioncut[iregion]+(const char*)TCut(cutmass&&cutnJets),40,0,200,dolog,"p_{T}^{#mu2} [GeV] (MET>100GeV)","Lepton/SubLeadingMuonPt_met100_"+regionname[iregion],doPF,StoreRatioPlots);
    
    make_kin_plot("phi1","met[4]>100&&id1==0&&"+regioncut[iregion]+(const char*)TCut(cutmass&&cutnJets),32,-3.2,3.2,nolog,"#phi^{e1} (MET>100GeV)","Lepton/LeadingElectronPhi_met100_"+regionname[iregion],doPF,StoreRatioPlots);
    make_kin_plot("phi1","met[4]>100&&id1==1&&"+regioncut[iregion]+(const char*)TCut(cutmass&&cutnJets),32,-3.2,3.2,nolog,"#phi^{#mu1} (MET>100GeV)","Lepton/LeadingMuonPhi_met100_"+regionname[iregion],doPF,StoreRatioPlots);
    make_kin_plot("phi2","met[4]>100&&id2==0&&"+regioncut[iregion]+(const char*)TCut(cutmass&&cutnJets),32,-3.2,3.2,nolog,"#phi^{e2} (MET>100GeV)","Lepton/SubLeadingElectronPhi_met100_"+regionname[iregion],doPF,StoreRatioPlots);
    make_kin_plot("phi2","met[4]>100&&id2==1&&"+regioncut[iregion]+(const char*)TCut(cutmass&&cutnJets),32,-3.2,3.2,nolog,"#phi^{#mu2} (MET>100GeV)","Lepton/SubLeadingMuonPhi_met100_"+regionname[iregion],doPF,StoreRatioPlots);
    
    make_kin_plot("eta1","met[4]>100&&id1==0&&"+regioncut[iregion]+(const char*)TCut(cutmass&&cutnJets),30,-3.0,3.0,nolog,"#eta^{e1} (MET>100GeV)","Lepton/LeadingElectronEta_met100_"+regionname[iregion],doPF,StoreRatioPlots);
    make_kin_plot("eta1","met[4]>100&&id1==1&&"+regioncut[iregion]+(const char*)TCut(cutmass&&cutnJets),30,-3.0,3.0,nolog,"#eta^{#mu1} (MET>100GeV)","Lepton/LeadingMuonEta_met100_"+regionname[iregion],doPF,StoreRatioPlots);
    make_kin_plot("eta2","met[4]>100&&id2==0&&"+regioncut[iregion]+(const char*)TCut(cutmass&&cutnJets),30,-3.0,3.0,nolog,"#eta^{e2} (MET>100GeV)","Lepton/SubLeadingElectronEta_met100_"+regionname[iregion],doPF,StoreRatioPlots);
    make_kin_plot("eta2","met[4]>100&&id2==1&&"+regioncut[iregion]+(const char*)TCut(cutmass&&cutnJets),30,-3.0,3.0,nolog,"#eta^{#mu2} (MET>100GeV)","Lepton/SubLeadingMuonEta_met100_"+regionname[iregion],doPF,StoreRatioPlots);
    
    
  }

}

void make_comp_plot( string var, string xlabel, string filename, float jzbcut, string mcjzb, string datajzb,
                     int nbins, float xmin, float xmax, bool log,
                     float ymin=0, float ymax=0, bool leftJustified=false ) {
  flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak ---- WATCH OUT: the argument in the function changed!

  TCut weightbackup=cutWeight;//backing up the correct weight (restoring below!)
  if(weightbackup==TCut("1.0")||weightbackup==TCut("1")) write_warning(__FUNCTION__,"WATCH OUT THE WEIGHT HAS POSSIBLY NOT BEEN RESET!!!! PLEASE CHANGE LINE "+any2string(__LINE__));
  //if(var=="numVtx") cutWeight=TCut("1.0");
  TCut jzbData[]= { TCut(TString(datajzb+">"+any2string(jzbcut))),TCut(TString(datajzb+"<-"+any2string(jzbcut))) };
  TCut jzbMC[]  = { TCut(TString(mcjzb+">"+any2string(jzbcut))),TCut(TString(mcjzb+"<-"+any2string(jzbcut))) };

  flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak ---- below: the next ~20 lines changed!
  int nRegions=4;
  if(!PlottingSetup::RestrictToMassPeak) {
    nRegions=2;
  }
  
  string sRegions[] = { "SFZP","OFZP","SFSB","OFSB" };
  TCut kRegions[] = { cutOSSF&&cutnJets&&cutmass,     cutOSOF&&cutnJets&&cutmass,
                      cutOSSF&&cutnJets&&sidebandcut, cutOSOF&&cutnJets&&sidebandcut };
  
  //find ymax
  TH1F *Refdatahisto = allsamples.Draw("datahisto",   var,nbins,xmin,xmax,xlabel,"events",kRegions[0]&&jzbData[0],data,luminosity);
  ymax=int((Refdatahisto->GetMaximum()+2*TMath::Sqrt(Refdatahisto->GetMaximum()))*1.2+0.5);
  delete Refdatahisto;
                      
  for ( int iregion=0; iregion<nRegions; ++iregion ) 
    for ( int ijzb=0; ijzb<2; ++ijzb ) {
      TCanvas *ccomp = new TCanvas("ccomp","Comparison plot",600,400);
      ccomp->SetLogy(log);
      TH1F *datahisto = allsamples.Draw("datahisto",   var,nbins,xmin,xmax,xlabel,"events",kRegions[iregion]&&jzbData[ijzb],data,luminosity);
      TH1F *lm3histo  = signalsamples.Draw("lm3histo",   var,nbins,xmin,xmax,xlabel,"events",kRegions[iregion]&&jzbMC[ijzb],data,luminosity,signalsamples.FindSample("LM3"));
      THStack mcstack = allsamples.DrawStack("mcstack",var,nbins,xmin,xmax,xlabel,"events",kRegions[iregion]&&jzbMC[ijzb],  mc,  luminosity);
      datahisto->SetMarkerSize(DataMarkerSize);
      if (ymax>ymin) datahisto->SetMaximum(ymax);
      lm3histo->SetLineStyle(2);
      datahisto->Draw("e1");
      mcstack.Draw("same");
      datahisto->Draw("same,e1");
      lm3histo->Draw("hist,same");
      TLegend *kinleg = allsamples.allbglegend((sRegions[iregion]+(ijzb?"neg":"pos")).c_str());
      if ( leftJustified ) {
        Float_t w = kinleg->GetX2()-kinleg->GetX1();
        kinleg->SetX1(0.2);
        kinleg->SetX2(0.2+w);
      }
      kinleg->AddEntry(lm3histo,"LM3","l");
      kinleg->Draw();
      TText* write_variable = write_text(0.99,0.01,var);
      write_variable->SetTextAlign(31);
      write_variable->SetTextSize(0.02);
      ccomp->RedrawAxis();
      CompleteSave(ccomp,"compare/JZBcut_at_"+any2string(jzbcut)+"/"+filename+"/"+filename+sRegions[iregion]+(ijzb?"neg":"pos"));
      delete datahisto;
      delete ccomp;
      delete lm3histo;
    }
  cutWeight=weightbackup;
}


void region_comparison_plots(string mcjzb, string datajzb, vector<float> jzb_cuts) {
  dout << "Creating comparison plots for signal and control regions" << endl;
  // Compare a few quantities in the signal region and all 7 control regions

//  switch_overunderflow(true);  // switching overflow/underflow bins on
    
  
    flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak ---- the arguments changed
  for(int ijzb=0;ijzb<(int)jzb_cuts.size();ijzb++) {
    float jzbcut=jzb_cuts[ijzb]; // Comparison plots are done for this JZB cut
    float mll_low=50;float mll_high=170;
    if(!PlottingSetup::RestrictToMassPeak) {
        mll_high=300;
        mll_low=20;
    }
    make_comp_plot("pfJetGoodPt[0]/pfJetGoodPt[1]","pt_{j}^{1}/pt_{j}^{2}","j1j2ratio",jzbcut,mcjzb,datajzb,100,0,10,true);
    make_comp_plot("TMath::Abs(pfJetDphiMet[0])","|#Delta#phi(jet1,MET)|","dphiJ1MET",jzbcut,mcjzb,datajzb,32,0,3.2,false,0,0,true);

    make_comp_plot("mll","m_{ll} [GeV]","mll",jzbcut,mcjzb,datajzb,56,mll_low,mll_high,false,0,16.);
    make_comp_plot("met[4]","pfMET [GeV]","pfmet",jzbcut,mcjzb,datajzb,18,0,360,false,0,16.);
    make_comp_plot("pfJetGoodNum40","#(jets)","njets",jzbcut,mcjzb,datajzb,10,0,10, false,0,35.);
    make_comp_plot("pfJetGoodNumBtag","#(b-jets)","nBjets",jzbcut,mcjzb,datajzb,10,0,10, false,0,35.);
    make_comp_plot("pt","Z p_{T} [GeV]","Zpt",jzbcut,mcjzb,datajzb,26,0,525,false,0.,21.);
    make_comp_plot("numVtx","#(prim. vertices)","nvtx",jzbcut,mcjzb,datajzb,40,0.,40.,false,0,16.);
    make_comp_plot("TMath::Abs(dphi)","#Delta#phi(leptons)","dphilep",jzbcut,mcjzb,datajzb,10,0.,3.1415,false,0,16.,true);
    make_comp_plot("TMath::Abs(dphi_sumJetVSZ[1])","#Delta#phi(Z,jets)","dphiZjets",jzbcut,mcjzb,datajzb,10,0.,3.1415,false,0,16.,true);
    make_comp_plot("eta1","#eta_1","eta1",jzbcut,mcjzb,datajzb,10,0.,2.5,false,0,16.);
    make_comp_plot("eta2","#eta_2","eta2",jzbcut,mcjzb,datajzb,10,0.,2.5,false,0,16.);
  }

  switch_overunderflow(false); // switching overflow/underflow bins off
}


void do_kinematic_PF_plots(string mcjzb, string datajzb)
{
  do_kinematic_plots(mcjzb,datajzb,true);
}

void signal_bg_comparison()
{
  TCanvas *can = new TCanvas("can","Signal Background Comparison Canvas");
  can->SetLogy(1);

  int sbg_nbins=130;
  float sbg_min=-500; //-110;
  float sbg_max=800; //jzbHigh;
  
  float simulatedlumi=luminosity;//in pb please - adjust to your likings
  
  TH1F *JZBplotZJETs = allsamples.Draw("JZBplotZJETs",jzbvariablemc,sbg_nbins,sbg_min,sbg_max, "JZB [GeV]", "events",cutmass&&cutOSSF&&cutnJets,mc,simulatedlumi,allsamples.FindSample("/DY"));
  TH1F *JZBplotLM4;
  if(PlottingSetup::RestrictToMassPeak) JZBplotLM4 = allsamples.Draw("JZBplotLM4",jzbvariablemc,sbg_nbins,sbg_min,sbg_max, "JZB [GeV]", "events",cutmass&&cutOSSF&&cutnJets,mc,simulatedlumi,allsamples.FindSample("LM4"));
  else JZBplotLM4 = allsamples.Draw("JZBplotLM4",jzbvariablemc,sbg_nbins,sbg_min,sbg_max, "JZB [GeV]", "events",cutmass&&cutOSSF&&cutnJets,mc,simulatedlumi,allsamples.FindSample("LM3"));
  TH1F *JZBplotTtbar = allsamples.Draw("JZBplotTtbar",jzbvariablemc,sbg_nbins,sbg_min,sbg_max, "JZB [GeV]", "events",cutmass&&cutOSSF&&cutnJets,mc,simulatedlumi,allsamples.FindSample("TTJets"));
  
  JZBplotTtbar->SetLineColor(allsamples.GetColor("TTJet"));
  JZBplotZJETs->SetFillColor(allsamples.GetColor("DY"));
  JZBplotZJETs->SetLineColor(kBlack);
  JZBplotLM4->SetLineStyle(2);
  JZBplotZJETs->SetMaximum(JZBplotZJETs->GetMaximum()*5);
  JZBplotZJETs->SetMinimum(1);

  JZBplotTtbar->SetMaximum(JZBplotZJETs->GetMaximum());
  JZBplotTtbar->SetMinimum(0.01);
  JZBplotTtbar->SetFillColor(allsamples.GetColor("TTJets"));
  JZBplotTtbar->DrawClone("histo");
  JZBplotZJETs->Draw("histo,same");
  JZBplotTtbar->SetFillColor(0);
  JZBplotTtbar->DrawClone("histo,same");
  JZBplotTtbar->SetFillColor(allsamples.GetColor("TTJets"));
  JZBplotLM4->Draw("histo,same");


  TLegend *signal_bg_comparison_leg2 =  make_legend("",0.55,0.75,false);
  signal_bg_comparison_leg2->AddEntry(JZBplotZJETs,"Z+Jets","f");
  signal_bg_comparison_leg2->AddEntry(JZBplotTtbar,"t#bar{t}","f");
  if(PlottingSetup::RestrictToMassPeak) signal_bg_comparison_leg2->AddEntry(JZBplotLM4,"LM4","f");
  else signal_bg_comparison_leg2->AddEntry(JZBplotLM4,"LM3","f");
  signal_bg_comparison_leg2->Draw();
  DrawMCPrelim(simulatedlumi);
  CompleteSave(can,"jzb_bg_vs_signal_distribution");

  // Define illustrative set of SMS points
  TCut kSMS1("MassGlu==250&&MassLSP==75");
  TCut kSMS2("MassGlu==800&&MassLSP==200");
  TCut kSMS3("MassGlu==1050&&MassLSP==850");
  TCut kSMS4("MassGlu==1200&&MassLSP==100");
  
  //If the scan samples haven't been loaded yet, this is a good point to load them (all of them!)
  if((scansample.collection).size()<2) define_SMS_sample(false, allsamples, signalsamples, scansample, true); // loading ALL zones for the scans, not only the basic one.

  
  TH1F *JZBplotSMS1 = scansample.Draw("JZBplotSMS1",jzbvariablemc,sbg_nbins,sbg_min,sbg_max, "JZB [GeV]", "events",cutmass&&cutOSSF&&cutnJets&&kSMS1,mc,simulatedlumi,scansample.FindSample("t"));
  JZBplotSMS1->Scale(JZBplotLM4->Integral()/JZBplotSMS1->Integral());

  TH1F *JZBplotSMS2 = scansample.Draw("JZBplotSMS2",jzbvariablemc,sbg_nbins,sbg_min,sbg_max, "JZB [GeV]", "events",cutmass&&cutOSSF&&cutnJets&&kSMS2,mc,simulatedlumi,scansample.FindSample("t"));
  JZBplotSMS2->Scale(JZBplotLM4->Integral()/JZBplotSMS2->Integral());

  TH1F *JZBplotSMS3 = scansample.Draw("JZBplotSMS3",jzbvariablemc,sbg_nbins,sbg_min,sbg_max, "JZB [GeV]", "events",cutmass&&cutOSSF&&cutnJets&&kSMS3,mc,simulatedlumi,scansample.FindSample("t"));
  JZBplotSMS3->Scale(JZBplotLM4->Integral()/JZBplotSMS3->Integral());
  
  TH1F *JZBplotSMS4 = scansample.Draw("JZBplotSMS4",jzbvariablemc,sbg_nbins,sbg_min,sbg_max, "JZB [GeV]", "events",cutmass&&cutOSSF&&cutnJets&&kSMS4,mc,simulatedlumi,scansample.FindSample("t"));
  JZBplotSMS4->Scale(JZBplotLM4->Integral()/JZBplotSMS4->Integral());

  // Draw all plots overlaid
  JZBplotTtbar->SetMaximum(JZBplotZJETs->GetMaximum());
  JZBplotTtbar->SetMinimum(0.01);
  JZBplotTtbar->SetFillColor(allsamples.GetColor("TTJets"));
  JZBplotTtbar->DrawClone("histo");
  JZBplotZJETs->Draw("histo,same");
  JZBplotTtbar->SetFillColor(0);
  JZBplotTtbar->DrawClone("histo,same");
  JZBplotTtbar->SetFillColor(allsamples.GetColor("TTJets"));

  JZBplotSMS1->SetLineColor(kRed+1);
  JZBplotSMS2->SetLineColor(kBlue+1);
  JZBplotSMS3->SetLineColor(kRed+1);
  JZBplotSMS4->SetLineColor(kBlue+1);
  JZBplotSMS3->SetLineStyle(2);
  JZBplotSMS4->SetLineStyle(2);
  
  JZBplotSMS1->Draw("histo,same");
  JZBplotSMS2->Draw("histo,same");
  JZBplotSMS3->Draw("histo,same");
  JZBplotSMS4->Draw("histo,same");
  JZBplotLM4->SetLineColor(kGreen);JZBplotLM4->Draw("histo,same");
  TLegend *signal_bg_comparison_leg6 = make_legend("",0.55,0.55,false);
  signal_bg_comparison_leg6->AddEntry(JZBplotZJETs,"Z+Jets","f");
  signal_bg_comparison_leg6->AddEntry(JZBplotTtbar,"t#bar{t}","f");
  signal_bg_comparison_leg6->AddEntry(JZBplotSMS1,"","");
  signal_bg_comparison_leg6->AddEntry(JZBplotSMS1,"SMS parameters","");
  signal_bg_comparison_leg6->AddEntry(JZBplotSMS1,"(250,75) [GeV]","f");
  signal_bg_comparison_leg6->AddEntry(JZBplotSMS2,"(800,200) [GeV]","f");
  signal_bg_comparison_leg6->AddEntry(JZBplotSMS3,"(1050,850) [GeV]","f");
  signal_bg_comparison_leg6->AddEntry(JZBplotSMS4,"(1200,100) [GeV]","f");
  signal_bg_comparison_leg6->AddEntry(JZBplotLM4,"LM4","f");
  signal_bg_comparison_leg6->Draw();
  DrawMCPrelim(simulatedlumi);
  CompleteSave(can,"jzb_bg_vs_signal_distribution_SMS__summary");
    
  while((scansample.collection).size() > 1) scansample.RemoveLastSample();
    
}

vector<TF1*> do_cb_fit_to_plot(TH1F *histo, float Sigma, float doingfitacrosstheboard=false) {
  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);
  
  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->SetLineStyle(2);
  
  BpredFuncN->SetParameters(BpredFunc->GetParameters());
  BpredFuncN->SetLineColor(kBlue);
  BpredFuncN->SetLineStyle(2);
  
  vector<TF1*> functions;
  functions.push_back(BpredFuncN);
  functions.push_back(BpredFunc);
  functions.push_back(BpredFuncP);
  return functions;
}


TF1* do_logpar_fit_to_plot(TH1F *osof) {
  TCanvas *logpar_fit_can = new TCanvas("logpar_fit_can","Fit canvas for LogPar");
  TF1 *logparfunc = new TF1("logparfunc",LogParabola,0,300,3);
  TF1 *logparfunc2 = new TF1("logparfunc2",LogParabola,0,(osof->GetXaxis()->GetBinLowEdge(osof->GetNbinsX())+osof->GetXaxis()->GetBinWidth(osof->GetNbinsX())),3);
  TF1 *logparfuncN = new TF1("logparfuncN",LogParabolaN,0,(osof->GetXaxis()->GetBinLowEdge(osof->GetNbinsX())+osof->GetXaxis()->GetBinWidth(osof->GetNbinsX())),3);
  TF1 *logparfuncP = new TF1("logparfuncP",LogParabolaP,0,(osof->GetXaxis()->GetBinLowEdge(osof->GetNbinsX())+osof->GetXaxis()->GetBinWidth(osof->GetNbinsX())),3);
  osof->SetMinimum(0);
  osof->Fit(logparfunc,"QR");
  osof->Draw();
  logparfunc->SetLineWidth(2);
  logparfunc2->SetParameters(logparfunc->GetParameters());
  logparfuncN->SetParameters(logparfunc->GetParameters());
  logparfuncP->SetParameters(logparfunc->GetParameters());
  stringstream fitinfo;
  fitinfo << "#Chi^{2} / NDF : " << logparfunc->GetChisquare() << " / " << logparfunc->GetNDF();
  TText *writefitinfo = write_text(0.8,0.8,fitinfo.str());
  writefitinfo->SetTextSize(0.03);
  DrawPrelim();
  writefitinfo->Draw();
  logparfunc->Draw("same");
  logparfunc2->Draw("same");
  logparfuncN->SetLineStyle(2);
  logparfuncP->SetLineStyle(2);
  logparfuncN->Draw("same");
  logparfuncP->Draw("same");
  CompleteSave(logpar_fit_can,"MakingOfBpredFunction/Bpred_Data_LogPar_Fit_To_TTbarPred");
  delete logpar_fit_can;
  return logparfunc2;
}

vector<TF1*> do_extended_fit_to_plot(TH1F *prediction, TH1F *Tprediction, TH1F *ossf, TH1F *osof,int isdata) {
  /* there are mainly two background contributions: Z+Jets (a) and ttbar (b). So:
  a) The case is clear - we take the OSSF prediction - OSOF prediction, and fit a crystal ball function to it. We then extract the CB parameters.
  b) For ttbar, we use the OSOF distribution and look at the [10,100] GeV JZB range, and fit our log parabola. We then extract the LP parameters.
  Once we have these two components, we use the combined parameters to get the final function and we're done.
  */
  //Step 1: take the OSSF prediction - OSOF prediction, and fit a crystal ball function to it
  TH1F *step1cb = (TH1F*)ossf->Clone("step1cb");
  step1cb->Add(osof,-1);
  vector<TF1*> functions = do_cb_fit_to_plot(step1cb,PlottingSetup::JZBPeakWidthData);
  TF1 *zjetscrystalball = functions[1];
  
  //Step 2: use the OSOF distribution and look at the [10,100] GeV JZB range, and fit our log parabola
//  TH1F *ttbarprediction=(TH1F*)prediction->Clone("ttbarprediction");
//  ttbarprediction->Add(ossf,-1);//without the Z+Jets estimate, this is really just the ttbar estimate!
//      the line above is not necessary anymore as we're now looking at a prediction without Z+Jets, and not multiplied with (1.0/3)
  TF1 *ttbarlogpar = do_logpar_fit_to_plot(Tprediction);
  flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
  if(PlottingSetup::RestrictToMassPeak) ttbarlogpar->SetParameter(0,1.0/3*ttbarlogpar->GetParameter(0));//correcting for the fact that we didn't multiply with (1.0/3);


  TF1 *ttbarlogparP = new TF1("ttbarlogparP",LogParabolaP,0,(prediction->GetXaxis()->GetBinLowEdge(prediction->GetNbinsX())+prediction->GetXaxis()->GetBinWidth(prediction->GetNbinsX())),8);
  TF1 *ttbarlogparN = new TF1("ttbarlogparN",LogParabolaN,0,(prediction->GetXaxis()->GetBinLowEdge(prediction->GetNbinsX())+prediction->GetXaxis()->GetBinWidth(prediction->GetNbinsX())),8);
  
  //and now fuse the two!
  TF1 *kmlp = new TF1("kmlp", CrystalBallPlusLogParabola, 0,(prediction->GetXaxis()->GetBinLowEdge(prediction->GetNbinsX())+prediction->GetXaxis()->GetBinWidth(prediction->GetNbinsX())),8);
  TF1 *kmlpP= new TF1("kmlpP",CrystalBallPlusLogParabolaP,0,(prediction->GetXaxis()->GetBinLowEdge(prediction->GetNbinsX())+prediction->GetXaxis()->GetBinWidth(prediction->GetNbinsX())),8);
  TF1 *kmlpN= new TF1("kmlpN",CrystalBallPlusLogParabolaN,0,(prediction->GetXaxis()->GetBinLowEdge(prediction->GetNbinsX())+prediction->GetXaxis()->GetBinWidth(prediction->GetNbinsX())),8);
  double kmlp_pars[10];
  for(int i=0;i<5;i++) kmlp_pars[i]=zjetscrystalball->GetParameter(i);
  for(int i=0;i<3;i++) kmlp_pars[5+i]=ttbarlogpar->GetParameter(i);
  ttbarlogparP->SetParameters(ttbarlogpar->GetParameters());
  ttbarlogparN->SetParameters(ttbarlogpar->GetParameters());
  kmlp->SetParameters(kmlp_pars);
  prediction->Fit(kmlp,"Q");//fitting the final result (done this in the past but kicked it)
/*
if you want to start from scratch (without the partial fitting and only fitting the whole thing, some good start values could be : 
  */
  kmlp_pars[0]=kmlp->GetParameter(0);
  kmlp_pars[1]=3.6198;
  kmlp_pars[2]=16.4664;
  kmlp_pars[3]=1.92253;
  kmlp_pars[4]=3.56099;
  kmlp_pars[5]=5.83;
  kmlp_pars[6]=0.000757479;
  kmlp_pars[7]=95.6157;
  kmlp_pars[8]=0;
  kmlp_pars[9]=0;
  kmlp->SetParameters(kmlp_pars);
/**/
  prediction->Fit(kmlp,"Q");//fitting the final result (done this in the past but kicked it)

  kmlpP->SetParameters(kmlp->GetParameters());
  kmlpN->SetParameters(kmlp->GetParameters());
  
  // now that we're done, let's save all of this so we can have a look at it afterwards.
  TCanvas *can = new TCanvas("can","Prediction Fit Canvas");
  can->SetLogy(1);
  prediction->SetMarkerColor(kRed);
  prediction->Draw();
  
  kmlp->SetLineColor(TColor::GetColor("#04B404"));
  kmlpP->SetLineColor(TColor::GetColor("#04B404"));
  kmlpN->SetLineColor(TColor::GetColor("#04B404"));
  kmlp->Draw("same");
  kmlpN->SetLineStyle(2);
  kmlpP->SetLineStyle(2);
  kmlpN->Draw("same");
  kmlpP->Draw("same");
  
  ttbarlogpar->SetLineColor(TColor::GetColor("#CC2EFA"));//purple
  ttbarlogpar->Draw("same");
  ttbarlogparP->SetLineColor(TColor::GetColor("#CC2EFA"));//purple
  ttbarlogparN->SetLineColor(TColor::GetColor("#CC2EFA"));//purple
  ttbarlogparP->SetLineStyle(2);
  ttbarlogparN->SetLineStyle(2);
  ttbarlogparP->Draw("same");
  ttbarlogparN->Draw("same");
  
  functions[0]->Draw("same");functions[1]->Draw("same");functions[2]->Draw("same");

  TLegend *analyticalBpredLEG = make_legend("",0.5,0.55);
  analyticalBpredLEG->AddEntry(prediction,"predicted","p");
  analyticalBpredLEG->AddEntry(functions[1],"Crystal Ball fit","l");
  analyticalBpredLEG->AddEntry(functions[0],"1#sigma Crystal Ball fit","l");
  analyticalBpredLEG->AddEntry(ttbarlogparN,"TTbar fit","l");
  analyticalBpredLEG->AddEntry(ttbarlogpar,"1#sigma TTbar fit","l");
  analyticalBpredLEG->AddEntry(kmlp,"Combined function","l");
  analyticalBpredLEG->AddEntry(kmlpN,"1#sigma combined function","l");
  analyticalBpredLEG->Draw("same");
  
  if(isdata==0) CompleteSave(can,"MakingOfBpredFunction/Bpred_MC_Analytical_Function_Composition");
  if(isdata==1) CompleteSave(can,"MakingOfBpredFunction/Bpred_data_Analytical_Function_Composition");
  if(isdata==2) CompleteSave(can,"MakingOfBpredFunction/Bpred_MCBnS_Analytical_Function_Composition");
  delete can;
  
  //and finally: prep return functions
  vector<TF1*> return_functions;
  return_functions.push_back(kmlpN);
  return_functions.push_back(kmlp);
  return_functions.push_back(kmlpP);
  
  return_functions.push_back(ttbarlogparN);
  return_functions.push_back(ttbarlogpar);
  return_functions.push_back(ttbarlogparP);
  
  return_functions.push_back(functions[0]);
  return_functions.push_back(functions[1]);
  return_functions.push_back(functions[2]);
  
  return return_functions;
}

void do_prediction_plot(string jzb, TCanvas *globalcanvas, float high, int use_data, bool overlay_signal = false,string subdir="" )
{
//  switch_overunderflow(true);
  bool is_data=false;
  bool use_signal=false;
  if(use_data==1) is_data=true;
  if(use_data==2) use_signal=true;
  int nbins=50;//100;
  if(is_data) nbins=50;
  float low=0;
  float hi=500;
  
  TH1F *blankback = new TH1F("blankback","blankback",int(high/10),0,high);
  TH1F *RcorrJZBeemm   = allsamples.Draw("RcorrJZBeemm",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,is_data, luminosity,use_signal);
  TH1F *LcorrJZBeemm   = allsamples.Draw("LcorrJZBeemm",("-"+jzb).c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,is_data, luminosity,use_signal);
  TH1F *RcorrJZBem     = allsamples.Draw("RcorrJZBem",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,is_data, luminosity,use_signal);
  TH1F *LcorrJZBem     = allsamples.Draw("LcorrJZBem",("-"+jzb).c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,is_data, luminosity,use_signal);

  blankback->GetXaxis()->SetTitle(RcorrJZBeemm->GetXaxis()->GetTitle());
  blankback->GetYaxis()->SetTitle(RcorrJZBeemm->GetYaxis()->GetTitle());
  blankback->GetXaxis()->CenterTitle();
  blankback->GetYaxis()->CenterTitle();
  
  flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak ---- prediction changed.
  TH1F *RcorrJZBSBem;
  TH1F *LcorrJZBSBem;
  TH1F *RcorrJZBSBeemm;
  TH1F *LcorrJZBSBeemm;
  
  TH1F *RcorrJZBeemmNoS;

    //these are for the ratio
  
  TH1F *JRcorrJZBeemm   = allsamples.Draw("JRcorrJZBeemm",jzb.c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,is_data, luminosity,use_signal);
  TH1F *JLcorrJZBeemm   = allsamples.Draw("JLcorrJZBeemm",("-"+jzb).c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,is_data, luminosity,use_signal);
  TH1F *JRcorrJZBem     = allsamples.Draw("JRcorrJZBem",jzb.c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,is_data, luminosity,use_signal);
  TH1F *JLcorrJZBem     = allsamples.Draw("JLcorrJZBem",("-"+jzb).c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,is_data, luminosity,use_signal);
  
  TH1F *JRcorrJZBSBem;
  TH1F *JLcorrJZBSBem;
  TH1F *JRcorrJZBSBeemm;
  TH1F *JLcorrJZBSBeemm;
  
  if(use_data==2 || overlay_signal) RcorrJZBeemmNoS = allsamples.Draw("RcorrJZBeemmNoS",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,is_data, luminosity,false);

  
  if(PlottingSetup::RestrictToMassPeak) {
    RcorrJZBSBem   = allsamples.Draw("RcorrJZBSBem",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets,is_data, luminosity,use_signal);
    LcorrJZBSBem   = allsamples.Draw("LcorrJZBSBem",("-"+jzb).c_str(),nbins,low,hi, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets,is_data, luminosity,use_signal);
    RcorrJZBSBeemm = allsamples.Draw("RcorrJZBSBeemm",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets,is_data, luminosity,use_signal);
    LcorrJZBSBeemm = allsamples.Draw("LcorrJZBSBeemm",("-"+jzb).c_str(),nbins,low,hi, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets,is_data, luminosity,use_signal);
    
    //these are for the ratio
    JRcorrJZBSBem   = allsamples.Draw("JRcorrJZBSBem",jzb.c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets,is_data, luminosity,use_signal);
    JLcorrJZBSBem   = allsamples.Draw("JLcorrJZBSBem",("-"+jzb).c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets,is_data, luminosity,use_signal);
    JRcorrJZBSBeemm = allsamples.Draw("JRcorrJZBSBeemm",jzb.c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets,is_data, luminosity,use_signal);
    JLcorrJZBSBeemm = allsamples.Draw("JLcorrJZBSBeemm",("-"+jzb).c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets,is_data, luminosity,use_signal);
  }
  
  TH1F *lm4RcorrJZBeemm;
  if(overlay_signal || use_data == 2 || use_data == 1) lm4RcorrJZBeemm = allsamples.Draw("lm4RcorrJZBSBeemm",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,is_data, luminosity,allsamples.FindSample("LM"));
  
  flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak ---- prediction changed.
  
  TH1F *Zjetspred = (TH1F*)LcorrJZBeemm->Clone("Zjetspred");
  TH1F *TTbarpred = (TH1F*)RcorrJZBem->Clone("TTbarpred");
    
  TH1F *Bpred = (TH1F*)LcorrJZBeemm->Clone("Bpred");
  TH1F *JBpred = (TH1F*)JLcorrJZBeemm->Clone("Bpred");
    
  TH1F *BpredSys = new TH1F("Bpredsys","Bpredsys",PlottingSetup::global_ratio_binning.size()-1,&PlottingSetup::global_ratio_binning[0]);
  ClearHisto(BpredSys);
    
  if(PlottingSetup::RestrictToMassPeak) {
    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.);
    
    TTbarpred->Scale(1.0/3);
    Zjetspred->Add(LcorrJZBem,-1.0/3.);
    Zjetspred->Add(LcorrJZBSBem,-1.0/3.);
    TTbarpred->Add(RcorrJZBSBem,1.0/3.);
    Zjetspred->Add(LcorrJZBSBeemm,-1.0/3.);
    TTbarpred->Add(RcorrJZBSBeemm,1.0/3.);
    
    //these are for the ratio
    JBpred->Add(JRcorrJZBem,1.0/3.);
    JBpred->Add(JLcorrJZBem,-1.0/3.);
    JBpred->Add(JRcorrJZBSBem,1.0/3.);
    JBpred->Add(JLcorrJZBSBem,-1.0/3.);
    JBpred->Add(JRcorrJZBSBeemm,1.0/3.);
    JBpred->Add(JLcorrJZBSBeemm,-1.0/3.);
      
    //Systematics:
    AddSquared(BpredSys,JLcorrJZBeemm,zjetsestimateuncertONPEAK*zjetsestimateuncertONPEAK);
    AddSquared(BpredSys,JRcorrJZBem,emuncertONPEAK*emuncertONPEAK*(1.0/9));
    AddSquared(BpredSys,JLcorrJZBem,emuncertONPEAK*emuncertONPEAK*(1.0/9));
    AddSquared(BpredSys,JRcorrJZBSBem,emsidebanduncertONPEAK*emsidebanduncertONPEAK*(1.0/9));
    AddSquared(BpredSys,JLcorrJZBSBem,emsidebanduncertONPEAK*emsidebanduncertONPEAK*(1.0/9));
    AddSquared(BpredSys,JRcorrJZBSBeemm,eemmsidebanduncertONPEAK*eemmsidebanduncertONPEAK*(1.0/9));
    AddSquared(BpredSys,JLcorrJZBSBeemm,eemmsidebanduncertONPEAK*eemmsidebanduncertONPEAK*(1.0/9));
  } else {
    Bpred->Add(RcorrJZBem,1.0);
    Bpred->Add(LcorrJZBem,-1.0);
    
    Zjetspred->Add(LcorrJZBem,-1.0);
    
    //these are for the ratio
    JBpred->Add(JRcorrJZBem,1.0);
    JBpred->Add(JLcorrJZBem,-1.0);
      
    //Systematics
    AddSquared(BpredSys,JLcorrJZBeemm,zjetsestimateuncertOFFPEAK*zjetsestimateuncertOFFPEAK);
    AddSquared(BpredSys,JRcorrJZBem,emuncertOFFPEAK*emuncertOFFPEAK);
    AddSquared(BpredSys,JLcorrJZBem,emuncertOFFPEAK*emuncertOFFPEAK);

  }
    
  SQRT(BpredSys);
  BpredSys->Divide(JBpred);

  
  flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak ---- prediction changed
  TH1F *Tpred = (TH1F*)RcorrJZBem->Clone("Bpred");
  Tpred->Add(LcorrJZBem,-1.0);
  if(PlottingSetup::RestrictToMassPeak) {
    Tpred->Add(RcorrJZBSBem,1.0);
    Tpred->Add(LcorrJZBSBem,-1.0);
    Tpred->Add(RcorrJZBSBeemm,1.0);
    Tpred->Add(LcorrJZBSBeemm,-1.0);
  }
  
  globalcanvas->cd();
  globalcanvas->SetLogy(1);

  RcorrJZBeemm->SetMarkerStyle(20);
  RcorrJZBeemm->GetXaxis()->SetRangeUser(0,high);
  RcorrJZBeemm->SetMinimum(0.1);

  Bpred->SetLineColor(kRed);
  Bpred->SetStats(0);
  
  int versok=false;
  if(gROOT->GetVersionInt()>=53000) versok=true;


  if ( overlay_signal || use_data==2 ) lm4RcorrJZBeemm->SetLineColor(TColor::GetColor("#088A08"));
  RcorrJZBeemm->SetMarkerSize(DataMarkerSize);
  
  TLegend *legBpred = make_legend("",0.6,0.55);
  TLegend *legBpred2 = make_legend("",0.6,0.55);


  vector<TF1*> analytical_function;
  TPad *kinpad = new TPad("kinpad","kinpad",0,0,1,1);
  kinpad->cd();
  kinpad->SetLogy(1);

  string Bpredsaveas="Bpred_Data";
  blankback->SetMaximum(5*RcorrJZBeemm->GetMaximum());
  blankback->SetMinimum(0.1);
  if(use_data!=1) blankback->SetMinimum(0.1);
  blankback->Draw();
  if(use_data==1)
  {
    //Bpred->SetLineWidth(3); //paper style.overruled.
    //lm4RcorrJZBeemm->SetLineWidth(3); //paper style.overruled.
    analytical_function = do_extended_fit_to_plot(Bpred,Tpred,LcorrJZBeemm,LcorrJZBem,is_data);
    kinpad->cd();//necessary because the extended fit function creates its own canvas
    RcorrJZBeemm->Draw("e1x0,same");

    Bpred->Draw("hist,same");
    //analytical_function[0]->Draw("same"); analytical_function[1]->Draw("same");analytical_function[2]->Draw("same");
    RcorrJZBeemm->Draw("e1x0,same");//HAVE IT ON TOP!
    lm4RcorrJZBeemm->Draw("hist,same");
    legBpred->AddEntry(RcorrJZBeemm,"observed","p");
    legBpred->AddEntry(Bpred,"predicted","l");
//    legBpred->AddEntry(analytical_function[1],"predicted fit","l");
//    legBpred->AddEntry(analytical_function[2],"stat. uncert.","l");
    legBpred->AddEntry(lm4RcorrJZBeemm,(allsamples.collection[allsamples.FindSample("LM")[0]].samplename).c_str(),"l");
    legBpred->Draw();
    DrawPrelim();

    //this plot shows what the prediction is composed of
    TPad *predcomppad = new TPad("predcomppad","predcomppad",0,0,1,1);
    float CurrentBpredLineWidth=Bpred->GetLineWidth();
    Bpred->SetLineWidth(2);
    predcomppad->cd();
    predcomppad->SetLogy(1);
    
    TH1F *jzbnegative = (TH1F*)LcorrJZBeemm->Clone("jzbnegative");
    TH1F *sidebandsemu = (TH1F*)Bpred->Clone("sidebandsemu");
    sidebandsemu->Add(jzbnegative,-1);
    
    jzbnegative->SetFillColor(allsamples.GetColor((allsamples.FindSample("DY"))[0]));
    jzbnegative->SetLineColor(allsamples.GetColor((allsamples.FindSample("DY"))[0]));
    sidebandsemu->SetLineColor(allsamples.GetColor((allsamples.FindSample("TTJets"))[0]));
    sidebandsemu->SetFillColor(allsamples.GetColor((allsamples.FindSample("TTJets"))[0]));
    
    THStack predcomposition("predcomposition","prediction composition");
    predcomposition.Add(sidebandsemu);
    predcomposition.Add(jzbnegative);
    blankback->Draw();
    RcorrJZBeemm->Draw("e1x0,same");
    predcomposition.Draw("histo,same");//
    Bpred->Draw("hist,same");
//    analytical_function[0]->Draw("same"); analytical_function[1]->Draw("same");analytical_function[2]->Draw("same");
    RcorrJZBeemm->Draw("e1x0,same");//HAVE IT ON TOP!
//    lm4RcorrJZBeemm->SetLineColor(kOrange+1);
    lm4RcorrJZBeemm->SetLineWidth(2);
    //lm4RcorrJZBeemm->SetLineWidth(2); // paper style. overruled.
    lm4RcorrJZBeemm->Draw("histo,same");
    DrawPrelim();
    TLegend *speciallegBpred = make_legend("",0.45,0.55);
    //TLegend *speciallegBpred = make_legend("",0.35,0.55); // paper style. overruled.
    speciallegBpred->AddEntry(RcorrJZBeemm,"Data","pl");
    speciallegBpred->AddEntry(Bpred,"Total background","l");
    speciallegBpred->AddEntry(jzbnegative,"JZB<0 (data)","f");
    if(PlottingSetup::RestrictToMassPeak) speciallegBpred->AddEntry(sidebandsemu,"Sidebands/e#mu (data)","f");
    else speciallegBpred->AddEntry(sidebandsemu,"e#mu (data)","f");
//    speciallegBpred->AddEntry(lm4RcorrJZBeemmC,"LM4","l");
    speciallegBpred->AddEntry(lm4RcorrJZBeemm,"LM4","l");
    speciallegBpred->Draw();
    save_with_ratio(JRcorrJZBeemm,JBpred,predcomppad,subdir+"Bpred_Data_____PredictionComposition",true,true,"data/pred",BpredSys);
    
    TCanvas *specialcanv = new TCanvas("specialcanv","specialcanv");
    specialcanv->SetLogy(1);
//    THStack kostack = allsamples.DrawStack("RcorrJZBeemm",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,!is_data, luminosity,use_signal);
    blankback->Draw();
//    kostack.Draw("same");
    predcomposition.Draw();
    Bpred->Draw("hist,same");
    //analytical_function[0]->Draw("same"); analytical_function[1]->Draw("same");analytical_function[2]->Draw("same");
    RcorrJZBeemm->Draw("e1x0,same");//HAVE IT ON TOP!
    legBpred->Draw();
    DrawPrelim();
    CompleteSave(specialcanv,subdir+"Bpred_Data_____PredictionCompositioninMC");
    Bpred->SetLineWidth((int)CurrentBpredLineWidth);
    
    
    //for(int i=1;i<=Bpred->GetNbinsX();i++) cout << Bpred->GetBinLowEdge(i) << ";" << Bpred->GetBinLowEdge(i)+Bpred->GetBinWidth(i) << ";;" << RcorrJZBeemm->GetBinContent(i) << ";" << LcorrJZBeemm->GetBinContent(i) << ";" << RcorrJZBem->GetBinContent(i) << ";" << LcorrJZBem->GetBinContent(i) << endl;
    
    delete speciallegBpred;
    delete Zjetspred;
    delete TTbarpred;
    
    kinpad->cd();
  }
  if(use_data==0) {
    RcorrJZBeemm->Draw("e1x0,same");
    //Bpred->SetLineWidth(3); // paper style. overruled.
    Bpred->Draw("hist,same");
    RcorrJZBeemm->Draw("e1x0,same");//HAVE IT ON TOP!
    legBpred->AddEntry(RcorrJZBeemm,"MC true","p");
    legBpred->AddEntry(Bpred,"MC predicted","l");    
    if(versok) legBpred->AddEntry((TObject*)0,"",""); // Just for alignment // causes seg fault on root v5.18
    if(versok) legBpred->AddEntry((TObject*)0,"",""); // causes seg fault on root v5.18
    if ( overlay_signal ) legBpred->AddEntry(lm4RcorrJZBeemm,"LM4","l");
    legBpred->Draw();
    DrawMCPrelim();
    Bpredsaveas="Bpred_MC";
//    CompleteSave(globalcanvas,"Bpred_MC"); // done below in save_with_ratio
  }
  if(use_data==2) {
    RcorrJZBeemm->Draw("e1x0,same");
    //Bpred->SetLineWidth(3); // paper style. overruled.
    Bpred->Draw("hist,same");
    RcorrJZBeemm->Draw("e1x0,same");//HAVE IT ON TOP!
    legBpred->AddEntry(RcorrJZBeemm,"MC true","p");
    legBpred->AddEntry(Bpred,"MC predicted","l");
    legBpred2->AddEntry(RcorrJZBeemm,"MC true","p");
    legBpred2->AddEntry(Bpred,"MC predicted","l");
    {
      if(versok) legBpred->AddEntry((TObject*)0,"",""); // Just for alignment // causes seg fault on root v5.18 --> now only allowed for root >=v5.30
      if(versok) legBpred->AddEntry((TObject*)0,"",""); // causes seg fault on root v5.18 --> now only allowed for root >=v5.30
      legBpred->Draw();
      DrawMCPrelim();
      Bpredsaveas="Bpred_MCwithS";
//      CompleteSave(globalcanvas,"Bpred_MCwithS"); // done below in save_with_ratio
    }
    {
      //lm4RcorrJZBeemm->SetLineWidth(3); //paper style. overruled.
      //RcorrJZBeemmNoS->SetLineWidth(3); //paper style. overruled.
      //lm4RcorrJZBeemm->SetLineStyle(2); //paper style. overruled.
      //RcorrJZBeemmNoS->SetLineStyle(3); //paper style. overruled.
      //lm4RcorrJZBeemm->SetLineColor(kOrange+1); //paper style. overruled.
      
      RcorrJZBeemmNoS->SetLineStyle(2);
      legBpred2->AddEntry(RcorrJZBeemmNoS,"MC B","l");
      legBpred2->AddEntry(lm4RcorrJZBeemm,"MC S","l");
      legBpred2->Draw();
      RcorrJZBeemmNoS->SetLineColor(TColor::GetColor("#61210B"));
      RcorrJZBeemmNoS->Draw("histo,same");
      RcorrJZBeemm->Draw("e1x0,same");//HAVE IT ON TOP!
      lm4RcorrJZBeemm->Draw("histo,same");
      DrawMCPrelim();
      Bpredsaveas="Bpred_MCwithS__plus";
//      CompleteSave(globalcanvas,"Bpred_MCwithS__plus"); // done below in save_with_ratio
    }
  }


  //3rd last argument: do special bpred ratio, 2nd last argument: extended range!, last: y-axis title
  string ytitle("ratio");
  if ( use_data==1 ) ytitle = "data/pred";
  //save_with_ratio(JRcorrJZBeemm,JBpred,kinpad,Bpredsaveas,true,use_data!=1,ytitle);
  save_with_ratio(JRcorrJZBeemm,JBpred,kinpad,subdir+Bpredsaveas,true,false,ytitle,BpredSys);//not extending the y range anymore up to 4

  
  flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
  // The part below is meaningless for the offpeak analysis (it's a comparison of the different estimates but there is but one estimate!)
  if(PlottingSetup::RestrictToMassPeak) {
    TH1F *Bpredem = (TH1F*)LcorrJZBeemm->Clone("Bpredem");
    Bpredem->Add(RcorrJZBem);
    Bpredem->Add(LcorrJZBem,-1);
    TH1F *BpredSBem = (TH1F*)LcorrJZBeemm->Clone("BpredSBem");
    BpredSBem->Add(RcorrJZBSBem);
    Bpred->Add(LcorrJZBSBem,-1);
    TH1F *BpredSBeemm = (TH1F*)LcorrJZBeemm->Clone("BpredSBeemm");  
    BpredSBeemm->Add(RcorrJZBSBeemm);
    BpredSBeemm->Add(LcorrJZBSBeemm,-1.0);
    globalcanvas->cd();
    globalcanvas->SetLogy(1);
    
    RcorrJZBeemm->SetMarkerStyle(20);
    RcorrJZBeemm->GetXaxis()->SetRangeUser(0,high);
    blankback->Draw();
    RcorrJZBeemm->Draw("e1x0,same");
    RcorrJZBeemm->SetMarkerSize(DataMarkerSize);
    
    Bpredem->SetLineColor(kRed+1);
    Bpredem->SetStats(0);
    Bpredem->Draw("hist,same");
    
    BpredSBem->SetLineColor(kGreen+2);//TColor::GetColor("#0B6138"));
    BpredSBem->SetLineStyle(2);
    BpredSBem->Draw("hist,same");
    
    BpredSBeemm->SetLineColor(kBlue+1);
    BpredSBeemm->SetLineStyle(3);
    BpredSBeemm->Draw("hist,same");
    RcorrJZBeemm->Draw("e1x0,same");//HAVE IT ON TOP!
    
    TLegend *legBpredc = make_legend("",0.6,0.55);
    if(use_data==1)
    {
      legBpredc->AddEntry(RcorrJZBeemm,"observed","p");
      legBpredc->AddEntry(Bpredem,"OFZP","l");
      legBpredc->AddEntry(BpredSBem,"OFSB","l");
      legBpredc->AddEntry(BpredSBeemm,"SFSB","l");
      legBpredc->Draw();
      CompleteSave(globalcanvas,subdir+"Bpred_Data_comparison");
    }
    if(use_data==0) {
      legBpredc->AddEntry(RcorrJZBeemm,"MC true","p");
      legBpredc->AddEntry(Bpredem,"MC OFZP","l");
      legBpredc->AddEntry(BpredSBem,"MC OFSB","l");
      legBpredc->AddEntry(BpredSBeemm,"MC SFSB","l");
      legBpredc->Draw();
      CompleteSave(globalcanvas,subdir+"Bpred_MC_comparison");
    }
    if(use_data==2) {
      legBpredc->AddEntry(RcorrJZBeemm,"MC true","p");
      legBpredc->AddEntry(Bpredem,"MC OFZP","l");
      legBpredc->AddEntry(BpredSBem,"MC OFSB","l");
      legBpredc->AddEntry(BpredSBeemm,"MC SFSB","l");
      if ( overlay_signal ) legBpred->AddEntry(lm4RcorrJZBeemm,"LM4","l");
      legBpredc->Draw();
      CompleteSave(globalcanvas,subdir+"Bpred_MCwithS_comparison");
    }
  }
  
  TFile *f = new TFile("tester.root","RECREATE");
  RcorrJZBeemm->Write();
  Bpred->Write();
  f->Close();
  
  delete RcorrJZBeemm;
  delete LcorrJZBeemm;
  delete RcorrJZBem;
  delete LcorrJZBem;
  
  delete JRcorrJZBeemm;
  delete JLcorrJZBeemm;
  delete JRcorrJZBem;
  delete JLcorrJZBem;

  delete blankback;
  
  delete BpredSys;
  if(PlottingSetup::RestrictToMassPeak) {
    delete RcorrJZBSBem;
    delete LcorrJZBSBem;
    delete RcorrJZBSBeemm;
    delete LcorrJZBSBeemm;

    delete JRcorrJZBSBem;
    delete JLcorrJZBSBem;
    delete JRcorrJZBSBeemm;
    delete JLcorrJZBSBeemm;
  }
  if(overlay_signal || use_data==2) delete lm4RcorrJZBeemm;
  switch_overunderflow(false);
}

void do_prediction_plots(string mcjzb, string datajzb, float DataSigma, float MCSigma, bool overlay_signal ) {
  TCanvas *globalcanvas = new TCanvas("globalcanvas","Prediction Canvas");
  do_prediction_plot(datajzb,globalcanvas,jzbHigh ,data,overlay_signal);
  if ( !PlottingSetup::Approved ) {
    do_prediction_plot(mcjzb,globalcanvas,jzbHigh ,mc,overlay_signal);
    do_prediction_plot(mcjzb,globalcanvas,jzbHigh ,mcwithsignal,overlay_signal);
  } else {
    write_info(__FUNCTION__,"You set approved to true, therefore not producing prediction/observation plots for MC with and without signal.");
  }
}

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;
  TH1F *LcorrJZBSBem;
  TH1F *RcorrJZBSBeemm;
  TH1F *LcorrJZBSbeemm;
  
  if(PlottingSetup::RestrictToMassPeak) {
    RcorrJZBSBem   = allsamples.Draw("RcorrJZBSbem",jzb.c_str(),binning, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets,do_data, luminosity,dosignal);
    LcorrJZBSBem   = allsamples.Draw("LcorrJZBSbem",("-"+jzb).c_str(),binning, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets,do_data, luminosity,dosignal);
    RcorrJZBSBeemm = allsamples.Draw("RcorrJZBSbeemm",jzb.c_str(),binning, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets,do_data, luminosity,dosignal);
    LcorrJZBSbeemm = allsamples.Draw("LcorrJZBSbeemm",("-"+jzb).c_str(),binning, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets,do_data, luminosity,dosignal);
  }
  

  flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
  TH1F *Bpred = (TH1F*)LcorrJZBeemm->Clone("Bpred");
  if(PlottingSetup::RestrictToMassPeak) {
    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);
  } else {
    Bpred->Add(RcorrJZBem,1.0);
    Bpred->Add(LcorrJZBem,-1.0);
  }

  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,false);
  
  
  JZBratio->SetTitle("");
  JZBratio->GetYaxis()->SetRangeUser(0.0,9.0);
//  if(is_data==1) JZBratio->GetXaxis()->SetRangeUser(0,jzbHigh );
  
  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]");
  if ( high>0 ) JZBratio->GetXaxis()->SetRangeUser(0.0,high);
  JZBratio->GetYaxis()->SetNdivisions(519);
  JZBratio->GetYaxis()->SetRangeUser(0.0,4.0);
  JZBratio->GetYaxis()->CenterTitle();
  JZBratio->GetXaxis()->CenterTitle();
  JZBratio->SetMarkerSize(DataMarkerSize);
  JZBratio->Draw("AP");
  /////----------------------------
  TPaveText *writeresult = new TPaveText(0.15,0.78,0.49,0.91,"blNDC");
  writeresult->SetFillStyle(4000);
  writeresult->SetFillColor(kWhite);
  writeresult->SetTextFont(42);
  writeresult->SetTextSize(0.03);
  writeresult->SetTextAlign(12);
  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.5",0,1000);
  TF1 *bottomline =  new TF1("","0.5",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",0.6,0.55,false);//this line is just to have the default CMS Preliminary (...) on the canvas as well.
  if(is_data==1) DrawPrelim();
  else DrawMCPrelim();
  TLegend *leg = new TLegend(0.55,0.75,0.89,0.89);
  leg->SetTextFont(42);
  leg->SetTextSize(0.04);
//   if(is_data==1) leg->SetHeader("Ratio (data)");
//   else leg->SetHeader("Ratio (MC)");

  TString MCtitle("MC ");
  if (is_data==1) MCtitle = "";

  leg->SetFillStyle(4000);
  leg->SetFillColor(kWhite);
  leg->SetTextFont(42);
//  leg->AddEntry(topline,"+20\% sys envelope","l");
  leg->AddEntry(JZBratio,MCtitle+"obs / "+MCtitle+"pred","p");
  leg->AddEntry(oneline,"ratio = 1","l");
//  leg->AddEntry(pol0d,"fit in [0,30] GeV","l");
//  leg->AddEntry(bottomline,"#pm50% envelope","l");


  //leg->Draw("same"); // no longer drawing legend

  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,vector<float> ratio_binning) {
  TCanvas *globalc = new TCanvas("globalc","Ratio Plot Canvas");
  globalc->SetLogy(0);
  
  do_ratio_plot(mc,ratio_binning,mcjzb,globalc, jzbHigh  );
  do_ratio_plot(data,ratio_binning,datajzb,globalc, jzbHigh  );
  do_ratio_plot(mcwithsignal,ratio_binning,mcjzb,globalc, jzbHigh  );
}

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();
}


void lepton_comparison_plots() {
  Float_t ymin = 1.e-5, ymax = 0.25;
  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("/DY"));
  TH1F *mmmc = allsamples.Draw("mmmc","mll",50,50,150, "mll [GeV]", "events", cutOSSF&&cutnJets&&"(id1==1)",mc, luminosity,allsamples.FindSample("/DY"));
  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, "lepton_comparison/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, "lepton_comparison/mll_effratio_data");
  
  TH1F *jeed = allsamples.Draw("jeed",jzbvariabledata,    int((jzbHigh+150)/5),-150,jzbHigh , "JZB [GeV]", "events", cutOSSF&&cutnJets&&"(id1==0)",data, luminosity);
  TH1F *jmmd = allsamples.Draw("jmmd",jzbvariabledata,    int((jzbHigh+150)/5),-150,jzbHigh , "JZB [GeV]", "events", cutOSSF&&cutnJets&&"(id1==1)",data, luminosity);
  TH1F *jeemmd = allsamples.Draw("jeemmd",jzbvariabledata,int((jzbHigh+150)/5),-150,jzbHigh , "JZB [GeV]", "events", cutOSSF&&cutnJets,data, luminosity);
  dout << "ee : " << jeed->GetMean() << "+/-" << jeed->GetMeanError() << endl;
  dout << "ee : " << jmmd->GetMean() << "+/-" << jmmd->GetMeanError() << endl;
  dout << "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());
  TH1* njeemmd = jeemmd->DrawNormalized();
  njeemmd->SetMinimum(ymin);
  njeemmd->SetMaximum(ymax);

  jeed->DrawNormalized("histo,same");
  jmmd->DrawNormalized("histo,same");
  jeemmd->DrawNormalized("same");
  TLegend *jleg2 = make_legend(" ");
  jleg2->AddEntry(jeemmd,"ee and #mu#mu","p");
  jleg2->AddEntry(jeed,"ee","l");
  jleg2->AddEntry(jmmd,"#mu#mu","l");
  jleg2->Draw();
  CompleteSave(can,"lepton_comparison/jzb_effratio_data");

  TPad *eemmpad = new TPad("eemmpad","eemmpad",0,0,1,1);
  eemmpad->cd();
  eemmpad->SetLogy(1);
  jeed->Draw("histo");
  jmmd->Draw("histo,same");
  TLegend *eemmlegend = make_legend(" ");
  eemmlegend->AddEntry(jeed,"ee","l");
  eemmlegend->AddEntry(jmmd,"#mu#mu","l");
  eemmlegend->Draw();
  DrawPrelim();
  save_with_ratio(jeed,jmmd,eemmpad->cd(),"lepton_comparison/jzb_Comparing_ee_mm_data");
  
  TH1F *zjeed = allsamples.Draw("zjeed",jzbvariablemc,    int((jzbHigh+150)/5),-150,jzbHigh , "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets&&"(id1==0)",mc,  luminosity,allsamples.FindSample("/DY"));
  TH1F *zjmmd = allsamples.Draw("zjmmd",jzbvariablemc,    int((jzbHigh+150)/5),-150,jzbHigh , "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets&&"(id1==1)",mc,  luminosity,allsamples.FindSample("/DY"));
  TH1F *zjeemmd = allsamples.Draw("zjeemmd",jzbvariablemc,int((jzbHigh+150)/5),-150,jzbHigh , "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,            mc,  luminosity,allsamples.FindSample("/DY"));
  dout << "Z+Jets ee : " << zjeed->GetMean() << "+/-" << zjeed->GetMeanError() << endl;
  dout << "Z+Jets ee : " << zjmmd->GetMean() << "+/-" << zjmmd->GetMeanError() << endl;
  dout << "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());

  TH1* nzjeemmd = zjeemmd->DrawNormalized();
  nzjeemmd->SetMinimum(ymin);
  nzjeemmd->SetMaximum(ymax);
  zjeed->DrawNormalized("histo,same");
  zjmmd->DrawNormalized("histo,same");
  zjeemmd->DrawNormalized("same");
  TLegend *zjleg2 = make_legend("Z+jets MC");
  zjleg2->AddEntry(jeemmd,"ee and #mu#mu","p");
  zjleg2->AddEntry(jeed,"ee","l");
  zjleg2->AddEntry(jmmd,"#mu#mu","l");
  zjleg2->Draw();
  CompleteSave(can,"lepton_comparison/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,"lepton_comparison/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,int((jzbHigh+110)/5),-110,jzbHigh , "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data, luminosity);
  TH1F *jzbemd = allsamples.Draw("jzbemd",jzbvariabledata,int((jzbHigh+110)/5),-110,jzbHigh , "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();
  leg4->AddEntry(jzbld,"SFZP","p");
  leg4->AddEntry(jzbemd,"OFZP","p");
  leg4->AddEntry((TObject*)0,"",""); //causes segmentation violation
  leg4->AddEntry((TObject*)0,"",""); //causes segmentation violation
  leg4->Draw();
  CompleteSave(can,"lepton_comparison/jzb_eemumu_emu_data");
  
  TH1F *ttbarjzbld  = allsamples.Draw("ttbarjzbld",jzbvariablemc,int((jzbHigh+150)/5),-150,jzbHigh, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data, luminosity,allsamples.FindSample("TTJet"));
  TH1F *ttbarjzbemd = allsamples.Draw("ttbarjzbemd",jzbvariablemc,int((jzbHigh+150)/5),-150,jzbHigh, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data, luminosity,allsamples.FindSample("TTJet"));
  ttbarjzbld->SetLineColor(allsamples.GetColor("TTJet"));
  ttbarjzbemd->SetLineColor(allsamples.GetColor("TTJet"));
  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,"lepton_comparison/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;
}

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


void draw_pure_jzb_histo(TCut cut,string datavariable, string mcvariable, string savename, TCanvas *can,vector<float> binning) {
  TPad *jzbpad = new TPad("jzbpad","jzbpad",0,0,1,1);
  jzbpad->cd();
  jzbpad->SetLogy(1);
  string xlabel="JZB [GeV]";

  TH1F *datahisto = allsamples.Draw("datahisto",datavariable,binning, xlabel, "events",cut,data,luminosity);
  THStack mcstack = allsamples.DrawStack("mcstack",mcvariable,binning, xlabel, "events",cut,mc,luminosity);
  
  datahisto->SetMinimum(0.1);
  datahisto->SetMarkerSize(DataMarkerSize);
  datahisto->Draw("e1");
  mcstack.Draw("same");
  datahisto->Draw("same,e1");
  
  TLegend *leg;
  if (!PlottingSetup::RestrictToMassPeak) {
    if(is_OF(cut)) leg = allsamples.allbglegend("Opposite flavor",datahisto);
    else leg = allsamples.allbglegend("Same flavor",datahisto);
  } else { 
    if(is_OF(cut) && is_ZP(cut)) leg = allsamples.allbglegend("OFZP",datahisto);
    else if(!is_OF(cut) &&  is_ZP(cut)) leg = allsamples.allbglegend("SFZP",datahisto);
    else if( is_OF(cut) && !is_ZP(cut)) leg = allsamples.allbglegend("OFSB",datahisto);
    else if(!is_OF(cut) && !is_ZP(cut)) leg = allsamples.allbglegend("SFSB",datahisto);
    else {
      std::cerr << "Unable to decode cut: " << cut.GetTitle() << std::endl;
      exit(-1);
    }
  }
  leg->Draw();
  string write_cut = decipher_cut(cut,"");
  TText *writeline1 = write_cut_on_canvas(write_cut.c_str());
  writeline1->SetTextSize(0.035);
  writeline1->Draw();
  if(!Contains(savename,"Dibosons")) save_with_ratio(datahisto,mcstack,jzbpad->cd(),("jzb/"+savename));
  else save_with_ratio(datahisto,mcstack,jzbpad->cd(),savename);
  TPad *jzbpad2 = new TPad("jzbpad2","jzbpad2",0,0,1,1);
  jzbpad2->cd();
  jzbpad2->SetLogy(1);
  datahisto->GetXaxis()->SetRangeUser(0,binning[binning.size()-1]);
  datahisto->SetMinimum(0.1);
  datahisto->SetMarkerSize(DataMarkerSize);
  datahisto->Draw("e1");
  mcstack.Draw("same");
  datahisto->Draw("same,e1");
  leg->SetHeader("");
  leg->Draw();
  writeline1->SetTextSize(0.035);
  writeline1->Draw();
  DrawPrelim();
  if(!Contains(savename,"Dibosons")) save_with_ratio(datahisto,mcstack,jzbpad2->cd(),("jzb/PositiveSideOnly/"+savename+""));
  else save_with_ratio(datahisto,mcstack,jzbpad2->cd(),(savename+"__PosOnly"));
  datahisto->Delete();
  mcstack.Delete();
}

Double_t GausR(Double_t *x, Double_t *par) {
  return gRandom->Gaus(x[0],par[0]);
}

void produce_stretched_jzb_plots(string mcjzb, string datajzb,vector<float> ratio_binning) {
  TCanvas *dican = new TCanvas("dican","JZB Plots Canvas");
  float max=jzbHigh ;
  float min=-120;
  int nbins=(int)((max-min)/5.0); // we want 5 GeV/bin
  int coarserbins=int(nbins/2.0);
  int rebinnedbins=int(nbins/4.0);
  
  vector<float>binning;vector<float>coarse_binning;vector<float>coarsest_binning;
  for(int i=0;i<=nbins;i++)binning.push_back(min+i*(max-min)/((float)nbins));
  for(int i=0;i<=coarserbins;i++)coarse_binning.push_back(min+i*(max-min)/((float)coarserbins));
  for(int i=0;i<=rebinnedbins;i++)coarsest_binning.push_back(min+i*(max-min)/((float)rebinnedbins));

  draw_pure_jzb_histo(cutOSSF&&cutnJets&&cutmass,datajzb,mcjzb,"Dibosons/jzb_OS_SFZP",dican,binning);
  draw_pure_jzb_histo(cutOSOF&&cutnJets&&cutmass,datajzb,mcjzb,"Dibosons/jzb_OS_OFZP",dican,binning);
  draw_pure_jzb_histo(cutOSSF&&cutnJets&&cutmass&&"id1==0",datajzb,mcjzb,"Dibosons/ee/jzb_OS_SFZP",dican,binning);
  draw_pure_jzb_histo(cutOSSF&&cutnJets&&cutmass&&"id1==1",datajzb,mcjzb,"Dibosons/mm/jzb_OS_SFZP",dican,binning);
  draw_pure_jzb_histo(cutOSOF&&cutnJets&&cutmass&&"id1==0",datajzb,mcjzb,"Dibosons/ee/jzb_OS_OFZP",dican,binning);
  draw_pure_jzb_histo(cutOSOF&&cutnJets&&cutmass&&"id1==1",datajzb,mcjzb,"Dibosons/mm/jzb_OS_OFZP",dican,binning);
  if(PlottingSetup::RestrictToMassPeak) draw_pure_jzb_histo(cutOSSF&&cutnJets&&sidebandcut,datajzb,mcjzb,"Dibosons/jzb_OS_SFSB",dican,binning);
  if(PlottingSetup::RestrictToMassPeak) draw_pure_jzb_histo(cutOSOF&&cutnJets&&sidebandcut,datajzb,mcjzb,"Dibosons/jzb_OS_OFSB",dican,binning);
  
  draw_pure_jzb_histo(cutOSSF&&cutnJets&&cutmass,datajzb,mcjzb,"Dibosons/jzb_OS_SFZP_coarse",dican,coarse_binning);
  draw_pure_jzb_histo(cutOSOF&&cutnJets&&cutmass,datajzb,mcjzb,"Dibosons/jzb_OS_OFZP_coarse",dican,coarse_binning);
  if(PlottingSetup::RestrictToMassPeak) draw_pure_jzb_histo(cutOSSF&&cutnJets&&sidebandcut,datajzb,mcjzb,"Dibosons/jzb_OS_SFSB_coarse",dican,coarse_binning);
  if(PlottingSetup::RestrictToMassPeak) draw_pure_jzb_histo(cutOSOF&&cutnJets&&sidebandcut,datajzb,mcjzb,"Dibosons/jzb_OS_OFSB_coarse",dican,coarse_binning);
  
  delete dican;
}
  

void diboson_plots(string mcjzb, string datajzb,vector<float> ratio_binning) {
  vector<int> SamplesToBeModified = allsamples.FindSampleBySampleName("WW/WZ/ZZ");
  
  if(SamplesToBeModified.size()==0 || SamplesToBeModified[0]==-1) {
    write_error(__FUNCTION__,"Could not find any diboson samples - aborting diboson plots");
    return;
  }
  
  float stretchfactor = 100.0;
  vector<string> labels;
  
  
  dout << "Going to increase the cross section for diboson samples ... " << endl;
  for(int i=0;i<(int)SamplesToBeModified.size();i++) {
    float origxs=(allsamples.collection)[SamplesToBeModified[i]].xs;
    (allsamples.collection)[SamplesToBeModified[i]].xs=origxs*stretchfactor;
    dout << "     Increased xs for sample " << (allsamples.collection)[SamplesToBeModified[i]].filename << " from " << origxs << " to " << (allsamples.collection)[SamplesToBeModified[i]].xs << " (by a factor of " << stretchfactor << ")" << endl;
    labels.push_back((allsamples.collection)[SamplesToBeModified[i]].samplename);
    (allsamples.collection)[SamplesToBeModified[i]].samplename=any2string(int(stretchfactor))+" x "+(allsamples.collection)[SamplesToBeModified[i]].samplename;
    dout << "         (also renamed it to " << (allsamples.collection)[SamplesToBeModified[i]].samplename << " )" << endl;
  }
  
  dout << "Going to produce JZB plots" << endl;
  produce_stretched_jzb_plots(mcjzb,datajzb,ratio_binning);
  TCanvas *gloca = new TCanvas("gloca","gloca");
  
  dout << "Going to produce prediction plots" << endl;
  do_prediction_plot(mcjzb, gloca, PlottingSetup::jzbHigh, 0, false,"Dibosons/Bpred/" ); // do only MC plots, no signal
  do_prediction_plot(mcjzb, gloca, PlottingSetup::jzbHigh, 0, false,"Dibosons/Bpred/" ); // do MC plots with signal
  delete gloca;

  dout << "Going to reset the cross section for diboson samples ... " << endl;
  for(int i=0;i<(int)SamplesToBeModified.size();i++) {
    float Upxs=(allsamples.collection)[SamplesToBeModified[i]].xs;
    (allsamples.collection)[SamplesToBeModified[i]].xs=(allsamples.collection)[SamplesToBeModified[i]].xs*(1.0/stretchfactor);
    string Upname=(allsamples.collection)[SamplesToBeModified[i]].samplename;
    (allsamples.collection)[SamplesToBeModified[i]].samplename=labels[i];
    dout << "     Reset xs for sample " << (allsamples.collection)[SamplesToBeModified[i]].samplename << " from " << Upxs << " to " << (allsamples.collection)[SamplesToBeModified[i]].xs << " (by a factor of " << stretchfactor << ") and reset the correct name (from " << Upname << ")" << endl;
    
  }
  
}


void draw_normalized_data_vs_data_histo(TCut cut1, TCut cut2, string variable, string legentry1, string legentry2, string savename, TCanvas *can,vector<float> binning) {
  TPad *jzbpad = new TPad("jzbpad","jzbpad",0,0,1,1);
  jzbpad->cd();
  jzbpad->SetLogy(1);
  string xlabel="JZB [GeV]";

  TH1F *datahisto1 = allsamples.Draw("datahisto1",variable,binning, xlabel, "events",cut1,data,luminosity);
  datahisto1->SetLineColor(kRed);
  datahisto1->SetMarkerColor(kRed);
  TH1F *datahisto2 = allsamples.Draw("datahisto2",variable,binning, xlabel, "events",cut2,data,luminosity);
  datahisto2->SetLineColor(kBlue);
  datahisto2->SetMarkerColor(kBlue);
  
  datahisto2->SetMarkerSize(DataMarkerSize);
  datahisto1->DrawNormalized("e1");
  datahisto2->DrawNormalized("histo,same");
  datahisto1->DrawNormalized("same,e1");
  
  TLegend *leg = make_legend();
  leg->AddEntry(datahisto1,legentry1.c_str());
  leg->AddEntry(datahisto2,legentry2.c_str());
  leg->Draw();
  
  save_with_ratio(datahisto1,datahisto2,jzbpad->cd(),("jzb/"+savename));
  
  datahisto1->Delete();
  datahisto2->Delete();
}


void jzb_plots(string mcjzb, string datajzb,vector<float> ratio_binning) {
  TCanvas *can = new TCanvas("can","JZB Plots Canvas");
  float max=jzbHigh ;
  float min=-120;
  int nbins=(int)((max-min)/5.0); // we want 5 GeV/bin
  int coarserbins=int(nbins/2.0);
  int rebinnedbins=int(nbins/4.0);
  
  vector<float>binning;vector<float>coarse_binning;vector<float>coarsest_binning;
  for(int i=0;i<=nbins;i++)binning.push_back(min+i*(max-min)/((float)nbins));
  for(int i=0;i<=coarserbins;i++)coarse_binning.push_back(min+i*(max-min)/((float)coarserbins));
  for(int i=0;i<=rebinnedbins;i++)coarsest_binning.push_back(min+i*(max-min)/((float)rebinnedbins));

  if ( !PlottingSetup::Approved ) {
    draw_pure_jzb_histo(cutOSSF&&cutnJets&&cutmass,datajzb,mcjzb,"jzb_OS_SFZP",can,binning);
    draw_pure_jzb_histo(cutOSOF&&cutnJets&&cutmass,datajzb,mcjzb,"jzb_OS_OFZP",can,binning);
    draw_pure_jzb_histo(cutOSSF&&cutnJets&&cutmass&&"id1==0",datajzb,mcjzb,"ee/jzb_OS_SFZP",can,binning);
    draw_pure_jzb_histo(cutOSSF&&cutnJets&&cutmass&&"id1==1",datajzb,mcjzb,"mm/jzb_OS_SFZP",can,binning);
    draw_pure_jzb_histo(cutOSOF&&cutnJets&&cutmass&&"id1==0",datajzb,mcjzb,"ee/jzb_OS_OFZP",can,binning);
    draw_pure_jzb_histo(cutOSOF&&cutnJets&&cutmass&&"id1==1",datajzb,mcjzb,"mm/jzb_OS_OFZP",can,binning);
    flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
    if(PlottingSetup::RestrictToMassPeak) draw_pure_jzb_histo(cutOSSF&&cutnJets&&sidebandcut,datajzb,mcjzb,"jzb_OS_SFSB",can,binning);
    if(PlottingSetup::RestrictToMassPeak) draw_pure_jzb_histo(cutOSOF&&cutnJets&&sidebandcut,datajzb,mcjzb,"jzb_OS_OFSB",can,binning);
    draw_normalized_data_vs_data_histo(cutOSOF&&cutnJets&&cutmass&&"id1==0",cutOSOF&&cutnJets&&cutmass&&"id1==1",datajzb,"ee","mm","jzb_ee_vs_mm",can,binning);
    draw_normalized_data_vs_data_histo(cutOSOF&&cutnJets&&cutmass&&"id1==0",cutOSOF&&cutnJets&&cutmass&&"id1==1",datajzb,"ee","mm","jzb_ee_vs_mm_coarse",can,coarse_binning);
    draw_normalized_data_vs_data_histo(cutOSOF&&cutnJets&&cutmass&&"id1==0",cutOSOF&&cutnJets&&cutmass&&"id1==1",datajzb,"ee","mm","jzb_ee_vs_mm_coarsest",can,coarsest_binning);

  }
  
  draw_pure_jzb_histo(cutOSSF&&cutnJets&&cutmass,datajzb,mcjzb,"jzb_OS_SFZP_coarse",can,coarse_binning);
  if ( !PlottingSetup::Approved ) {
    draw_pure_jzb_histo(cutOSOF&&cutnJets&&cutmass,datajzb,mcjzb,"jzb_OS_OFZP_coarse",can,coarse_binning);
    flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
    if(PlottingSetup::RestrictToMassPeak) draw_pure_jzb_histo(cutOSSF&&cutnJets&&sidebandcut,datajzb,mcjzb,"jzb_OS_SFSB_coarse",can,coarse_binning);
    if(PlottingSetup::RestrictToMassPeak) draw_pure_jzb_histo(cutOSOF&&cutnJets&&sidebandcut,datajzb,mcjzb,"jzb_OS_OFSB_coarse",can,coarse_binning);
  }
  delete can;
}


void calculate_all_yields(string mcdrawcommand,vector<float> jzb_cuts) {
  dout << "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);
    dout << "Above " << allbgs->GetBinLowEdge(i) << " GeV/c : " << cumulative << endl;
  }
}


TCut give_jzb_expression(string mcjzb, float jzbcut, string posneg="pos") {
  stringstream res;
  res << "(" << mcjzb;
  if(posneg=="pos") res << ">";
  else res << "<-";
  res << jzbcut << ")";
  return TCut(res.str().c_str());
}

string sigdig(float number, int nsigdig=3, float min=0) {
  //this function tries to extract n significant digits, and if the number is below (min), it returns "<min"
  if(number<min) return "< "+any2string(min);
  stringstream sol;
  sol << setprecision(nsigdig) << number;
  
  return sol.str();
}
  
string jzb_tex_command(string region, string posneg) {
  if(posneg=="pos") posneg="POS";
  else posneg="NEG";
  stringstream texcommand;
  texcommand<<"\\"<<region <<"JZB"<<posneg;
  return texcommand.str();
}
//  \SFZPJZBPOS
//  Sample &  \OFZPJZBPOS  & \OFZPJZBNEG &  \SFZPJZBPOS & \SFZPJZBNEG & \OFSBJZBPOS  & \OFSBJZBNEG  & \SFSBJZBPOS  & \SFSBJZBNEG  \\\hline
  
void compute_MC_yields(string mcjzb,vector<float> jzb_cuts) {
  dout << "Calculating background yields in MC:" << endl;
  
  TCanvas *yica = new TCanvas("yica","yield canvas");
  
  int nRegions=4;
  if(!PlottingSetup::RestrictToMassPeak) nRegions=2;
  string tsRegions[] = {"SFZP","OFZP","SFSB","OFSB"};
  string posneg[] = {"pos","neg"};
  TCut tkRegions[] = {cutOSSF&&cutnJets&&cutmass,cutOSOF&&cutnJets&&cutmass,cutOSSF&&cutnJets&&sidebandcut,cutOSOF&&cutnJets&&sidebandcut};

  for(int ijzb=0;ijzb<(int)jzb_cuts.size();ijzb++) {
    TCut jzbMC[]  = { give_jzb_expression(mcjzb,jzb_cuts[ijzb],"pos"), give_jzb_expression(mcjzb,jzb_cuts[ijzb],"neg") };
    dout << "_________________________________________________________" << endl;
    dout << "Table for JZB> " << jzb_cuts[ijzb] << endl;
    for(int isample=0;isample<(int)(allsamples.collection).size();isample++) {
      if(!(allsamples.collection)[isample].is_data) dout << (allsamples.collection)[isample].samplename << "  &  ";
      else dout << "Sample & ";
      for(int iregion=0;iregion<nRegions;iregion++) {
        for(int ipos=0;ipos<2;ipos++) {
          if((allsamples.collection)[isample].is_data) dout << jzb_tex_command(tsRegions[iregion],posneg[ipos]) << "  &  ";
          else {
            vector<int> specific;specific.push_back(isample);
            TH1F *shisto = allsamples.Draw("shisto","mll",1,0,500,"tester","events",tkRegions[iregion]&&jzbMC[ipos],mc,luminosity,specific);
            dout << sigdig(shisto->Integral(),3,0.05) <<"  &  ";
            delete shisto;
          }
      }//end of ipos
    }//end of iregion
    dout << " \\\\" << endl;
  }//end of isample
  }//end of ijzb
  dout << "  \\hline" << endl;
  
  delete yica;
}

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,jzbHigh, "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) {
    dout << "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);
    dout << "--> Quality of Z+Jets / TTbar fit : chi2/ndf = " << ZplusJetsplusTTbar->GetChisquare() << "/" << ZplusJetsplusTTbar->GetNDF() << endl;
    ZplusJetsplusTTbar->Draw("same");
    ZplusJetsplusTTbarP->SetParameters(ZplusJetsplusTTbar->GetParameters());
    ZplusJetsplusTTbarN->SetParameters(ZplusJetsplusTTbar->GetParameters());
    dout << "TTbar contribution found (after fitting) : " << ZplusJetsplusTTbar->GetParameter(5) << endl;
    float factor = ZplusJetsplusTTbar->GetParameter(5) / logparc->GetParameter(0);
    dout << "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);
  }
  dout << "Statistics about our fits: " << endl; 
  dout << "Z+Jets shape: Chi2/ndf = " << zjetsfunc->GetChisquare() << "/" << ossfleft->GetNbinsX() << endl;
  dout << "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 threejetswith30gev=0;
/*  
  int twojetswith50gev=1;
  int electrons_only=0;
  int mu_only=1;
 
  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);
}

float find_one_correction_factor(string FindKeyword, bool dodata, TCut SpecialCut, string SaveAs) {
  TCanvas *cancorr = new TCanvas("cancorr","Canvas for Response Correction");
  cancorr->SetLogz();
  cancorr->SetRightMargin(0.13);
  TCut zptforresponsepresentation(Restrmasscut&&SpecialCut&&passtrig);

  if(PlottingSetup::DoBTag) zptforresponsepresentation=zptforresponsepresentation&&PlottingSetup::bTagRequirement;
  TH2F *niceresponseplotd = new TH2F("niceresponseplotd","",100,0,600,100,0,5);
  niceresponseplotd->Sumw2();
  TH2F* emuResponse = (TH2F*)niceresponseplotd->Clone("emuResponse");
  vector<int> SampleIndices=allsamples.FindSample(FindKeyword);
  for(int iSample=0;iSample<(int)SampleIndices.size();iSample++) {
    if((allsamples.collection)[SampleIndices[iSample]].is_data && !dodata) continue;
    if((allsamples.collection)[SampleIndices[iSample]].is_data ==false && dodata) continue;
    
    dout << "   Response correction : Using sample " << (allsamples.collection)[SampleIndices[iSample]].filename << " for " << FindKeyword << endl;
    (allsamples.collection)[SampleIndices[iSample]].events->Draw("sumJetPt[1]/pt:pt>>+niceresponseplotd",(zptforresponsepresentation&&cutOSSF)*cutWeight);
    (allsamples.collection)[SampleIndices[iSample]].events->Draw("sumJetPt[1]/pt:pt>>+emuResponse",(zptforresponsepresentation&&cutOSOF)*cutWeight);
  }
  niceresponseplotd->Add(emuResponse,-1);
  
  niceresponseplotd->SetStats(0);
  niceresponseplotd->GetXaxis()->SetTitle("Z p_{T} [GeV]");
  niceresponseplotd->GetYaxis()->SetTitle("Response");
  niceresponseplotd->GetXaxis()->CenterTitle();
  niceresponseplotd->GetYaxis()->CenterTitle();
  niceresponseplotd->Draw("COLZ");
  TProfile * profd = (TProfile*)niceresponseplotd->ProfileX();
  profd->Rebin(2);
  profd->SetMarkerSize(DataMarkerSize);
  profd->Fit("pol0","","same,e1",100,400);
  DrawPrelim();
  string stitle="Data";
  if(!Contains(FindKeyword,"Data")) stitle="MC";
  TText* title = write_text(0.5,0.7,stitle.c_str());
  title->SetTextAlign(12);
  title->Draw();
  TF1 *datapol=(TF1*)profd->GetFunction("pol0");
  float correction=datapol->GetParameter(0);
  stringstream resstring;
  resstring<<"Response: "<<std::setprecision(2)<<100*correction<<" %";
  TText* restitle = write_text(0.5,0.65,resstring.str());
  restitle->SetTextAlign(12);
  restitle->SetTextSize(0.03);
  restitle->Draw();
  CompleteSave(cancorr,"ResponseCorrection/Response_Correction_Illustration_New_"+SaveAs);
  delete cancorr;
  delete niceresponseplotd;
  delete profd;
  return correction;
}

void find_correction_factors(string &jzbvardata,string &jzbvarmc) {
  
  dout << "Computing response corrections: " << endl;
  //Step 1 : Get results
  float datacorrection=find_one_correction_factor("Data",true,"","Data");
  float mccorrection=find_one_correction_factor("DY",false,"","MC");
  
  float dataEEcorrection=find_one_correction_factor("Data",true,"id1==0","Data_ee");
  float mcEEcorrection=find_one_correction_factor("DY",false,"id1==0","MC_ee");
  
  float dataMMcorrection=find_one_correction_factor("Data",true,"id1==1","Data_mm");
  float mcMMcorrection=find_one_correction_factor("DY",false,"id1==1","MC_mm");
  
  cout << "Corrections : " << endl;
  cout << "   Data : " << datacorrection << endl;
  cout << "     ee (" << dataEEcorrection << ") , mm (" << dataMMcorrection << ")" << endl;
  cout << "   MC : " << mccorrection << endl;
  cout << "     ee (" << mcEEcorrection << ") , mm (" << mcMMcorrection << ")" << endl;
  
  //Step 2: Processing the result and making it into something useful :-)
  stringstream jzbvardatas;
  jzbvardatas << "(";
  
  if(dataEEcorrection>=1) jzbvardatas<<"((id1==0&&id1==id2)*(jzb[1]-" << dataEEcorrection-1 << "*pt))";
  if(dataEEcorrection<1)  jzbvardatas<<"((id1==0&&id1==id2)*(jzb[1]+" << 1-dataEEcorrection << "*pt))";
  
  if(dataMMcorrection>=1) jzbvardatas<<"+((id1==1&&id1==id2)*(jzb[1]-" << dataMMcorrection-1 << "*pt))";
  if(dataMMcorrection<1)  jzbvardatas<<"+((id1==1&&id1==id2)*(jzb[1]+" << 1-dataMMcorrection << "*pt))";
  
  float averagecorrection=(dataMMcorrection+dataEEcorrection)/2.0;
  
  if(datacorrection>=1) jzbvardatas<<"+((id1!=id2)*(jzb[1]-" << datacorrection-1 << "*pt))";
  if(datacorrection<1) jzbvardatas<<"+((id1!=id2)*(jzb[1]+" << 1-datacorrection << "*pt))";
  
  jzbvardatas << ")";
  jzbvardata=jzbvardatas.str();
  
  stringstream jzbvarmcs;
  jzbvarmcs << "(";
  
  if(mcEEcorrection>=1) jzbvarmcs<<"((id1==0&&id1==id2)*(jzb[1]-" << mcEEcorrection-1 << "*pt))";
  if(mcEEcorrection<1)  jzbvarmcs<<"((id1==0&&id1==id2)*(jzb[1]+" << 1-mcEEcorrection << "*pt))";
  
  if(mcMMcorrection>=1) jzbvarmcs<<"+((id1==1&&id1==id2)*(jzb[1]-" << mcMMcorrection-1 << "*pt))";
  if(mcMMcorrection<1)  jzbvarmcs<<"+((id1==1&&id1==id2)*(jzb[1]+" << 1-mcMMcorrection << "*pt))";
  
  float averagemccorrection=(mcMMcorrection+mcEEcorrection)/2.0;
  
  if(mccorrection>=1) jzbvarmcs<<"+((id1!=id2)*(jzb[1]-" << mccorrection-1 << "*pt))";
  if(mccorrection<1) jzbvarmcs<<"+((id1!=id2)*(jzb[1]+" << 1-mccorrection << "*pt))";
  
  jzbvarmcs << ")";
  jzbvarmc=jzbvarmcs.str();

  dout << "JZB Z pt correction summary : " << endl;
  dout << "  Data: The response is " << datacorrection << "  --> jzb variable is now : " << jzbvardata << endl;
  dout << "  MC  : The response is " << mccorrection << "  --> jzb variable is now : " << jzbvarmc << endl;
  
}

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

void save_template(string mcjzb, string datajzb,vector<float> jzb_cuts,float MCPeakError,float DataPeakError, vector<float> jzb_shape_limit_bins) {
  dout << "Saving configuration template!" << endl;
  ofstream configfile;
  configfile.open("../DistributedModelCalculations/last_configuration.C");
  configfile<<"#include <iostream>\n";
  configfile<<"#include <vector>\n";
  configfile<<"#ifndef SampleClassLoaded\n";
  configfile<<"#include \"SampleClass.C\"\n";
  configfile<<"#endif\n";
  configfile<<"#define SetupLoaded\n";
  configfile<<"#ifndef ResultLibraryClassLoaded\n";
  configfile<<"#include \"ResultLibraryClass.C\"\n";
  configfile<<"#endif\n";

  configfile<<"\nusing namespace std;\n\n";
  
  configfile<<"namespace PlottingSetup { \n";
  configfile<<"string datajzb=\"datajzb_ERROR\";\n";
  configfile<<"string mcjzb=\"mcjzb_ERROR\";\n";
  configfile<<"vector<float>jzb_cuts;\n";
  configfile<<"vector<float>jzb_shape_limit_bins;\n";
  configfile<<"float MCPeakError=-999;\n";
  configfile<<"float DataPeakError=-999;\n";
  configfile<<"}\n\n";

  configfile<<"void read_config() {\n";
  configfile<<"datajzb=\""<<datajzb<<"\";\n";
  configfile<<"mcjzb=\""<<mcjzb<<"\";\n\n";
  configfile<<"\n\nMCPeakError="<<MCPeakError<<";\n";
  configfile<<"DataPeakError="<<DataPeakError<<";\n\n";
  for(int i=0;i<(int)jzb_cuts.size();i++) configfile<<"jzb_cuts.push_back("<<jzb_cuts[i]<<"); // JZB cut at " << jzb_cuts[i] << "\n";
  configfile<<"\n\n";
  for(int i=0;i<(int)Nobs.size();i++) configfile<<"Nobs.push_back("<<Nobs[i]<<"); // JZB cut at " << jzb_cuts[i] << "\n";
  for(int i=0;i<(int)Npred.size();i++) configfile<<"Npred.push_back("<<Npred[i]<<"); // JZB cut at " << jzb_cuts[i] << "\n";
  for(int i=0;i<(int)Nprederr.size();i++) configfile<<"Nprederr.push_back("<<Nprederr[i]<<"); // JZB cut at " << jzb_cuts[i] << "\n";
  configfile<<"\n\n";
  for(int i=0;i<(int)flippedNobs.size();i++) configfile<<"flippedNobs.push_back("<<flippedNobs[i]<<"); // JZB cut at " << jzb_cuts[i] << "\n";
  for(int i=0;i<(int)flippedNpred.size();i++) configfile<<"flippedNpred.push_back("<<flippedNpred[i]<<"); // JZB cut at " << jzb_cuts[i] << "\n";
  for(int i=0;i<(int)flippedNprederr.size();i++) configfile<<"flippedNprederr.push_back("<<flippedNprederr[i]<<"); // JZB cut at " << jzb_cuts[i] << "\n";
  for(int i=0;i<(int)jzb_shape_limit_bins.size();i++) configfile<<"jzb_shape_limit_bins.push_back("<<jzb_shape_limit_bins[i]<<"); // JZB shape bin boundary at " << jzb_shape_limit_bins[i] << "\n";
  configfile<<"\n\n";
  configfile<<"\n\n";
  configfile<<"luminosity="<<luminosity<<";\n";
  configfile<<"RestrictToMassPeak="<<RestrictToMassPeak<<";//defines the type of analysis we're running\n";
  
  configfile<<"\n\ncout << \"Configuration successfully loaded!\" << endl; \n \n } \n \n";
  
  configfile.close();

}

float get_nonzero_minimum(TH1F *histo) {
  float min=histo->GetMaximum();
  for(int ibin=1;ibin<=histo->GetNbinsX();ibin++) {
    float curcont=histo->GetBinContent(ibin);
    if(curcont<min&&curcont>0) min=curcont;
  }
  return min;
}

void draw_all_ttbar_histos(TCanvas *can, vector<TH1F*> histos, string drawoption="", float manualmin=-9) {
  can->cd();
  float min=1;
  float max=histos[0]->GetMaximum();
  if(manualmin>=0) min=manualmin;
  else {
    for(int i=1;i<(int)histos.size();i++) {
      float curmin=get_nonzero_minimum(histos[i]);
      float curmax=histos[i]->GetMaximum();
      if(curmin<min) min=curmin;
      if(curmax>max) max=curmax;
    }
  }
  histos[0]->GetYaxis()->SetRangeUser(min,4*max);
  histos[0]->Draw(drawoption.c_str());
  stringstream drawopt;
  drawopt << drawoption << ",same";
  for(int i=1;i<(int)histos.size();i++) {
    histos[i]->Draw(drawopt.str().c_str());
  }
}

void ttbar_sidebands_comparison(string mcjzb, vector<float> binning, string prestring) {
  //in the case of the on peak analysis, we compare the 3 control regions to the real value
  //in the case of the OFF peak analysis, we compare our control region to the real value
  TCut weightbackup=cutWeight;
  
  bool doPURW=false;
  
  
  if(!doPURW) {
    dout << "Not doing PU reweighting for ttbar closure test" << endl;
    cutWeight="1.0";
    // Do it without PU re-weighting
    float MCPeakNoPU=0,MCPeakErrorNoPU=0,DataPeakNoPU=0,DataPeakErrorNoPU=0,MCSigma=0,DataSigma=0;
    stringstream resultsNoPU;
    stringstream noPUdatajzb;
    stringstream noPUmcjzb;
    
    stringstream mcjzbnoPU;
    find_peaks(MCPeakNoPU,MCPeakErrorNoPU, DataPeakNoPU, DataPeakErrorNoPU,resultsNoPU,true,noPUdatajzb,noPUmcjzb);
    mcjzb = noPUmcjzb.str();
  }

      
  float simulatedlumi = luminosity; //in pb please - adjust to your likings

  TH1F *TZem  = systsamples.Draw("TZem",     mcjzb,binning,"JZB [GeV]","events",cutmass&&cutOSOF&&cutnJets,mc,simulatedlumi,systsamples.FindSample("TTJets"));
  TH1F *nTZem = systsamples.Draw("nTZem","-"+mcjzb,binning,"JZB [GeV]","events",cutmass&&cutOSOF&&cutnJets,mc,simulatedlumi,systsamples.FindSample("TTJets"));
  TH1F *TSem;
  TH1F *nTSem;
  TH1F *TZeemm  = systsamples.Draw("TZeemm",     mcjzb,binning,"JZB [GeV]","events",cutmass&&cutOSSF&&cutnJets,mc,simulatedlumi,systsamples.FindSample("TTJets"));
  TH1F *nTZeemm = systsamples.Draw("nTZeemm","-"+mcjzb,binning,"JZB [GeV]","events",cutmass&&cutOSSF&&cutnJets,mc,simulatedlumi,systsamples.FindSample("TTJets"));
  TH1F *TSeemm;
  TH1F *nTSeemm;
  
  if(PlottingSetup::RestrictToMassPeak) {
    TSem    = systsamples.Draw("TSem",       mcjzb,binning,"JZB [GeV]","events",sidebandcut&&cutOSOF&&cutnJets,mc,simulatedlumi,systsamples.FindSample("TTJets"));
    nTSem   = systsamples.Draw("nTSem",  "-"+mcjzb,binning,"JZB [GeV]","events",sidebandcut&&cutOSOF&&cutnJets,mc,simulatedlumi,systsamples.FindSample("TTJets"));
    TSeemm  = systsamples.Draw("TSeemm",     mcjzb,binning,"JZB [GeV]","events",sidebandcut&&cutOSSF&&cutnJets,mc,simulatedlumi,systsamples.FindSample("TTJets"));
    nTSeemm = systsamples.Draw("nTSeemm","-"+mcjzb,binning,"JZB [GeV]","events",sidebandcut&&cutOSSF&&cutnJets,mc,simulatedlumi,systsamples.FindSample("TTJets"));
  }

  TCanvas *tcan = new TCanvas("tcan","tcan");
  tcan->SetLogy(1);
  
  TZeemm->SetLineColor(kBlack);
  TZem->SetLineColor(kRed);
  TZeemm->SetMarkerColor(kBlack);
  TZem->SetMarkerColor(kRed);


  if(PlottingSetup::RestrictToMassPeak) {
        TSem->SetLineColor(TColor::GetColor("#00A616"));
        TSeemm->SetLineColor(kMagenta+2);
        TSem->SetMarkerColor(TColor::GetColor("#00A616"));
        TSeemm->SetMarkerColor(kMagenta+2);
        TSem->SetLineStyle(2);
        TSeemm->SetLineStyle(3);
  }
  
  vector<TH1F*> histos;
  TZem->GetXaxis()->SetRangeUser(-100,binning[binning.size()-1]);
  TZeemm->GetXaxis()->SetRangeUser(-100,binning[binning.size()-1]);
  histos.push_back(TZem);
  histos.push_back(TZeemm);
  if(PlottingSetup::RestrictToMassPeak) {
    TSeemm->GetXaxis()->SetRangeUser(-100,binning[binning.size()-1]);
    TSem->GetXaxis()->SetRangeUser(-100,binning[binning.size()-1]);
    histos.push_back(TSem);
    histos.push_back(TSeemm);
  }
  draw_all_ttbar_histos(tcan,histos,"histo",1);
  
  TLegend *leg = make_legend("MC t#bar{t}",0.6,0.65,false);
  leg->AddEntry(TZeemm,"SFZP","l");
  leg->AddEntry(TZem,"OFZP","l");
  if(PlottingSetup::RestrictToMassPeak) {
    leg->AddEntry(TSeemm,"SFSB","l");
    leg->AddEntry(TSem,"OFSB","l");
  }
  leg->Draw("same");
  DrawMCPrelim(simulatedlumi);
  CompleteSave(tcan,"Systematics/"+prestring+"/ttbar_shape_comparison");
  TH1F *TZemcopy = (TH1F*)TZem->Clone("TZemcopy");
  TH1F *TZeemmcopy = (TH1F*)TZeemm->Clone("TZeemmcopy");
  TH1F *TSeemmcopy;
  TH1F *TSemcopy;
  if(PlottingSetup::RestrictToMassPeak) {
    TSeemmcopy = (TH1F*)TSeemm->Clone("TSeemmcopy");
    TSemcopy = (TH1F*)TSem->Clone("TSemcopy");
  }

  TZem->Divide(TZeemm);
  TZem->SetMarkerStyle(21);
  if(PlottingSetup::RestrictToMassPeak) {
    TSem->Divide(TZeemm);
    TSeemm->Divide(TZeemm);
    TSem->SetMarkerStyle(24);
    TSeemm->SetMarkerStyle(32);
  }

  tcan->SetLogy(0);
  TZem->GetYaxis()->SetRangeUser(0,2.5);
  TZem->GetYaxis()->SetTitle("ratio");
  TZem->Draw();
  if(PlottingSetup::RestrictToMassPeak) {
    TSem->Draw("same");
    TSeemm->Draw("same");
  }
  
  float linepos=emuncertONPEAK;
  if(!PlottingSetup::RestrictToMassPeak) linepos=emuncertOFFPEAK;
  
  TLine *top = new TLine(binning[0],1.0+linepos,binning[binning.size()-1],1.0+linepos);
  TLine *center = new TLine(binning[0],1.0,binning[binning.size()-1],1.0);
  TLine *bottom = new TLine(binning[0],1.0-linepos,binning[binning.size()-1],1.0-linepos);
  
  /*write_warning(__FUNCTION__,"These two lines are to be removed!");
  TLine *topalt = new TLine(binning[0],1.0+0.1,binning[binning.size()-1],1.0+0.1);
  TLine *bottomalt = new TLine(binning[0],1.0-0.1,binning[binning.size()-1],1.0-0.1);
  topalt->SetLineColor(kRed);topalt->SetLineStyle(3);
  bottomalt->SetLineColor(kRed);bottomalt->SetLineStyle(3);
  topalt->Draw("same");bottomalt->Draw("same");*/
  
  
  top->SetLineColor(kBlue);top->SetLineStyle(2);
  bottom->SetLineColor(kBlue);bottom->SetLineStyle(2);
  center->SetLineColor(kBlue);
  
  top->Draw("same");
  center->Draw("same");
  bottom->Draw("same");
  
  TLegend *leg2 = make_legend("MC t#bar{t}",0.55,0.75,false);
  leg2->AddEntry(TZem,"OFZP / SFZP","ple");
  if(PlottingSetup::RestrictToMassPeak) {
    leg2->AddEntry(TSeemm,"SFSB / SFZP","ple");
    leg2->AddEntry(TSem,"OFSB / SFZP","ple");
  }
  leg2->AddEntry(bottom,"syst. envelope","l");
  leg2->SetX1(0.25);leg2->SetX2(0.6);
  leg2->SetY1(0.65);
  
  leg2->Draw("same");
  
  DrawMCPrelim(simulatedlumi);
  CompleteSave(tcan,"Systematics/"+prestring+"/ttbar_shape_comparison_ratio");
  
  if (0) { // Turn this off: we don't need this

    ///-------------- second part: only look at the quantity we actually care about!
    TH1F *leftsfzp  = (TH1F*)nTZeemm->Clone("leftsfzp");
    TH1F *rightsfzp = (TH1F*)TZeemmcopy->Clone("rightsfzp");
    rightsfzp->Add(leftsfzp,-1);
    TH1F *leftofzp  = (TH1F*)nTZem->Clone("leftofzp");
    TH1F *rightofzp = (TH1F*)TZemcopy->Clone("rightofzp");
    rightofzp->Add(leftofzp,-1);
    TH1F *leftofsb;
    TH1F *rightofsb;
    TH1F *leftsfsb;
    TH1F *rightsfsb;
    if(PlottingSetup::RestrictToMassPeak) {
      leftofsb  = (TH1F*)nTSem->Clone("leftofsb");
      rightofsb = (TH1F*)TSemcopy->Clone("rightofsb");
      rightofsb->Add(leftofsb,-1);
      leftsfsb  = (TH1F*)nTSeemm->Clone("leftsfsb");
      rightsfsb = (TH1F*)TSeemmcopy->Clone("rightsfsb");
      rightsfsb->Add(leftsfsb,-1);
    }

    tcan->SetLogy(1);
    rightsfzp->GetXaxis()->SetRangeUser(0,binning[binning.size()-1]);
    rightsfzp->GetYaxis()->SetTitle("#deltaJZB / 25 GeV");
    rightsfzp->Draw("histo");
    rightofzp->Draw("histo,same");
    if(PlottingSetup::RestrictToMassPeak) {
      rightofsb->Draw("histo,same");
      rightsfsb->Draw("histo,same");
    }
    DrawMCPrelim(simulatedlumi);

    TLegend *legA = make_legend("MC t#bar{t}",0.6,0.65,false);
    legA->AddEntry(rightsfzp,"SFZP","l");
    legA->AddEntry(rightofzp,"OFZP","l");
    if(PlottingSetup::RestrictToMassPeak) {
      legA->AddEntry(rightofsb,"SFSB","l");
      legA->AddEntry(rightsfsb,"OFSB","l");
    }
    legA->Draw();

    CompleteSave(tcan,"Systematics/"+prestring+"/ttbar_deltajzb_comparison");

    tcan->SetLogy(0);
    rightofzp->Divide(rightsfzp);
    rightofzp->GetXaxis()->SetRangeUser(0.0,binning[binning.size()-1]);
    rightofzp->GetYaxis()->SetRangeUser(0.0,2.5);
    rightofzp->GetYaxis()->SetTitle("#deltaJZB ratio");
    rightofzp->Draw();
    if(PlottingSetup::RestrictToMassPeak) {
      rightofsb->Divide(rightsfzp);
      rightsfsb->Divide(rightsfzp);
      rightofsb->Draw("same");
      rightsfsb->Draw("same");
    }

    TLine *top2 = new TLine(0.0,1.0+linepos,binning[binning.size()-1],1.0+linepos);
    TLine *center2 = new TLine(0.0,1.0,binning[binning.size()-1],1.0);
    TLine *bottom2 = new TLine(0.0,1.0-linepos,binning[binning.size()-1],1.0-linepos);
  
    top2->SetLineColor(kBlue);top2->SetLineStyle(2);
    bottom2->SetLineColor(kBlue);bottom2->SetLineStyle(2);
    center2->SetLineColor(kBlue);
  
    top2->Draw("same");
    center2->Draw("same");
    bottom2->Draw("same");

    rightofzp->SetMarkerStyle(21);
    if(PlottingSetup::RestrictToMassPeak) {
      rightofsb->SetMarkerStyle(24);
      rightsfsb->SetMarkerStyle(32);
    }

    TLegend *leg3 = make_legend("MC t#bar{t}",0.55,0.75,false);
    leg3->AddEntry(rightofzp,"OFZP / SFZP","ple");
    if(PlottingSetup::RestrictToMassPeak) {
      leg3->AddEntry(rightsfsb,"SFSB / SFZP","ple");
      leg3->AddEntry(rightofsb,"OFSB / SFZP","ple");
    }
    leg3->AddEntry(bottom,"syst. envelope","l");
    leg3->SetX1(0.25);leg3->SetX2(0.6);
    leg3->SetY1(0.65);
  
    leg3->Draw("same");
  
    DrawMCPrelim(simulatedlumi);
    CompleteSave(tcan,"Systematics/"+prestring+"/ttbar_deltajzb_comparison_ratio");

  }

  delete TZem;
  delete nTZem;
  delete TZeemm;
  delete nTZeemm;
  if(PlottingSetup::RestrictToMassPeak) {
        delete TSem;
        delete nTSem;
        delete TSeemm;
        delete nTSeemm;
  }

  delete tcan;
  cutWeight=weightbackup;
}

void ttbar_sidebands_comparison(string mcjzb, vector<float> jzb_binning) {
  vector<float> nicer_binning;
  
/*  nicer_binning.push_back(-400);
  nicer_binning.push_back(-250);
  nicer_binning.push_back(-200);
  nicer_binning.push_back(-150);
  nicer_binning.push_back(-100);
  nicer_binning.push_back(-50);
  nicer_binning.push_back(-20);
  
  nicer_binning.push_back(0);
  nicer_binning.push_back(20);
  nicer_binning.push_back(50);
  nicer_binning.push_back(100);
  nicer_binning.push_back(150);
  nicer_binning.push_back(200);
  nicer_binning.push_back(250);
  nicer_binning.push_back(400);*/
  
  nicer_binning.push_back(-100);
  nicer_binning.push_back(-50);
  nicer_binning.push_back(-25);
  nicer_binning.push_back(0);
  nicer_binning.push_back(25);
  nicer_binning.push_back(50);
  nicer_binning.push_back(75);
  nicer_binning.push_back(100);
  nicer_binning.push_back(125);
  nicer_binning.push_back(150);
  //nicer_binning.push_back(175);
  nicer_binning.push_back(200);
  nicer_binning.push_back(250);
  nicer_binning.push_back(300);
  nicer_binning.push_back(400);
  
  ttbar_sidebands_comparison(mcjzb,nicer_binning, "ttbar/");
}


void zjets_prediction_comparison(string mcjzbWithPUa) {
  TCanvas *zcan = new TCanvas("zcan","zcan");  
//  zcan->SetLogy(1);
  TCut weightbackup=cutWeight;
  
  bool UsePURW=false;
  
  
  string mcjzb;
  if(UsePURW) {
    mcjzb=mcjzbWithPUa;
  } else {
    // Do it without PU re-weighting
    float MCPeakNoPU=0,MCPeakErrorNoPU=0,DataPeakNoPU=0,DataPeakErrorNoPU=0,MCSigma=0,DataSigma=0;
    stringstream resultsNoPU;
    stringstream noPUdatajzb;
    stringstream noPUmcjzb;
    
    find_peaks(MCPeakNoPU,MCPeakErrorNoPU, DataPeakNoPU, DataPeakErrorNoPU,resultsNoPU,false,noPUdatajzb,noPUmcjzb);
    dout << "The peak corrected JZB expression for MC without pileup is : " <<  noPUmcjzb.str() << endl;
    
    mcjzb = noPUmcjzb.str();
    
    cutWeight="1.0";
  }

  
  vector<float> binning;
  binning.push_back(0);
  binning.push_back(10);
  binning.push_back(20);
  binning.push_back(40);
  binning.push_back(60);
//   binning.push_back(50);
//   binning.push_back(60);
//   binning.push_back(70);
//   binning.push_back(80);
//   binning.push_back(90);
  binning.push_back(100);

  float simulatedlumi = luminosity;//in pb please - adjust to your likings
  
  TCut kPos((mcjzb+">0").c_str());
  TCut kNeg((mcjzb+"<0").c_str());
  string var( "abs("+mcjzb+")" );
  
  TCut notTau("abs(genMID1)!=15");
  TCut ONLYTau("mll>20");
  

  TH1F *hJZBpos = systsamples.Draw("hJZBpos",var,binning, "JZB [GeV]", "events",cutmass&&cutOSSF&&cutnJets&&notTau&&kPos,mc,simulatedlumi,systsamples.FindSample("/DY"));
  TH1F *hJZBneg = systsamples.Draw("hJZBneg",var,binning, "JZB [GeV]", "events",cutmass&&cutOSSF&&cutnJets&&notTau&&kNeg,mc,simulatedlumi,systsamples.FindSample("/DY"));
  
  hJZBpos->SetLineColor(kBlack);
  hJZBneg->SetLineColor(kRed);
  
  hJZBpos->SetMinimum(1.0);
  hJZBpos->Draw("e1");
  hJZBneg->Draw("same,hist");
  hJZBpos->Draw("same,e1"); // So it's on top...
  
  TLegend *leg = make_legend("MC Z+jets #rightarrow ee,#mu#mu",0.55,0.75,false);
  leg->AddEntry(hJZBpos,"Observed","pe");
  leg->AddEntry(hJZBneg,"Predicted","l");
  leg->Draw("same");
  DrawMCPrelim(simulatedlumi);
  CompleteSave(zcan,"Systematics/ZJets/zjets_eemm_prediction");
  
  TH1F* hratio = (TH1F*)hJZBpos->Clone("hratio");
  hratio->Divide(hJZBneg);
  
  for(int i=1;i<=hJZBpos->GetNbinsX();i++) {
    cout << "Positive: " << hJZBpos->GetBinContent(i) << "   vs    Negative : " << hJZBneg->GetBinContent(i) << "   (ratio : " << hJZBpos->GetBinContent(i) / hJZBneg->GetBinContent(i) << endl;
  }
  
  zcan->SetLogy(0);
  hratio->GetYaxis()->SetRangeUser(0,2.0);
  hratio->GetYaxis()->SetTitle("Observed/Predicted");
  hratio->Draw("e1");
  
  TLine *top = new TLine(binning[0],1.25,binning[binning.size()-1],1.25);
  TLine *center = new TLine(binning[0],1.0,binning[binning.size()-1],1.0);
  TLine *bottom = new TLine(binning[0],0.75,binning[binning.size()-1],0.75);
  

  top->SetLineColor(kBlue);top->SetLineStyle(2);
  bottom->SetLineColor(kBlue);bottom->SetLineStyle(2);
  center->SetLineColor(kBlue);
  
  top->Draw("same");
  center->Draw("same");
  bottom->Draw("same");
  
  TLegend *leg2 = make_legend("MC Z+jets #rightarrow ee,#mu#mu",0.25,0.75,false);
  leg2->AddEntry(hratio,"obs / pred","pe");
  leg2->AddEntry(bottom,"syst. envelope","l");
  leg2->Draw("same");
  DrawMCPrelim(simulatedlumi);
  CompleteSave(zcan,"Systematics/ZJets/zjets_eemm_prediction_ratio");
  
  TCut reducedNJets(cutnJets);
  
  TH1F *TAUhJZBpos = systsamples.Draw("TAUhJZBpos",var,binning, "JZB [GeV]", "events",cutmass&&cutOSSF&&reducedNJets&&ONLYTau&&kPos,mc,simulatedlumi,systsamples.FindSample("/DY"));
  TH1F *LcorrJZBeemm     = systsamples.Draw("LcorrJZBeemm",var,binning, "JZB [GeV]", "events",cutmass&&cutOSSF&&reducedNJets&&ONLYTau&&kNeg,mc,simulatedlumi,systsamples.FindSample("/DY"));
  TH1F *RcorrJZBem       = systsamples.Draw("RcorrJZBem",var,binning, "JZB [GeV]", "events",cutmass&&cutOSOF&&reducedNJets&&ONLYTau&&kPos,mc,simulatedlumi,systsamples.FindSample("/DY"));
  TH1F *LcorrJZBem       = systsamples.Draw("LcorrJZBem",var,binning, "JZB [GeV]", "events",cutmass&&cutOSOF&&reducedNJets&&ONLYTau&&kNeg,mc,simulatedlumi,systsamples.FindSample("/DY"));

  TH1F *RcorrJZBSBem;
  TH1F *LcorrJZBSBem;
  TH1F *RcorrJZBSBeemm;
  TH1F *LcorrJZBSBeemm;

  if(PlottingSetup::RestrictToMassPeak) {
    RcorrJZBSBem   = systsamples.Draw("RcorrJZBSBem",var,binning, "JZB [GeV]", "events",sidebandcut&&cutOSOF&&reducedNJets&&ONLYTau&&kPos,mc,simulatedlumi,systsamples.FindSample("/DY"));
    LcorrJZBSBem   = systsamples.Draw("LcorrJZBSBem",var,binning, "JZB [GeV]", "events",sidebandcut&&cutOSOF&&reducedNJets&&ONLYTau&&kNeg,mc,simulatedlumi,systsamples.FindSample("/DY"));
    RcorrJZBSBeemm = systsamples.Draw("RcorrJZBSBeemm",var,binning, "JZB [GeV]", "events",sidebandcut&&cutOSSF&&reducedNJets&&ONLYTau&&kPos,mc,simulatedlumi,systsamples.FindSample("/DY"));
    LcorrJZBSBeemm = systsamples.Draw("LcorrJZBSBeemm",var,binning, "JZB [GeV]", "events",sidebandcut&&cutOSSF&&reducedNJets&&ONLYTau&&kNeg,mc,simulatedlumi,systsamples.FindSample("/DY"));
  }
  
  TH1F *Bpred = (TH1F*)LcorrJZBeemm->Clone("Bpred");
  if(PlottingSetup::RestrictToMassPeak) {
    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.);
  } else {
    Bpred->Add(RcorrJZBem,1.0);
    Bpred->Add(LcorrJZBem,-1.0);
  }
  
  Bpred->SetLineColor(kRed);

  TAUhJZBpos->SetLineColor(kBlack);
  Bpred->SetLineColor(kRed);
  
  TAUhJZBpos->SetMinimum(1.0);
  TAUhJZBpos->Draw("e1");
  Bpred->Draw("same,hist");
  TAUhJZBpos->Draw("same,e1");
  
  TLegend *TAUleg = make_legend("MC Z+jets #rightarrow ee,#mu#mu,#tau#tau",0.55,0.75,false);
  TAUleg->AddEntry(TAUhJZBpos,"Observed","pe");
  TAUleg->AddEntry(Bpred,"Predicted","l");
  TAUleg->Draw("same");
  DrawMCPrelim(simulatedlumi);
  CompleteSave(zcan,"Systematics/ZJets/zjets_eemumutautau_prediction");
  
  TH1F* TAUhratio = (TH1F*)TAUhJZBpos->Clone("TAUhratio");
  TAUhratio->Divide(Bpred);
  
  for(int i=1;i<=TAUhJZBpos->GetNbinsX();i++) {
    cout << "ee/mm/tautau observed: " << TAUhJZBpos->GetBinContent(i) << "   vs    predicted : " << Bpred->GetBinContent(i) << "   (ratio : " << TAUhJZBpos->GetBinContent(i) / Bpred->GetBinContent(i) << endl;
  }
  
  zcan->SetLogy(0);
  TAUhratio->GetYaxis()->SetRangeUser(0,2.5);
  TAUhratio->GetYaxis()->SetTitle("Observed/Predicted");
  TAUhratio->Draw("e1");
  
  top->Draw("same");
  center->Draw("same");
  bottom->Draw("same");
  
  TLegend *TAUleg2 = make_legend("MC Z+jets #rightarrow #tau#tau",0.25,0.75,false);
  TAUleg2->AddEntry(TAUhratio,"obs / pred","pe");
  TAUleg2->AddEntry(bottom,"syst. envelope","l");
  TAUleg2->Draw("same");
  DrawMCPrelim(simulatedlumi);
  CompleteSave(zcan,"Systematics/ZJets/zjets_eemumutautau_prediction_ratio");
  
  delete Bpred;
  delete TAUhJZBpos;
  delete LcorrJZBeemm;
  delete RcorrJZBem;
  delete LcorrJZBem;
  
  if(PlottingSetup::RestrictToMassPeak) {
    delete RcorrJZBSBem;
    delete LcorrJZBSBem;
    delete RcorrJZBSBeemm;
    delete LcorrJZBSBeemm;
  }
  
  
  delete zcan;
  cutWeight=weightbackup;

}


void sideband_assessment(string datajzb, float min=30.0, float max=50.0) {
  tout << endl << endl;
  stringstream bordercut;
  bordercut << "(TMath::Abs(" << datajzb << ")<" << max << ")&&(TMath::Abs(" << datajzb << ")>" << min << ")";
  tout << bordercut.str().c_str() << endl;
  TH1F *ofsb    = allsamples.Draw("ofsb",datajzb,100, 0,100, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets&&bordercut.str().c_str(),data, luminosity);
  TH1F *ofzp    = allsamples.Draw("ofzp",datajzb,100, 0,100, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets&&bordercut.str().c_str(),data, luminosity);
  float OFSB = ofsb->Integral();
  float OFZP = ofzp->Integral();

  tout << "\\begin{table}[hbtp]" << endl;
  tout << "\\renewcommand{\\arraystretch}{1.3}" << endl;
  tout << "\\begin{center}" << endl;
  tout << "\\caption{Total number of events observed in the auxiliary region  $|JZB|>"<<min<<"\\GeV \\cup |JZB|<"<<max<<"\\GeV$ for the {\\OFZP} and {\\OFSB}.}\\label{tab:auxcounts}" << endl;
  tout << "\\begin{tabular}{l|cc}" << endl;
  tout << "\\hline" << endl;
  tout << "& {\\OFZP} & {\\OFSB} \\\\\\hline" << endl;
  tout << "\\#(events) & "<<OFZP<<" & "<<OFSB<<"\\\\ \\hline" << endl;
  tout << "\\end{tabular}" << endl;
  tout << "\\end{center}" << endl;
  tout << "\\end{table}" << endl;


}

void make_table(samplecollection &coll, string jzbexpr, bool is_data, vector<float> jzb_cuts, string subselection="none") {
  
  vector<float> jzbcutprediction;
  vector<float> metcutprediction;

  vector<float> jzbcutobservation;
  vector<float> metcutobservation;
  
  TCanvas *cannie = new TCanvas("cannie","cannie");
  
  for(int icut=0;icut<(int)jzb_cuts.size();icut++) {
    float currcut=jzb_cuts[icut];
    int nbins=1;float low=currcut;
    vector<int> mysample;
    if(subselection!="none") mysample=coll.FindSample(subselection);
    TH1F *RcorrJZBeemm   = coll.Draw("RcorrJZBeemm",jzbexpr,1,currcut,14000, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,is_data, luminosity, mysample);
    TH1F *LcorrJZBeemm   = coll.Draw("LcorrJZBeemm",("-("+jzbexpr+")").c_str(),nbins,low,14000, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,is_data, luminosity, mysample);
    TH1F *RcorrJZBem     = coll.Draw("RcorrJZBem",jzbexpr.c_str(),nbins,low,14000, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,is_data, luminosity, mysample);
    TH1F *LcorrJZBem     = coll.Draw("LcorrJZBem",("-("+jzbexpr+")").c_str(),nbins,low,14000, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,is_data, luminosity, mysample);
    
    TH1F *RcorrJZBSBem;
    TH1F *LcorrJZBSBem;
    TH1F *RcorrJZBSBeemm;
    TH1F *LcorrJZBSBeemm;
    
    if(PlottingSetup::RestrictToMassPeak) {
      RcorrJZBSBem   = coll.Draw("RcorrJZBSBem",jzbexpr.c_str(),nbins,low,14000, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets,is_data, luminosity, mysample);
      LcorrJZBSBem   = coll.Draw("LcorrJZBSBem",("-("+jzbexpr+")").c_str(),nbins,low,14000, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets,is_data, luminosity, mysample);
      RcorrJZBSBeemm = coll.Draw("RcorrJZBSBeemm",jzbexpr.c_str(),nbins,low,14000, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets,is_data, luminosity, mysample);
      LcorrJZBSBeemm = coll.Draw("LcorrJZBSBeemm",("-("+jzbexpr+")").c_str(),nbins,low,14000, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets,is_data, luminosity, mysample);
    }
    
    TH1F *Bpred = (TH1F*)LcorrJZBeemm->Clone("Bpred");
    if(PlottingSetup::RestrictToMassPeak) {
      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.);
    } else {
      Bpred->Add(RcorrJZBem,1.0);
      Bpred->Add(LcorrJZBem,-1.0);
    }
    
    jzbcutobservation.push_back(RcorrJZBeemm->Integral());
    jzbcutprediction.push_back(Bpred->Integral());
    
    delete RcorrJZBeemm;
    delete LcorrJZBeemm;
    delete RcorrJZBem;
    delete LcorrJZBem;
    delete RcorrJZBSBem;
    delete LcorrJZBSBem;
    delete RcorrJZBSBeemm;
    delete LcorrJZBSBeemm;
    
    TH1F *MetObs         = coll.Draw("MetObs",("met[4]"),nbins,low,14000, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,is_data, luminosity, mysample);
    TH1F *MetPred        = coll.Draw("MetPred",("met[4]"),nbins,low,14000, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,is_data, luminosity, mysample);
    
    metcutobservation.push_back(MetObs->Integral());
    metcutprediction.push_back(MetPred->Integral());
    delete MetObs;
    delete MetPred;
  }//end of cut loop
    
  //prediction part
  if(is_data) cout << "Data prediction & ";
  if(subselection!="none") cout << subselection << " prediction &";
  for(int ij=0;ij<(int)jzb_cuts.size();ij++) cout << jzbcutprediction[ij] << " vs " << metcutprediction[ij] << " & ";
  
  cout << endl;
  //observation part
  if(is_data) cout << "Data observation & ";
  if(subselection!="none") cout << subselection << " observation &";
  for(int ij=0;ij<(int)jzb_cuts.size();ij++) cout << jzbcutobservation[ij] << " vs " << metcutobservation[ij] << " & ";
  cout << endl;
  cout << "_________________________________________________________________" << endl;
  delete cannie;
}

void met_jzb_cut(string datajzb, string mcjzb, vector<float> jzb_cut) {
  /*we want a table like this: 
    __________________     50   |   100  | ...
    | Data prediction |  a vs b | a vs b | ...
    | Data observed   |  a vs b | a vs b | ...
    --------------------------------------
    --------------------------------------
    | LM4 prediction  |  a vs b | a vs b | ...
    | LM4 observed    |  a vs b | a vs b | ...
    | LM8 prediction  |  a vs b | a vs b | ...
    | LM8 observed    |  a vs b | a vs b | ...
    
    where a is the result for a jzb cut at X, and b is the result for a met cut at X
*/
  make_table(allsamples,datajzb, true,jzb_cut,"none");
  make_table(signalsamples,mcjzb, false,jzb_cut,"LM4");
  make_table(signalsamples,mcjzb, false,jzb_cut,"LM8");
}
  

//________________________________________________________________________________________
// JZB Efficiency plot (efficiency of passing reco. JZB cut as function of generator JZB cut)
void JZBSelEff(string mcjzb, TTree* events, string informalname, vector<float> jzb_bins) {

  float min = 0, max = 400;
  float xbins[] = {min,10,20,30,40,50,60,70,80,90,100,110,120,140,150,160,170,180,190,200,210,220,250,300,max,max+100};
  int nbins = sizeof(xbins)/sizeof(float)-1;
  int markers[] = { 20, 26, 21, 24, 22, 25, 28 };

  
  TH1F* heff  = new TH1F("heff", "JZB eff ; generator JZB [GeV]; efficiency",nbins,xbins);
  TH1F* hgen  = new TH1F("hgen", "JZB gen ; generator JZB [GeV]; efficiency",nbins,xbins);
  TH1F* hreco = new TH1F("hreco","JZB reco ; generator JZB [GeV]; efficiency",nbins,xbins);

  TCut kgen(genMassCut&&"genZPt>0&&genNjets>2&&abs(genMID)==23"&&cutOSSF);
  TCut kreco(cutmass);

  TF1* func = new TF1("func","0.5*[2]*(TMath::Erf((x-[1])/[0])+1)",min,max);
  func->SetParNames("epsilon","x_{1/2}","sigma");
  func->SetParameter(0,50.);
  func->SetParameter(1,0.);
  func->SetParameter(2,1.);
  gStyle->SetOptStat(0);
  gStyle->SetOptFit(0);

  TCanvas *can = new TCanvas("can","Canvas for JZB Efficiency",600,600);
  can->SetGridx(1);
  can->SetGridy(1);
  can->SetLeftMargin(0.16);
  can->SetRightMargin(0.05);
  TLegend *leg = make_legend("",0.6,0.2,false,0.89,0.5);
  leg->SetBorderSize(1);
  leg->SetLineColor(kBlack);
  leg->SetTextFont(62);

  for ( int icut=0; icut<(int)jzb_bins.size(); ++icut ) {

    ostringstream selection;
    selection << mcjzb << ">" << jzb_bins[icut];
    TCut ksel(selection.str().c_str());
    events->Draw("genJZB>>hgen", kgen&&kreco,      "goff");
    events->Draw("genJZB>>hreco",kgen&&kreco&&ksel,"goff");

    // Loop over steps to get efficiency curve
    for ( Int_t iBin = 0; iBin<nbins-1; ++iBin ) {
      Float_t eff = hreco->GetBinContent(iBin+1)/hgen->GetBinContent(iBin+1);
      heff->SetBinContent(iBin+1,eff);
      heff->SetBinError(iBin+1,TMath::Sqrt(eff*(1-eff)/hgen->GetBinContent(iBin+1)));
    }
  
    heff->GetXaxis()->SetRangeUser(min, max);
//    heff->GetXaxis()->SetLabelSize(0.05); // paper style. overruled.
//    heff->GetYaxis()->SetLabelSize(0.05); // paper style. overruled.
//    heff->GetXaxis()->SetTitleSize(0.06); // paper style. overruled.
//    heff->GetYaxis()->SetTitleSize(0.06); // paper style. overruled.
    heff->SetMarkerStyle(markers[icut]);
    heff->Fit("func","Q+","same");

    // Print values
    dout << "+++ For " << selection.str() << std::endl;
    for ( int i=0; i<func->GetNpar(); ++i )
      dout << "   " << func->GetParName(i) << " " << func->GetParameter(i) << " \\pm " << func->GetParError(i) << std::endl;
    char hname[256]; sprintf(hname,"heff%d",icut);

    // Store plot
    TH1F* h = (TH1F*)heff->Clone(hname);
    h->SetNdivisions(505,"X");
    if ( icut) h->Draw("same");
    else h->Draw();
    char htitle[256]; sprintf(htitle,"JZB > %3.0f GeV", jzb_bins[icut]);
    leg->AddEntry(h,htitle,"p");

  }

  leg->Draw();
  DrawMCPrelim(0.0);
  CompleteSave(can, "Systematics/jzb_efficiency_curve"+informalname );

  delete hgen;
  delete hreco;
  delete heff;
}

//________________________________________________________________________________________
// Calls the above function for each signal sample
void plot_jzb_sel_eff(string mcjzb, samplecollection &signalsamples, vector<float> bins ) 
{  
  for (int isignal=0; isignal<(int)signalsamples.collection.size();isignal++) {
    dout << "JZB selection efficiency curve: " << std::endl;
    JZBSelEff(mcjzb,(signalsamples.collection)[isignal].events,(signalsamples.collection)[isignal].samplename,bins); // Only for some selected samples
  }
}

void qcd_plots(string datajzb, string mcjzb, vector<float> bins) {
  // What this function aims to do: 
  // Illustrate cut flow for QCD (requiring only one lepton, requiring etc.)
  // Illustrate how little QCD is left over! i.e. make some pointless JZB plot with only QCD to visualize the fact that there's not much really.
  TCanvas *can = new TCanvas("can","can");
  TPad *kinpad = new TPad("kinpad","kinpad",0,0,1,1);
  kinpad->cd();
  
  string jzb=mcjzb;
  
  float hi=400;
  bool use_signal=false;
  bool use_data=false;
  
  bool is_data=false;
  int nbins=50;//100;
  float low=0;
  float high=500;
  int versok=false;
  if(gROOT->GetVersionInt()>=53000) versok=true;

  TH1F *blankback = new TH1F("blankback","blankback",int(high/10),0,high);
  TH1F *RcorrJZBeemm   = qcdsamples.Draw("RcorrJZBeemm",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,is_data, luminosity,use_signal);
  TH1F *LcorrJZBeemm   = qcdsamples.Draw("LcorrJZBeemm",("-"+jzb).c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,is_data, luminosity,use_signal);
  TH1F *RcorrJZBem     = qcdsamples.Draw("RcorrJZBem",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,is_data, luminosity,use_signal);
  TH1F *LcorrJZBem     = qcdsamples.Draw("LcorrJZBem",("-"+jzb).c_str(),nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,is_data, luminosity,use_signal);
  blankback->GetXaxis()->SetTitle(RcorrJZBeemm->GetXaxis()->GetTitle());
  blankback->GetYaxis()->SetTitle(RcorrJZBeemm->GetYaxis()->GetTitle());
  blankback->GetXaxis()->CenterTitle();
  blankback->GetYaxis()->CenterTitle();
  
  flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak ---- prediction changed.
  TH1F *RcorrJZBSBem;
  TH1F *LcorrJZBSBem;
  TH1F *RcorrJZBSBeemm;
  TH1F *LcorrJZBSBeemm;
  
//  TH1F *RcorrJZBeemmNoS;

    //these are for the ratio
  TH1F *JRcorrJZBeemm   = qcdsamples.Draw("JRcorrJZBeemm",jzb.c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,is_data, luminosity,use_signal);
  TH1F *JLcorrJZBeemm   = qcdsamples.Draw("JLcorrJZBeemm",("-"+jzb).c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,is_data, luminosity,use_signal);
  TH1F *JRcorrJZBem     = qcdsamples.Draw("JRcorrJZBem",jzb.c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,is_data, luminosity,use_signal);
  TH1F *JLcorrJZBem     = qcdsamples.Draw("JLcorrJZBem",("-"+jzb).c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,is_data, luminosity,use_signal);
  
  TH1F *JRcorrJZBSBem;
  TH1F *JLcorrJZBSBem;
  TH1F *JRcorrJZBSBeemm;
  TH1F *JLcorrJZBSBeemm;
  
  if(PlottingSetup::RestrictToMassPeak) {
    RcorrJZBSBem   = qcdsamples.Draw("RcorrJZBSBem",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets,is_data, luminosity,use_signal);
    LcorrJZBSBem   = qcdsamples.Draw("LcorrJZBSBem",("-"+jzb).c_str(),nbins,low,hi, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets,is_data, luminosity,use_signal);
    RcorrJZBSBeemm = qcdsamples.Draw("RcorrJZBSBeemm",jzb.c_str(),nbins,low,hi, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets,is_data, luminosity,use_signal);
    LcorrJZBSBeemm = qcdsamples.Draw("LcorrJZBSBeemm",("-"+jzb).c_str(),nbins,low,hi, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets,is_data, luminosity,use_signal);

    //these are for the ratio
    JRcorrJZBSBem   = qcdsamples.Draw("JRcorrJZBSBem",jzb.c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets,is_data, luminosity,use_signal);
    JLcorrJZBSBem   = qcdsamples.Draw("JLcorrJZBSBem",("-"+jzb).c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets,is_data, luminosity,use_signal);
    JRcorrJZBSBeemm = qcdsamples.Draw("JRcorrJZBSBeemm",jzb.c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets,is_data, luminosity,use_signal);
    JLcorrJZBSBeemm = qcdsamples.Draw("JLcorrJZBSBeemm",("-"+jzb).c_str(),PlottingSetup::global_ratio_binning, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets,is_data, luminosity,use_signal);

  }
  
  flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak ---- prediction changed.
  
  TH1F *Zjetspred = (TH1F*)LcorrJZBeemm->Clone("Zjetspred");
  TH1F *TTbarpred = (TH1F*)RcorrJZBem->Clone("TTbarpred");
    
  TH1F *Bpred = (TH1F*)LcorrJZBeemm->Clone("Bpred");
  TH1F *JBpred = (TH1F*)JLcorrJZBeemm->Clone("Bpred");
  if(PlottingSetup::RestrictToMassPeak) {
    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.);
    
    TTbarpred->Scale(1.0/3);
    Zjetspred->Add(LcorrJZBem,-1.0/3.);
    Zjetspred->Add(LcorrJZBSBem,-1.0/3.);
    TTbarpred->Add(RcorrJZBSBem,1.0/3.);
    Zjetspred->Add(LcorrJZBSBeemm,-1.0/3.);
    TTbarpred->Add(RcorrJZBSBeemm,1.0/3.);
    
    //these are for the ratio
    JBpred->Add(JRcorrJZBem,1.0/3.);
    JBpred->Add(JLcorrJZBem,-1.0/3.);
    JBpred->Add(JRcorrJZBSBem,1.0/3.);
    JBpred->Add(JLcorrJZBSBem,-1.0/3.);
    JBpred->Add(JRcorrJZBSBeemm,1.0/3.);
    JBpred->Add(JLcorrJZBSBeemm,-1.0/3.);
  } else {
    Bpred->Add(RcorrJZBem,1.0);
    Bpred->Add(LcorrJZBem,-1.0);
    
    Zjetspred->Add(LcorrJZBem,-1.0);
    
    //these are for the ratio
    JBpred->Add(JRcorrJZBem,1.0);
    JBpred->Add(JLcorrJZBem,-1.0);
  }
    
  
  flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak ---- prediction changed
  
  TLegend *legBpred = make_legend("",0.6,0.55,false);
  RcorrJZBeemm->Draw("e1x0,same");
  Bpred->Draw("hist,same");
  RcorrJZBeemm->Draw("e1x0,same");//HAVE IT ON TOP!
  legBpred->AddEntry(RcorrJZBeemm,"MC observed","p");
  legBpred->AddEntry(Bpred,"MC predicted","l");
  if(versok) legBpred->AddEntry((TObject*)0,"",""); // Just for alignment // causes seg fault on root v5.18
  if(versok) legBpred->AddEntry((TObject*)0,"",""); // causes seg fault on root v5.18
  legBpred->Draw();
  DrawMCPrelim();

  //3rd last argument: do special bpred ratio, 2nd last argument: extended range!, last: y-axis title
  string ytitle("ratio");
  if ( use_data==1 ) ytitle = "data/pred";
  save_with_ratio(JRcorrJZBeemm,JBpred,kinpad,"QCD/Bpred",true,false,ytitle); 
  
  TH1F *allevents = qcdsamples.Draw("allevents","pfJetGoodNum",1,0,100, "internal code", "events", "" ,mc, luminosity);
  TH1F *ossf = qcdsamples.Draw("ossf","pfJetGoodNum",1,0,100, "internal code", "events", cutOSSF ,mc, luminosity);
  TH1F *osof = qcdsamples.Draw("osof","pfJetGoodNum",1,0,100, "internal code", "events", cutOSOF ,mc, luminosity);
  TH1F *njossf = qcdsamples.Draw("njossf","pfJetGoodNum",1,0,100, "internal code", "events", cutnJets&&cutOSSF ,mc, luminosity);
  TH1F *njosof = qcdsamples.Draw("njosof","pfJetGoodNum",1,0,100, "internal code", "events", cutnJets&&cutOSOF ,mc, luminosity);
  
  dout << "______________________________________________" << endl;
  dout << "QCD contribution: " << endl;
  dout << "Total number of events: " << allevents->Integral() << endl;
  dout << "OSSF events: " << ossf->Integral() << endl;
  dout << "OSOF events: " << osof->Integral() << endl;
  dout << "OSSF events with >=3 jets:" << njossf->Integral() << endl;
  dout << "OSOF events with >=3 jets:" << njosof->Integral() << endl;
  dout << "(note that no mass requirement has been imposed)" << endl;
  
  dout << "______________________________________________" << endl;
  dout << "How QCD shows up in the different regions: " << endl;
  dout << "OSSF: " << endl;
  if(PlottingSetup::RestrictToMassPeak) {
    dout << "     Z window: \t" << RcorrJZBeemm->Integral() << " (JZB>0) , " << LcorrJZBeemm->Integral() << " (JZB<0) --> total: " << RcorrJZBeemm->Integral() + LcorrJZBeemm->Integral() << endl;
    dout << "     sideband: \t" << RcorrJZBSBeemm->Integral() << " (JZB>0) , " << LcorrJZBSBeemm->Integral() << " (JZB<0) --> total: " << RcorrJZBSBeemm->Integral() + LcorrJZBSBeemm->Integral() << endl;
  } else {
    dout << "     " << RcorrJZBeemm->Integral() << " (JZB>0) , " << LcorrJZBeemm->Integral() << " (JZB<0) --> total: " << RcorrJZBeemm->Integral() + LcorrJZBeemm->Integral() << endl;
  }
  dout << "OSOF: " << endl;
  if(PlottingSetup::RestrictToMassPeak) {
    dout << "     Z window: \t" << RcorrJZBem->Integral() << " (JZB>0) , " << LcorrJZBem->Integral() << " (JZB<0) --> total: " << RcorrJZBem->Integral() + LcorrJZBem->Integral() << endl;
    dout << "     sideband: \t" << RcorrJZBSBem->Integral() << " (JZB>0) , " << LcorrJZBSBem->Integral() << " (JZB<0) --> total: " << RcorrJZBSBem->Integral() + LcorrJZBSBem->Integral() << endl;
  } else {
    dout << "     Z window: \t" << RcorrJZBem->Integral() << " (JZB>0) , " << LcorrJZBem->Integral() << " (JZB<0) --> total: " << RcorrJZBem->Integral() + LcorrJZBem->Integral() << endl;
  }
  dout << "Therefore: " << endl;
  if(PlottingSetup::RestrictToMassPeak) {
    dout << "    Prediction increases by : " << LcorrJZBeemm->Integral() << " + (1.0/3)*(" << RcorrJZBSBeemm->Integral() <<"-"<< LcorrJZBSBeemm->Integral() << ") (SFSB) ";
    dout << " + (1.0/3)*(" << RcorrJZBem->Integral() <<"-"<< LcorrJZBem->Integral() << ") (OFZP) ";
    dout << " + (1.0/3)*(" << RcorrJZBSBem->Integral() <<"-"<< LcorrJZBSBem->Integral() << ") (OFSB) ";
    dout << " = " << LcorrJZBeemm->Integral() + (1.0/3)*(RcorrJZBSBeemm->Integral() - LcorrJZBSBeemm->Integral() + RcorrJZBem->Integral() - LcorrJZBem->Integral() + RcorrJZBSBem->Integral() - LcorrJZBSBem->Integral()) << endl;
  } else {
    dout << "    Prediction increases by : " << LcorrJZBeemm->Integral();
    dout << " + (" << RcorrJZBem->Integral() <<"-"<< LcorrJZBem->Integral() << ") (OFZP) ";
    dout << " = " << LcorrJZBeemm->Integral() + RcorrJZBem->Integral() - LcorrJZBem->Integral() << endl;
  }
  dout << "    Observation increases by : " << RcorrJZBeemm->Integral() << endl;
  
  dout << endl;
  for(int i=0;i<(int)bins.size();i++) {
    dout << " JZB > " << bins[i] << " : " << endl;
    dout << "    Observation increases by : " << RcorrJZBeemm->Integral(RcorrJZBeemm->FindBin(bins[i]),RcorrJZBeemm->GetNbinsX()) << endl;
    if(PlottingSetup::RestrictToMassPeak) {
        dout << "    Prediction increases by : " << LcorrJZBeemm->Integral(RcorrJZBeemm->FindBin(bins[i]),RcorrJZBeemm->GetNbinsX()) + (1.0/3)*(RcorrJZBSBeemm->Integral(RcorrJZBeemm->FindBin(bins[i]),RcorrJZBeemm->GetNbinsX()) - LcorrJZBSBeemm->Integral(RcorrJZBeemm->FindBin(bins[i]),RcorrJZBeemm->GetNbinsX()) + RcorrJZBem->Integral(RcorrJZBeemm->FindBin(bins[i]),RcorrJZBeemm->GetNbinsX()) - LcorrJZBem->Integral(RcorrJZBeemm->FindBin(bins[i]),RcorrJZBeemm->GetNbinsX()) + RcorrJZBSBem->Integral(RcorrJZBeemm->FindBin(bins[i]),RcorrJZBeemm->GetNbinsX()) - LcorrJZBSBem->Integral(RcorrJZBeemm->FindBin(bins[i]),RcorrJZBeemm->GetNbinsX())) << endl;
    } else {
        dout << "    Prediction increases by : " << LcorrJZBeemm->Integral(RcorrJZBeemm->FindBin(bins[i]),RcorrJZBeemm->GetNbinsX()) + RcorrJZBem->Integral(RcorrJZBeemm->FindBin(bins[i]),RcorrJZBeemm->GetNbinsX()) - LcorrJZBem->Integral(RcorrJZBeemm->FindBin(bins[i]),RcorrJZBeemm->GetNbinsX()) << endl;
    }
  }
      
  delete can;
  delete allevents;
  if(ossf) delete ossf;
  if(RcorrJZBem) delete RcorrJZBem;
  if(LcorrJZBem) delete LcorrJZBem;
  if(RcorrJZBeemm) delete RcorrJZBeemm;
  if(LcorrJZBeemm) delete LcorrJZBeemm;
  if(PlottingSetup::RestrictToMassPeak&&RcorrJZBSBem) delete RcorrJZBSBem;
  if(PlottingSetup::RestrictToMassPeak&&LcorrJZBSBem) delete LcorrJZBSBem;
  if(PlottingSetup::RestrictToMassPeak&&RcorrJZBSBeemm) delete RcorrJZBSBeemm;
  if(PlottingSetup::RestrictToMassPeak&&LcorrJZBSBeemm) delete LcorrJZBSBeemm;
}

void check_ptsanity() {
  TCanvas *ptsancan = new TCanvas("ptsancan","ptsancan",600,1800);
  TH1F *individualpt1histos[allsamples.collection.size()];
  TH1F *individualpt2histos[allsamples.collection.size()];
  TH1F *fpt1 = new TH1F("fpt1","fpt1",50,0,50);
  fpt1->GetYaxis()->SetRangeUser(0,1);
  fpt1->GetXaxis()->SetTitle("p_{T,1}");
  fpt1->GetXaxis()->CenterTitle();

  TH1F *fpt2 = new TH1F("fpt2","fpt2",50,0,50);
  fpt2->GetXaxis()->SetTitle("p_{T,2}");
  fpt2->GetXaxis()->CenterTitle();
  
  ptsancan->Divide(1,3);
  ptsancan->cd(1);
  float maxpt1entry=0;
  float maxpt2entry=0;
  
  TLegend *leg = make_legend();
  leg->SetX1(0.0);
  leg->SetY1(0.0);
  leg->SetX2(1.0);
  leg->SetY2(1.0);
  
  
  for(int isample=0;isample<(int)allsamples.collection.size();isample++) {
    string nowname=(allsamples.collection)[isample].filename;
    cout << "Drawing: " << nowname << " (sample " << isample+1 << " / " << allsamples.collection.size() << ")" << endl;
    individualpt1histos[isample] = allsamples.Draw(GetNumericHistoName(),"pt1",50,0,50, "p_{T,1}", "events",cutOSSF&&cutnJets,mc,luminosity,allsamples.FindSample(nowname));
    individualpt2histos[isample] = allsamples.Draw(GetNumericHistoName(),"pt2",50,0,50, "p_{T,2}", "events",cutOSSF&&cutnJets,mc,luminosity,allsamples.FindSample(nowname));
    individualpt1histos[isample]->SetLineColor(isample+1);
    individualpt2histos[isample]->SetLineColor(isample+1);
    float currmaxpt1entry=individualpt1histos[isample]->GetMaximum()/individualpt1histos[isample]->Integral();
    float currmaxpt2entry=individualpt2histos[isample]->GetMaximum()/individualpt2histos[isample]->Integral();
    cout << "      pt 1 histo contains; " << individualpt1histos[isample]->Integral() << endl;
    cout << "      pt 2 histo contains; " << individualpt2histos[isample]->Integral() << endl;
    if(currmaxpt1entry>maxpt1entry)maxpt1entry=currmaxpt1entry; 
    if(currmaxpt2entry>maxpt2entry)maxpt2entry=currmaxpt2entry;
    leg->AddEntry(individualpt2histos[isample],((allsamples.collection)[isample].filename).c_str(),"f");
  }

  fpt1->GetYaxis()->SetRangeUser(0,maxpt1entry);
  fpt2->GetYaxis()->SetRangeUser(0,maxpt2entry);

  ptsancan->cd(1);
  fpt1->Draw();
  ptsancan->cd(2);
  fpt2->Draw();

  for(int isample=0;isample<(int)allsamples.collection.size();isample++) {
    ptsancan->cd(1);
    individualpt1histos[isample]->DrawNormalized("same,histo");
    ptsancan->cd(2);
    individualpt2histos[isample]->DrawNormalized("same,histo");
  }
  ptsancan->cd(3);
  leg->Draw();
  CompleteSave(ptsancan,"PtSanityCheck");
  
  delete ptsancan;
}

void do_mlls_plot(string mcjzb) {
  cout << "At this point we'd plot the mll distribution" << endl;
  TCanvas *sigcan = new TCanvas("sigcan","sigcan");
  for(int isig=0;isig<(int)(signalsamples.collection).size();isig++) {
    if(!(signalsamples.collection)[isig].events) continue;
    string nowname=(signalsamples.collection)[isig].filename;
    TH1F *mll = signalsamples.Draw("mllhisto","mll",150,0,150, "m_{ll}", "events",cutOSSF&&cutnJets,mc,luminosity,signalsamples.FindSample(nowname));
//    TH1F *mll = signalsamples.Draw("mllhisto","mll",150,0,150, "m_{ll}", "events","",mc,luminosity,signalsamples.FindSample(nowname));
    mll->SetLineColor(TColor::GetColor("#04B404"));
    stringstream poscutS;
    poscutS << "((" << mcjzb <<")>50)";
    TCut poscut(poscutS.str().c_str());
    TH1F *mllP = signalsamples.Draw("mllhistoP","mll",150,0,150, "m_{ll}", "events",cutOSSF&&cutnJets&&poscut,mc,luminosity,signalsamples.FindSample(nowname));
    mllP->SetLineColor(TColor::GetColor("#0040FF"));
    mll->Draw("histo");
    mllP->Draw("histo,same");
    TLegend *leg = make_legend();
    leg->SetY1(0.8);
    leg->AddEntry(mll,(signalsamples.collection)[isig].samplename.c_str(),"L");
    leg->AddEntry(mllP,((signalsamples.collection)[isig].samplename+", JZB>50").c_str(),"L");
    leg->Draw();
    TLine *lin = new TLine(71.2,0,71.2,mll->GetMaximum());
    TLine *lin2 = new TLine(111.2,0,111.2,mll->GetMaximum());
    lin->Draw("same");
    lin2->Draw("same");
    
    CompleteSave(sigcan,"MllShape/"+(signalsamples.collection)[isig].samplename);
    delete mll;
    delete mllP;
  }
}

void met_vs_jzb_plots(string datajzb, string mcjzb) {
  
  TCanvas *canmetjzb = new TCanvas("canmet","MET vs JZB canvas");
  canmetjzb->SetRightMargin(0.16);
  
  vector<string> findme;
  findme.push_back("DY");
  findme.push_back("TTJets");
  findme.push_back("LM");
  /*
  for(int ifind=0;ifind<(int)findme.size();ifind++) {
    vector<int> selsamples = allsamples.FindSample(findme[ifind]);
    TH2F *metvsjzb = new TH2F("metvsjzb","metvsjzb",200,0,100,400,-100,100);
    for(int isel=0;isel<(int)selsamples.size();isel++) {
      dout << "Producing MET:JZB plot ... working on sample: " << allsamples.collection[selsamples[isel]].filename << endl;
      allsamples.collection[selsamples[isel]].events->Draw("jzb[1]:met[4]>>+metvsjzb",cutmass&&cutOSSF&&cutnJets);
    }
    metvsjzb->Scale(allsamples.collection[selsamples[0]].weight);
    metvsjzb->SetStats(0);
    metvsjzb->GetXaxis()->SetTitle("MET (GeV)");
    metvsjzb->GetYaxis()->SetTitle("JZB (GeV)");
    metvsjzb->GetXaxis()->CenterTitle();
    metvsjzb->GetYaxis()->CenterTitle();
    metvsjzb->Draw("COLZ");
    TText* title = write_text(0.5,0.95,allsamples.collection[selsamples[0]].samplename);
    title->SetTextAlign(12);
    title->Draw();
    CompleteSave(canmetjzb,(string)"METvsJZBplots/"+findme[ifind]);
  }
  */
  
  dout << "About to produce MET plot for DY split up by JZB" << endl;
  
  int nbins=14;
  float low=0;
  float high=140;
  
  stringstream sLEFT;
  sLEFT << "((" << mcjzb << ")<0)";
  TCut LEFT(sLEFT.str().c_str());
  stringstream sRIGHT;
  sRIGHT << "((" << mcjzb << ")>0)";
  TCut RIGHT(sRIGHT.str().c_str());
  
  TH1F *metleft   = allsamples.Draw("metleft","met[4]",nbins,low,high, "MET [GeV]", "events", cutmass&&cutOSSF&&cutnJets&&LEFT,mc, luminosity, allsamples.FindSample("DYJets"));
  TH1F *metleftO  = allsamples.Draw("metleftO","met[4]",nbins,low,high, "MET [GeV]", "events", cutmass&&cutOSOF&&cutnJets&&LEFT,mc, luminosity, allsamples.FindSample("DYJets"));
  TH1F *metright  = allsamples.Draw("metright","met[4]",nbins,low,high, "MET [GeV]", "events", cutmass&&cutOSSF&&cutnJets&&RIGHT,mc, luminosity, allsamples.FindSample("DYJets"));
  TH1F *metrightO = allsamples.Draw("metrightO","met[4]",nbins,low,high, "MET [GeV]", "events", cutmass&&cutOSOF&&cutnJets&&RIGHT,mc, luminosity, allsamples.FindSample("DYJets"));
  
  
  TH1F *Bpred =  (TH1F*)metleft->Clone("Bpred");
  Bpred->Add(metleftO,-1);
  Bpred->Add(metrightO);
  TH1F *obs = (TH1F*)metright->Clone("obs");
  
  metleft->Add(metleftO,-1);
  metright->Add(metrightO,-1);
  
  metleft->SetLineColor(kRed);
  metright->SetLineColor(kBlack);
  TPad *metpad = new TPad("metpad","metpad",0,0,1,1);
  metpad->cd();
  metpad->SetLogy(1);
  metleft->Draw("histo");
  metright->Draw("same");
  TLegend *lg = make_legend();
  lg->SetX1(0.5);
  lg->SetY1(0.7);
  lg->AddEntry(metright,"JZB>0 (OSOF corrected)","P");
  lg->AddEntry(metleft,"JZB<0 (OSOF corrected)","L");
  lg->SetHeader("DY");
  
  lg->Draw();
  save_with_ratio(metright,metleft,metpad->cd(),"METvsJZBplots/ComparingLeftToRightinMETspectrum");
  
  TPad *metpad3 = new TPad("metpad3","metpad3",0,0,1,1);
  metpad3->cd();
  metpad3->SetLogy(1);
  Bpred->SetLineColor(kRed);
  Bpred->Draw("histo");
  obs->SetLineColor(kBlack);
  obs->Draw("same");
  TLegend *lg2 = make_legend();
  lg2->SetX1(0.5);
  lg2->SetY1(0.7);
  lg2->AddEntry(obs,"observed","P");
  lg2->AddEntry(Bpred,"predicted","L");
  lg2->SetHeader("DY");

  lg2->Draw();
  
  save_with_ratio(obs,Bpred,metpad3->cd(),"METvsJZBplots/ComparingPredObsinMET");
  
  TPad *metpad2 = new TPad("metpad2","metpad2",0,0,1,1);
  metpad2->cd();
  metpad2->SetLogy(1);
  
  TH1F *metlefta   = allsamples.Draw("metlefta","met[2]",nbins,low,high, "MET [GeV]", "events", cutmass&&cutOSSF&&cutnJets&&LEFT,mc, luminosity, allsamples.FindSample("DYJets"));
  TH1F *metleftOa  = allsamples.Draw("metleftOa","met[2]",nbins,low,high, "MET [GeV]", "events", cutmass&&cutOSOF&&cutnJets&&LEFT,mc, luminosity, allsamples.FindSample("DYJets"));
  TH1F *metrighta  = allsamples.Draw("metrighta","met[2]",nbins,low,high, "MET [GeV]", "events", cutmass&&cutOSSF&&cutnJets&&RIGHT,mc, luminosity, allsamples.FindSample("DYJets"));
  TH1F *metrightOa = allsamples.Draw("metrightOa","met[2]",nbins,low,high, "MET [GeV]", "events", cutmass&&cutOSOF&&cutnJets&&RIGHT,mc, luminosity, allsamples.FindSample("DYJets"));
  
  metlefta->Add(metleftOa,-1);
  metrighta->Add(metrightOa,-1);
  
  metlefta->SetLineColor(kRed);
  metpad2->cd();
  metlefta->Draw("histo");
  metrighta->Draw("same");
  lg->Draw();
  save_with_ratio(metrighta,metlefta,metpad2->cd(),"METvsJZBplots/ComparingLeftToRightinMET_type1_spectrum");
  
  delete Bpred;
  delete obs;
  
  float newhigh=300;
  int newNBins=30;
  
  TPad *metpad4 = new TPad("metpad4","metpad4",0,0,1,1);
  TH1F *Ametleft   = allsamples.Draw("Ametleft","met[4]",newNBins,low,newhigh, "MET [GeV]", "events", cutmass&&cutOSSF&&cutnJets&&LEFT,mc, luminosity);
  TH1F *AmetleftO  = allsamples.Draw("AmetleftO","met[4]",newNBins,low,newhigh, "MET [GeV]", "events", cutmass&&cutOSOF&&cutnJets&&LEFT,mc, luminosity);
  TH1F *Ametright  = allsamples.Draw("Ametright","met[4]",newNBins,low,newhigh, "MET [GeV]", "events", cutmass&&cutOSSF&&cutnJets&&RIGHT,mc, luminosity);
  TH1F *AmetrightO = allsamples.Draw("AmetrightO","met[4]",newNBins,low,newhigh, "MET [GeV]", "events", cutmass&&cutOSOF&&cutnJets&&RIGHT,mc, luminosity);
  
  TH1F *aBpred = (TH1F*)Ametleft->Clone("aBpred");
  aBpred->Add(AmetleftO,-1);
  aBpred->Add(AmetrightO);
  aBpred->SetLineColor(kRed);
  
  TH1F *aobs = (TH1F*)Ametright->Clone("aobs");
  metpad4->cd();
  metpad4->SetLogy(1);
  aobs->Draw();
  aBpred->Draw("histo,same");
  aobs->Draw("same");
  lg->SetHeader("All MC");
  lg->Draw();
  save_with_ratio(aobs,aBpred,metpad4->cd(),"METvsJZBplots/ComparingPredObsinMET_ALLSAMPLES");
  
  
  delete lg;
  delete canmetjzb;
  delete metleft;
  delete metleftO;
  delete metright;
  delete metrightO;
}
    

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");
  dout << "HELLO there!" << endl;
  
}
Revision:	1.48
Committed:	Wed Aug 15 09:18:07 2012 UTC (12 years, 8 months ago) by buchmann
Content type:	text/plain
Branch:	MAIN
Changes since 1.47:	+97 -21 lines
Log Message:	Added some kinematic plots as requested, and ZJets estimate study