Plotting/Modules/MetPlotting.C

#include <iostream>

using namespace std;

float Get_Met_Z_Prediction(TCut JetCut, float MetCut, int isdata, bool isDYonly, bool isAachen);

namespace MetPlotsSpace {
  float Zprediction_Uncertainty=0.2;
  float OFprediction_Uncertainty=0.07;
  float Zestimate__data=-1;
  float Zestimate__data_sys=-1;
  float Zestimate__data_stat=-1;
  float Zestimate__mc=-1;
  float Zestimate__mc_sys=-1;
  float Zestimate__mc_stat=-1;
  float Zestimate__dy=-1;
  float Zestimate__dy_sys=-1;
  float Zestimate__dy_stat=-1;
}

void ExperimentalMetPrediction(bool QuickRun, bool isAachen);
void ProvideEEOverMMEstimate(TCut); 

void makeOneRinoutPlot( TH2F* hrange, Int_t* bins, Int_t nBins, TString var, string name, bool doMC, TCut kCut = "" ) {

  Float_t systematics = 0.25;

  // mll settings
  Int_t nbins = 100;
  Float_t xmin = 20., xmax = 120.;
  TCanvas* mycan = new TCanvas("mycan","Canvas");
  mycan->SetLeftMargin(0.2);
  mycan->SetLogy(0);
  

  TCut kbase("pfJetGoodNum40>1&&pfJetGoodID[0]!=0"&&passtrig&&kCut);
  TCut kSF("id1==id2");
  TCut kOF("id1!=id2");

  // Reference: inclusive selection
  TCut kZP("pfJetGoodNum40==2");
  TH1F* h1, *h1OF;
  if ( !doMC ) {
    h1   = allsamples.Draw("h1",  "mll",nbins,xmin,xmax,"m_{ll}","events",kbase&&kZP&&kSF,data,luminosity);
    h1OF = allsamples.Draw("h1OF","mll",nbins,xmin,xmax,"m_{ll}","events",kbase&&kZP&&kOF,data,luminosity);
  } else {
    h1   = allsamples.Draw("h1",  "mll",nbins,xmin,xmax,"m_{ll}","events",kbase&&kZP&&kSF,mc,luminosity,allsamples.FindSample("Z_em"));
    h1OF = allsamples.Draw("h1OF","mll",nbins,xmin,xmax,"m_{ll}","events",kbase&&kZP&&kOF,mc,luminosity,allsamples.FindSample("Z_em"));
  }

  Int_t minBinSR = h1->FindBin(20.);
  Int_t maxBinSR = h1->FindBin(70.)-1;

  Int_t minBinZP = h1->FindBin(81.);
  Int_t maxBinZP = h1->FindBin(101.)-1;

  dout << "Integrating SR from " << h1->GetBinLowEdge(minBinSR) << " to " << h1->GetBinLowEdge(maxBinSR)+h1->GetBinWidth(maxBinSR) << std::endl;
  dout << "Integrating ZP from " << h1->GetBinLowEdge(minBinZP) << " to " << h1->GetBinLowEdge(maxBinZP)+h1->GetBinWidth(maxBinZP) << std::endl;

  // Subtract OF
  h1->Add(h1OF,-1);
  h1->SetLineColor(kRed);

  // Compute ratio
  Double_t yZP, eyZP, ySR, eySR;
  ySR = h1->IntegralAndError(minBinSR,maxBinSR,eySR);
  yZP = h1->IntegralAndError(minBinZP,maxBinZP,eyZP);
  dout << "Ratio: " << ySR/yZP << "+-" << computeRatioError(ySR,eySR,yZP,eyZP) << std::endl;

  std::stringstream twoJetsLegend;  
  twoJetsLegend << std::setprecision(1) << std::fixed;
  twoJetsLegend << "2-jets ratio: (" << ySR/yZP*100. << "#pm" << computeRatioError(ySR,eySR,yZP,eyZP)*100. << "#pm" << systematics*ySR/yZP*100. << ")%";

  TLine* line = new TLine(hrange->GetXaxis()->GetXmin(),ySR/yZP,hrange->GetXaxis()->GetXmax(),ySR/yZP);
  line->SetLineColor(kRed);
  TBox* errorBox = new TBox(hrange->GetXaxis()->GetXmin(),ySR/yZP*(1-systematics),hrange->GetXaxis()->GetXmax(),ySR/yZP*(1+systematics));
  errorBox->SetFillColor(kCyan);
  errorBox->SetFillStyle(1001);
  errorBox->SetLineColor(kWhite);

  TGraphErrors* ratio = new TGraphErrors(nBins);
  // Various cuts
  for ( int ibin = 0; ibin<nBins; ++ibin ) {
    std::stringstream cut;
    cut << var << ">=" << bins[ibin];
    if ( ibin+1<nBins ) cut << "&&" << var << "<" << bins[ibin+1];
    TCut kadd(cut.str().c_str());

    TH1F* h2, *h2OF;
    if ( !doMC ) {
      h2   = allsamples.Draw("h2",  "mll",nbins,xmin,xmax,"var","events",kbase&&kadd&&kSF,data,luminosity);
      h2OF = allsamples.Draw("h2OF","mll",nbins,xmin,xmax,"var","events",kbase&&kadd&&kOF,data,luminosity);
    } else {
      h2   = allsamples.Draw("h2",  "mll",nbins,xmin,xmax,"var","events",kbase&&kadd&&kSF,mc,luminosity,allsamples.FindSample("Z_em"));
      h2OF = allsamples.Draw("h2OF","mll",nbins,xmin,xmax,"var","events",kbase&&kadd&&kOF,mc,luminosity,allsamples.FindSample("Z_em"));
    }
    h2->Add(h2OF,-1);
    h2->SetLineColor(kBlue);

    ySR = h2->IntegralAndError(minBinSR,maxBinSR,eySR);
    yZP = h2->IntegralAndError(minBinZP,maxBinZP,eyZP);

    if ( ibin+1<nBins ) {
      ratio->SetPoint(ibin,(bins[ibin+1]+bins[ibin])/2.0,ySR/yZP);
      ratio->SetPointError(ibin,(bins[ibin+1]-bins[ibin])/2.0,computeRatioError(ySR,eySR,yZP,eyZP));    
    } else {
      Float_t width = ratio->GetErrorX(ibin-1);
      ratio->SetPoint(ibin,bins[ibin]+width,ySR/yZP);
      ratio->SetPointError(ibin,width,computeRatioError(ySR,eySR,yZP,eyZP));
    }
    
    dout << "Ratio " << cut.str() << ": " << ySR/yZP << "+-" << computeRatioError(ySR,eySR,yZP,eyZP) << std::endl;

    h2->Delete();
    h2OF->Delete();

  }

  std::stringstream syserrLegend;
  syserrLegend << std::setprecision(0) << std::fixed << systematics*100. << "% systematic unc.";

  hrange->GetYaxis()->SetTitleOffset(1.3);
  hrange->GetYaxis()->SetDecimals(kTRUE);
  hrange->Draw();
  errorBox->Draw();
  line->Draw();
  ratio->Draw("P");

  TLegend* legend = new TLegend(0.25,0.6,0.8,0.9);
  legend->SetFillStyle(0);
  legend->SetBorderSize(0);
  if ( doMC ) legend->AddEntry(ratio,"DY Z+jets MC","lp");
  else legend->AddEntry(ratio,"Data","lp");
  legend->AddEntry(line,twoJetsLegend.str().c_str(),"l");
  legend->AddEntry(errorBox,syserrLegend.str().c_str(),"f");
  legend->Draw();

  mycan->RedrawAxis();
  if (!doMC) DrawPrelim();
  else DrawMCPrelim();

  CompleteSave(mycan,"MetPlots/Zlineshape_vs_"+name+(doMC?"_mc_":""));

  h1->Delete();
  h1OF->Delete();
  delete mycan;

}

int zlineshapeMet(bool doMC=true, string suffix="", float ymax = 0.2, TCut kCut = "" ) {

  TH2F* hrange = new TH2F("hrange","Range ; MET [GeV] ; Ratio low mass / Z peak",2,-1,61,2,0,ymax);
  Int_t metBins[] = { 0, 10, 20, 30, 40, 50 };
  Int_t nMetBins = sizeof(metBins)/sizeof(Int_t);
  makeOneRinoutPlot( hrange, metBins, nMetBins, "met[4]/MetFactor", "met"+suffix, doMC, kCut );
  hrange->Delete();
  return 0;

}

int zlineshapeJets(bool doMC=true, string suffix="", float ymax = 0.2, TCut kCut = "" ) {
  TH2F* hrange = new TH2F("hrange","Range ; #(jets) ; Ratio low mass / Z peak",2,1.9,6.1,2,0,ymax);
  Int_t bins[] = { 2, 3, 4, 5 };

  Int_t nBins = sizeof(bins)/sizeof(Int_t);
  makeOneRinoutPlot( hrange, bins, nBins, "pfJetGoodNum40", "njets"+suffix, doMC, kCut );
  hrange->Delete();
  return 0;

}

int zlineshapes(string suffix = "", TCut cut="" ) {

  dout << "--- Calculating R_in/out" << std::endl;
  zlineshapeMet(false,suffix,0.2,cut);
  zlineshapeJets(false,suffix,0.2,cut);
  zlineshapeMet(true,suffix,0.2,cut);
  zlineshapeJets(true,suffix,0.2,cut);
  dout << "--- DONE (Calculating R_in/out)" << std::endl;

  return 0;
}

void ExtractScaleFactor(TH1F *mllSF,TH1F *mllOF, THStack* mcMllSF, THStack* mcMllOF, TH1F *prediction, TLegend *leg, string saveasSig, TBox *srbox) {
  Int_t minbin = mllSF->FindBin(20.);
  Int_t maxbin = mllSF->FindBin(70.-1);

  // Get yields in OF region
  Float_t iDataOF = mllOF->Integral();
  Float_t iDataOFSR = mllOF->Integral(minbin,maxbin);
  Float_t iMCOF = 0.0;
  Float_t iMCOFSR = 0.0;
  TIter nextOF(mcMllOF->GetHists());
  TH1F* h;
  while ( h = (TH1F*)nextOF() ) {
    iMCOF += h->Integral();
    iMCOFSR += h->Integral(minbin,maxbin);
  }
  Float_t scale = iDataOF/iMCOF;
 
  // Re-scale OF
  nextOF = TIter(mcMllOF->GetHists());
  
  while ( h = (TH1F*)nextOF() ) {
    h->Scale(scale);
  }
  
  nextOF = TIter(mcMllOF->GetHists());
  
  // Rescale SF and count in signal region
  dout << "Integrating from " << mllSF->GetBinLowEdge(minbin) << " to " << mllSF->GetBinLowEdge(maxbin)+mllSF->GetBinWidth(maxbin) << std::endl;

  Float_t iDataSFSR = mllSF->Integral(minbin,maxbin);
  Float_t iMCSFSR = 0.0;
  TIter nextSF = TIter(mcMllSF->GetHists());
  while ( h = (TH1F*)nextSF() ) {
    h->Scale(scale);
    iMCSFSR += h->Integral(minbin,maxbin);
  }
  
  nextSF = TIter(mcMllSF->GetHists());
  while ( h = (TH1F*)nextSF() ) {
    iMCSFSR += h->Integral(minbin,maxbin);
  }
  mcMllSF->Modified();
  
  TPad* rcan2 = new TPad("rcan2","rcan2",0,0,1,1);
  rcan2->cd();
  mllSF->Draw();
  mcMllSF->Draw("histo,same");
  prediction->Draw("histo,same");
  mllSF->Draw("same");
  DrawPrelim();
  stringstream leghead;
  leghead << "MC scaled by " << std::setprecision(2) << scale << "";
  dout << "SCALE: " << scale << endl;
  TH1F *histo = new TH1F("histo","histo",1,0,1);histo->SetLineColor(kWhite);
  leg->AddEntry(histo,leghead.str().c_str(),"l");
  leg->Draw();
  srbox->Draw();
  stringstream saveasSig2;
  saveasSig2 << saveasSig << "__mcScaled";
  rcan2->Update();
  Save_With_Ratio( mllSF, *mcMllSF, rcan2, saveasSig2.str() );
  
  // restore original stacks
  nextOF = TIter(mcMllOF->GetHists());
  
  while ( h = (TH1F*)nextOF() ) {
    h->Scale(1/scale);
  }
  
  nextSF = TIter(mcMllSF->GetHists());
  while ( h = (TH1F*)nextSF() ) {
    h->Scale(1/scale);
  }
  mcMllSF->Modified();
  mcMllOF->Modified();
  
}

TGraphErrors* MakeErrorGraph(TH1F *histo) { 
  
  float dx[histo->GetNbinsX()];
  float dy[histo->GetNbinsX()];
  float x[histo->GetNbinsX()];
  float y[histo->GetNbinsX()];
  for(int i=1;i<=histo->GetNbinsX();i++) {
    x[i-1]=histo->GetBinCenter(i);
    y[i-1]=histo->GetBinContent(i);
    if(i>1) dx[i-1]=(histo->GetBinCenter(i)-histo->GetBinCenter(i-1))/2.0;
    else dx[i-1]=(histo->GetBinCenter(i+1)-histo->GetBinCenter(i))/2.0;
    dy[i-1]=histo->GetBinError(i);
  }
  
  TGraphErrors *gr = new TGraphErrors(histo->GetNbinsX(),x,y,dx,dy);
  gr->SetFillColor(TColor::GetColor("#2E9AFE"));
  return gr;
}

TGraphErrors* MakeErrorGraphSystematicAndStatistical(TH1F *ofpred, TH1F *sfpred, TH1F *prediction, TH1F *SystHisto) { 
  
  float dx[ofpred->GetNbinsX()];
  float dy[ofpred->GetNbinsX()];
  float x[ofpred->GetNbinsX()];
  float y[ofpred->GetNbinsX()];
  for(int i=1;i<=ofpred->GetNbinsX();i++) {
    x[i-1]=prediction->GetBinCenter(i);
    y[i-1]=prediction->GetBinContent(i);
    if(i>1) dx[i-1]=(prediction->GetBinCenter(i)-prediction->GetBinCenter(i-1))/2.0;
    else dx[i-1]=(prediction->GetBinCenter(i+1)-prediction->GetBinCenter(i))/2.0;
    if(ofpred->GetBinCenter(i)>20 && ofpred->GetBinCenter(i)<70) {
      //need to increase uncertainty by 5% due to extrapolation
      dy[i-1] = (MetPlotsSpace::Zprediction_Uncertainty+0.05)*(MetPlotsSpace::Zprediction_Uncertainty+0.05)*sfpred->GetBinContent(i)*sfpred->GetBinContent(i); //systematic for Z+Jets prediction
    } else {
      dy[i-1] = MetPlotsSpace::Zprediction_Uncertainty*MetPlotsSpace::Zprediction_Uncertainty*sfpred->GetBinContent(i)*sfpred->GetBinContent(i); //systematic for Z+Jets prediction
    }
    dy[i-1]+= MetPlotsSpace::OFprediction_Uncertainty*MetPlotsSpace::OFprediction_Uncertainty* ofpred->GetBinContent(i) * ofpred->GetBinContent(i); //systematic for OF prediction
    float sys=sqrt(dy[i-1])/prediction->GetBinContent(i);
    if(prediction->GetBinContent(i)==0) sys=0.0;
    if(sys!=sys || sys<0) sys=0;
    SystHisto->SetBinContent(i,sys);
    dy[i-1]+= prediction->GetBinError(i) * prediction->GetBinError(i); // plus statistical!
    dy[i-1]=sqrt(dy[i-1]);
  }
  
  TGraphErrors *gr = new TGraphErrors(ofpred->GetNbinsX(),x,y,dx,dy);
  gr->SetFillColor(TColor::GetColor("#2E9AFE"));//blue
//  gr->SetFillColor(TColor::GetColor("#FF8000"));//orange
  return gr;
}

float GetYield(TH1F *histo, float min, float max) {
  float res=0.0;
  for(int i=1;i<=histo->GetNbinsX();i++) {
      if(histo->GetBinLowEdge(i)+histo->GetBinWidth(i)<=min) continue;// considered in the next bin
      if(histo->GetBinLowEdge(i)>=max) continue;//above the threshold, out!
//      cout << "         added yield for bin [" << histo->GetBinLowEdge(i) << "," << histo->GetBinLowEdge(i)+histo->GetBinWidth(i) << "] which is " << histo->GetBinContent(i) << "   in GetYield" << endl;
      res+=histo->GetBinContent(i);
  }
  return res;
}

void ProduceYields(float min, float max, TH1F *data, THStack *stack) {
  dout << "   *************** <MC YIELDS> ********* " << endl;
  dout << "   Considering " << min << " < mll < " << max << "           " << endl;
  dout << "   Data : " << GetYield(data,min,max)  << endl;
  TIter nextSF(stack->GetHists());
  TH1F* h;
  while ( h = (TH1F*)nextSF() ) {
    dout << "   " << h->GetName() << " : " << GetYield(h,min,max)  << endl;
  }
  dout << "   *************** </MC YIELDS> ********* " << endl;
}
 
  
void WriteYield(THStack *mc, float low, float high) {
  TH1F *h;
  TIter NextHisto(mc->GetHists());
  while ( h = (TH1F*)NextHisto() ) {
    float CurrBKG=0;
    float CurrBKGErr=0;
//    cout << "        DEBUG : now looking at histo " << h->GetName() << endl;
    for(int i=1;i<=h->GetNbinsX()+1;i++) {
      if(h->GetBinLowEdge(i)+h->GetBinWidth(i)<=low) continue;// considered in the next bin
      if(h->GetBinLowEdge(i)>=high) continue;//above the threshold, out!
      CurrBKG+=h->GetBinContent(i);
      CurrBKGErr=sqrt(CurrBKGErr*CurrBKGErr+h->GetBinError(i)*h->GetBinError(i));
//      cout << "        DEBUG :: Considering bin " << h->GetBinLowEdge(i) << "," << h->GetBinLowEdge(i)+h->GetBinWidth(i) << " for  mll in [" << low << "," << high << "] ; contains " << h->GetBinContent(i) << endl;
    }
    dout << "      " << h->GetName() << " : " << CurrBKG << " +/- " << CurrBKGErr << endl;
  }
}

  
void ProduceMetPlotsWithCut(bool isAachen, TCut cut, string name, float cutat, int njets, bool doMC = false, float ymax = 80 ) {
  
  dout << " *-*-*-*-*-*-*-*- " << __FUNCTION__ << " *-*-*-*-*-*-*-*- " << endl;
  
  ofstream NiceSummary;
  NiceSummary.open("Summary.txt",ios::app);
  
  NiceSummary << "  *********************************" << endl;
  time_t t = time(0);   // get time now
  struct tm * now = localtime( & t );
  NiceSummary << "run on : " << now->tm_hour << ":" << now->tm_min << " on " << now->tm_mday << "." << (now->tm_mon + 1) << "." << (now->tm_year + 1900) << endl;
  NiceSummary << "Cut : " << (const char*) cut << endl;
  NiceSummary << "Name : " << name << endl;
  bool UseSpecialZprediction=false;
  
  TText *sel = WriteSelection(njets);
  if(cutat==100 && name=="") {
    UseSpecialZprediction=true;
    bool ReRunEstimate=false;
    //need to check if the results have already been stored; if not, need to get the estimate!
    if(MetPlotsSpace::Zestimate__data<0) ReRunEstimate=true;
    if(MetPlotsSpace::Zestimate__data_stat<0) ReRunEstimate=true;
    if(MetPlotsSpace::Zestimate__data_sys<0) ReRunEstimate=true;
    if(MetPlotsSpace::Zestimate__mc<0) ReRunEstimate=true;
    if(MetPlotsSpace::Zestimate__mc_stat<0) ReRunEstimate=true;
    if(MetPlotsSpace::Zestimate__mc_sys<0) ReRunEstimate=true;
    if(MetPlotsSpace::Zestimate__dy<0) ReRunEstimate=true;
    if(MetPlotsSpace::Zestimate__dy_stat<0) ReRunEstimate=true;
    if(MetPlotsSpace::Zestimate__dy_sys<0) ReRunEstimate=true;
    dout << "****************** About to do Z prediction " << endl;
    if(ReRunEstimate) ExperimentalMetPrediction(true,isAachen);//doing quick run (i.e. only data)
    dout << "****************** Done predicting the Z " << endl;
  }


  TCanvas *tcan = new TCanvas("tcan","tcan");
  dout << "Doing met plots" << endl;
  stringstream MetBaseCuts;
  MetBaseCuts << "met[4]/MetFactor>" << cutat << "&&" << cut.GetTitle() << " && mll>20";
  stringstream snjets;
  snjets << njets;
  TCut MetBaseCut(MetBaseCuts.str().c_str());
  TCut nJetsSignal(PlottingSetup::basicqualitycut&&TCut(("pfJetGoodNum40>="+snjets.str()).c_str())&&TCut("mll>20"));
  TCut nJetsControl(PlottingSetup::basiccut&&"met[4]/MetFactor>100&&met[4]/MetFactor<150&&pfJetGoodID[0]!=0&&pfJetGoodNum40==2&&mll>20"); // Common CR (modulo lepton selection)
  //TCut nJetsControl(PlottingSetup::basiccut&&"met[4]/MetFactor>75&&met[4]/MetFactor<150&&pfJetGoodNumBtag30>0&&pfJetGoodID[0]!=0&&pfJetGoodNum40==2"); // Alternative CR

  //compute SF / OF rate in (CR1+CR2), should give 0.941 +/- 0.05

  // Create histograms
  //int nbins = 30;
  int nbins = 60-4;
  float xmin=20., xmax = 300.;
  
  TH1F *mllsigEE = allsamples.Draw("mllsigEE","mll",nbins,xmin,xmax,"m_{ee} [GeV]",    "events",cutOSSF&&MetBaseCut&&nJetsSignal&&TCut("id1==0"),data,PlottingSetup::luminosity);
  TH1F *mllsigMM = allsamples.Draw("mllsigMM","mll",nbins,xmin,xmax,"m_{#mu#mu} [GeV]","events",cutOSSF&&MetBaseCut&&nJetsSignal&&TCut("id1==1"),data,PlottingSetup::luminosity);
  TH1F *mllsconEE  = allsamples.Draw("mllsconEE","mll",nbins,xmin,xmax,"m_{ll} [GeV]",    "events",cutOSSF&&cut&&nJetsControl&&TCut("id1==0"),data,PlottingSetup::luminosity);
  TH1F *mllsconMM  = allsamples.Draw("mllsconMM","mll",nbins,xmin,xmax,"m_{ll} [GeV]",    "events",cutOSSF&&cut&&nJetsControl&&TCut("id1==1"),data,PlottingSetup::luminosity);
  TH1F *mllOsigEM  = allsamples.Draw("mllOsigEM", "mll",nbins,xmin,xmax,"m_{ll} [GeV]","events",cutOSOF&&MetBaseCut&&nJetsSignal&&TCut("id1==0"),data,PlottingSetup::luminosity);
  TH1F *mllOsigME  = allsamples.Draw("mllOsigME", "mll",nbins,xmin,xmax,"m_{ll} [GeV]","events",cutOSOF&&MetBaseCut&&nJetsSignal&&TCut("id1==1"),data,PlottingSetup::luminosity);
  TH1F *mllOsconEM = allsamples.Draw("mllOsconEM","mll",nbins,xmin,xmax,"m_{ll} [GeV]","events",TCut(cutOSOF&&cut&&nJetsControl&&TCut("id1==0")),data,PlottingSetup::luminosity);
  TH1F *mllOsconME = allsamples.Draw("mllOsconME","mll",nbins,xmin,xmax,"m_{ll} [GeV]","events",TCut(cutOSOF&&cut&&nJetsControl&&TCut("id1==1")),data,PlottingSetup::luminosity);
  TH1F *ptsig    = allsamples.Draw("ptsig",  "pt",40,xmin,400,"m_{T}^{ll} [GeV]","events",TCut(cutOSSF&&MetBaseCut&&nJetsSignal),data,PlottingSetup::luminosity);
  TH1F *ptOsig   = allsamples.Draw("ptOsig", "pt",40,xmin,400,"p_{T}^{ll} [GeV]","events",TCut(cutOSOF&&MetBaseCut&&nJetsSignal),data,PlottingSetup::luminosity);
  TH1F *mllscon = (TH1F*)mllsconEE->Clone("mllscon");
  mllscon->Add(mllsconMM);
  TH1F *mllOscon = (TH1F*)mllOsconEM->Clone("mllOscon");
  mllOscon->Add(mllOsconME);
  TH1F *mllOsig  = (TH1F*)mllOsigEM->Clone("mllOsig");
  mllOsig->Add(mllOsigME);
  
  TH1F* mllsig = (TH1F*)mllsigEE->Clone("mllsig");
  mllsig->Add(mllsigMM);
  mllsig->GetXaxis()->SetTitle("m_{ll} [GeV]");

  THStack *mcMllsig, *mcMllsigEE,*mcMllsigMM,*mcMllscon,*mcMllsconEE,*mcMllsconMM, *mcMllOsig, *mcMllOscon, *mcptsig, *mcptOsig;
  if ( doMC ) {
    name += "_mc";
    mcMllsig  = new THStack(allsamples.DrawStack("mcMllsig","mll",nbins,xmin,xmax,"m_{ll} [GeV]","events",TCut(cutOSSF&&MetBaseCut&&nJetsSignal),mc,PlottingSetup::luminosity));
    mcMllsigEE  = new THStack(allsamples.DrawStack("mcMllsigEE","mll",nbins,xmin,xmax,"m_{ee} [GeV]","events",TCut(cutOSSF&&MetBaseCut&&nJetsSignal&&"id1==0"),mc,PlottingSetup::luminosity));
    mcMllsigMM  = new THStack(allsamples.DrawStack("mcMllsigMM","mll",nbins,xmin,xmax,"m_{#mu#mu} [GeV]","events",TCut(cutOSSF&&MetBaseCut&&nJetsSignal&&"id1==1"),mc,PlottingSetup::luminosity));
    mcMllscon = new THStack(allsamples.DrawStack("mcMllscon","mll",nbins,xmin,xmax,"m_{ll} [GeV]","events",TCut(cutOSSF&&cut&&nJetsControl),mc,PlottingSetup::luminosity));
    mcMllOsig = new THStack(allsamples.DrawStack("mcMllOsig","mll",nbins,xmin,xmax,"m_{ll} [GeV]","events",TCut(cutOSOF&&MetBaseCut&&nJetsSignal),mc,PlottingSetup::luminosity));
    mcMllOscon= new THStack(allsamples.DrawStack("mcMllOscon","mll",nbins,xmin,xmax,"m_{ll} [GeV]","events",TCut(cutOSOF&&cut&&nJetsControl),mc,PlottingSetup::luminosity));
    mcptsig    = new THStack(allsamples.DrawStack("mcptsig",  "pt",40,xmin,400,"m_{T}^{ll} [GeV]","events",TCut(cutOSSF&&MetBaseCut&&nJetsSignal),mc,PlottingSetup::luminosity));
    mcptOsig   = new THStack(allsamples.DrawStack("mcptOsig", "pt",40,xmin,400,"p_{T}^{ll} [GeV]","events",TCut(cutOSOF&&MetBaseCut&&nJetsSignal),mc,PlottingSetup::luminosity));
  }
  
  
  if ( doMC ) {
    write_info(__FUNCTION__,"Writing yields in MC");
    dout << "Signal region, low mass " << endl;
    WriteYield(mcMllsig,20,70);
    WriteYield(mcMllOsig,20,70);
    dout << "Signal region, high mass " << endl;
    WriteYield(mcMllsig,120,1000);
    WriteYield(mcMllOsig,120,1000);
    THStack *sdata = new THStack();
    sdata->Add(mllsig);
    THStack *sodata = new THStack();
    sodata->Add(mllOsig);
    dout << "NOW ALL DATA **** " << endl;
    dout << "Signal region, low mass " << endl;
    dout << "SF : " << endl;
    WriteYield(sdata,20,70);
    dout << endl;
    dout << "OF : " << endl;
    WriteYield(sodata,20,70);
    dout << "Signal region, high mass " << endl;
    dout << "SF : " << endl;
    WriteYield(sdata,120,1000);
    dout << endl;
    dout << "OF : " << endl;
    WriteYield(sodata,120,1000);
  }
  
  if(doMC) {
    TH1F *mllsigt = allsamples.Draw("mllsigt","mll",300,0,300,"m_{#mu#mu} [GeV]","events",TCut(cutOSSF&&MetBaseCut&&nJetsSignal),data,PlottingSetup::luminosity);
    THStack *mcMllsigt  = new THStack(allsamples.DrawStack("mcMllsigt","mll",300,0,300,"m_{ll} [GeV]","events",TCut(cutOSSF&&MetBaseCut&&nJetsSignal),mc,PlottingSetup::luminosity));
    ProduceYields(20,70,mllsigt,mcMllsigt);
    ProduceYields(91-10,91+10,mllsigt,mcMllsigt);
    write_warning(__FUNCTION__,"WATCH OUT Z REGION HAS BEEN DEACTIVATED");
    cout << (const char*) TCut(cutOSSF&&MetBaseCut&&nJetsSignal&&"id1==0") << endl;
    delete mllsigt;
    delete mcMllsigt;
  }
    
  mllOsig->SetLineColor(kRed);
  mllOscon->SetLineColor(kRed);
  
  TH1F *zlineshape = allsamples.Draw("zlineshape","mll",nbins,xmin,xmax,"m_{ll} (GeV)","events",cutOSSF&&TCut("pfJetGoodNum40==2")&&cut,data,PlottingSetup::luminosity);
  TH1F *Ozlineshape = allsamples.Draw("Ozlineshape","mll",nbins,xmin,xmax,"m_{ll} (GeV)","events",cutOSOF&&TCut("pfJetGoodNum40==2")&&cut,data,PlottingSetup::luminosity);
  zlineshape->Add(Ozlineshape,-1);
  TH1F *zlineshapeControl = (TH1F*)zlineshape->Clone("zlineshapeControl");
//   TH1F *zlineshapeFINE = allsamples.Draw("zlineshapeFINE","mll",50*nbins,xmin,xmax,"m_{ll} (GeV)","events",cutOSSF&&TCut("pfJetGoodNum40==1")&&cut,data,PlottingSetup::luminosity);
//   
//   float scalefactor = Get_Met_Z_Prediction(nJetsSignal,cutat, data, false) / (zlineshapeFINE->Integral(zlineshapeFINE->FindBin(91.1-20),zlineshapeFINE->FindBin(91.1+20)));
//   float scalefactor_Control = Get_Met_Z_Prediction(nJetsControl,cutat, data, false) / (zlineshapeFINE->Integral(zlineshapeFINE->FindBin(91.1-20),zlineshapeFINE->FindBin(91.1+20)));
//   delete zlineshapeFINE;
  
  
  Int_t scaleBinLow = mllsig->FindBin(86);
  Int_t scaleBinHigh = mllsig->FindBin(94);
  float scalefactor = (mllsig->Integral(scaleBinLow,scaleBinHigh)-mllOsig->Integral(scaleBinLow,scaleBinHigh));
  scalefactor /= zlineshape->Integral(scaleBinLow,scaleBinHigh);
  
  float scalefactor_Control = (mllscon->Integral(scaleBinLow,scaleBinHigh)-mllOscon->Integral(scaleBinLow,scaleBinHigh));
  scalefactor_Control /= zlineshapeControl->Integral(scaleBinLow,scaleBinHigh);
    
  dout << "Bins for scaling : " << scaleBinLow << " : " << scaleBinHigh << endl;
  
  if(UseSpecialZprediction) {
    scaleBinLow = mllsig->FindBin(81);
    scaleBinHigh = mllsig->FindBin(101);
    scalefactor = MetPlotsSpace::Zestimate__data/ (zlineshape->Integral(scaleBinLow,scaleBinHigh));
    dout << "Dividing: " << MetPlotsSpace::Zestimate__data << " by " <<  (zlineshape->Integral(scaleBinLow,scaleBinHigh)) << endl;
    write_warning(__FUNCTION__,"Not using JZB prediction for control region!");
  }
  
  dout << "Scale factors    : " << scalefactor << " : " << scalefactor_Control << endl;
  if(UseSpecialZprediction) dout << " NOTE: Used JZB prediction for scaling! (Bins )" << scaleBinLow << " to " << scaleBinHigh << endl;
  
  zlineshape->Scale(scalefactor);
    
  zlineshape->SetLineColor(kBlue);
  zlineshape->SetLineStyle(2);
  
  zlineshapeControl->Scale(scalefactor_Control);
  zlineshapeControl->SetLineColor(kBlue);
  zlineshapeControl->SetLineStyle(2);
  
  
  if(UseSpecialZprediction) {
    //need to update each bin with correct stat uncert
    float relDYerr = MetPlotsSpace::Zestimate__data_stat/MetPlotsSpace::Zestimate__data;
    for(int iz=1;iz<=zlineshape->GetNbinsX();iz++) {
      float bincontent=zlineshape->GetBinContent(iz);
      float binerror=zlineshape->GetBinError(iz);
      float finalerr=0;
      if(bincontent>0) finalerr+= (binerror/bincontent) * (binerror/bincontent);
      if(MetPlotsSpace::Zestimate__data>0) finalerr+= relDYerr*relDYerr;
      finalerr=bincontent * TMath::Sqrt(finalerr);
      zlineshape->SetBinError(iz,finalerr);
    }
  }
  

  //**************************
  write_info(__FUNCTION__,"LOW MASS YIELDS : ");
  float yields_mllsf=0,yields_mllmm=0,yields_mllee=0,yields_mllof=0,yields_mllofme=0,yields_mllofem=0;
  for(int i=1;i<=mllsigEE->GetNbinsX();i++) {
    if(mllsigEE->GetBinCenter(i)>20 && mllsigEE->GetBinCenter(i)<70) {
      yields_mllee+=mllsigEE->GetBinContent(i);
      yields_mllmm+=mllsigMM->GetBinContent(i);
      yields_mllsf+=mllsigEE->GetBinContent(i);
      yields_mllsf+=mllsigMM->GetBinContent(i);
      yields_mllof+=mllOsig->GetBinContent(i);
      yields_mllofem+=mllOsigEM->GetBinContent(i);
      yields_mllofme+=mllOsigME->GetBinContent(i);
    }
  }
  dout << "Observed : " << yields_mllsf << " (" << yields_mllee << "ee , " << yields_mllmm << "mm )" << endl;
  dout << "Predicted: " << yields_mllof << " (from OF, " << yields_mllofem << " em, " << yields_mllofme << " me) " << endl;
  
  NiceSummary << "Low mass observed : " << yields_mllsf << " (" << yields_mllee << "ee , " << yields_mllmm << "mm )" << endl;
  NiceSummary << "Low mass predicted: " << yields_mllof << " (from OF, " << yields_mllofem << " em, " << yields_mllofme << " me) " << endl;

  
  //**************************
  write_info(__FUNCTION__,"CONTROL REGION YIELDS : ");
  float cr1yields_mllsf=0,cr1yields_mllmm=0,cr1yields_mllee=0,cr1yields_mllof=0,cr1yields_mllofem=0,cr1yields_mllofme=0;
  float cr2yields_mllsf=0,cr2yields_mllmm=0,cr2yields_mllee=0,cr2yields_mllof=0,cr2yields_mllofem=0,cr2yields_mllofme=0;
  float pyields_mllsf=0,pyields_mllmm=0,pyields_mllee=0,pyields_mllof=0,pyields_mllofem=0,pyields_mllofme=0,pyields_mlldy=0,pyields_mllpred;
  for(int i=1;i<=mllsigEE->GetNbinsX();i++) {
    if(mllsigEE->GetBinCenter(i)>20 && mllsigEE->GetBinCenter(i)<70) {
      cr1yields_mllee+=mllsconEE->GetBinContent(i);
      cr1yields_mllmm+=mllsconMM->GetBinContent(i);
      cr1yields_mllsf+=mllsconEE->GetBinContent(i);
      cr1yields_mllsf+=mllsconMM->GetBinContent(i);
      cr1yields_mllof+=mllOscon->GetBinContent(i);
      cr1yields_mllofem+=mllOsconEM->GetBinContent(i);
      cr1yields_mllofme+=mllOsconME->GetBinContent(i);
    }
    if(mllsigEE->GetBinCenter(i)>120) {
      cr2yields_mllee+=mllsconEE->GetBinContent(i);
      cr2yields_mllmm+=mllsconMM->GetBinContent(i);
      cr2yields_mllsf+=mllsconEE->GetBinContent(i);
      cr2yields_mllsf+=mllsconMM->GetBinContent(i);
      cr2yields_mllof+=mllOscon->GetBinContent(i);
      cr2yields_mllofem+=mllOsconEM->GetBinContent(i);
      cr2yields_mllofme+=mllOsconME->GetBinContent(i);
    }
    if(mllsigEE->GetBinCenter(i)>=80 && mllsigEE->GetBinCenter(i)<=100) {
      cout << "Bin center at " << mllsigEE->GetBinCenter(i) << " dy " << zlineshapeControl->GetBinContent(i) << " , of " << mllOscon->GetBinContent(i) << " , sf " << mllsconEE->GetBinContent(i)+mllsconMM->GetBinContent(i) << endl;
      pyields_mllee+=mllsconEE->GetBinContent(i);
      pyields_mllmm+=mllsconMM->GetBinContent(i);
      pyields_mllsf+=mllsconEE->GetBinContent(i);
      pyields_mllsf+=mllsconMM->GetBinContent(i);
      pyields_mllof+=mllOscon->GetBinContent(i);
      pyields_mllofem+=mllOsconEM->GetBinContent(i);
      pyields_mllofme+=mllOsconME->GetBinContent(i);
      pyields_mlldy+=zlineshapeControl->GetBinContent(i);
      pyields_mllpred+=mllOscon->GetBinContent(i)+zlineshapeControl->GetBinContent(i);
    }
  }
  
  dout << "     CR1::Observed      : " << cr1yields_mllsf << " (" << cr1yields_mllee << "ee , " << cr1yields_mllmm << "mm )" << endl;
  dout << "     CR1::Predicted(OF) : " << cr1yields_mllof << " (from OF, " << cr1yields_mllofem << " em, " << cr1yields_mllofme << " me) " << endl;
  dout << "     CR1::Ratio         : " << cr1yields_mllsf/cr1yields_mllof << " +/- " << (cr1yields_mllsf/cr1yields_mllof) * sqrt(1.0/cr1yields_mllsf + 1.0/cr1yields_mllof) << endl;
  
  dout << "     CR2::Observed      : " << cr2yields_mllsf << " (" << cr2yields_mllee << "ee , " << cr2yields_mllmm << "mm )" << endl;
  dout << "     CR2::Predicted(OF) : " << cr2yields_mllof << " (from OF, " << cr2yields_mllofem << " em, " << cr2yields_mllofme << " me) " << endl;
  dout << "     CR2::Ratio         : " << cr2yields_mllsf/cr2yields_mllof << " +/- " << (cr2yields_mllsf/cr2yields_mllof) * sqrt(1.0/cr2yields_mllsf + 1.0/cr2yields_mllof) << endl;

  dout << "     Peak::Observed     : " << pyields_mllsf << " (" << pyields_mllee << "ee , " << pyields_mllmm << "mm )" << endl;
  dout << "     Peak::Predicted    : " << pyields_mllpred << " (from OF, " << pyields_mllofem << " em, " << pyields_mllofme << " me --> " << pyields_mllof << " em/me and DY " << pyields_mlldy << " ) " << endl;

  NiceSummary << "     CR1::Observed      : " << cr1yields_mllsf << " (" << cr1yields_mllee << "ee , " << cr1yields_mllmm << "mm )" << endl;
  NiceSummary << "     CR1::Predicted(OF) : " << cr1yields_mllof << " (from OF) " << endl;
  NiceSummary << "     CR1::Ratio         : " << cr1yields_mllsf/cr1yields_mllof << " +/- " << (cr1yields_mllsf/cr1yields_mllof) * sqrt(1.0/cr1yields_mllsf + 1.0/cr1yields_mllof) << endl;
  
  NiceSummary << "     CR2::Observed      : " << cr2yields_mllsf << " (" << cr2yields_mllee << "ee , " << cr2yields_mllmm << "mm )" << endl;
  NiceSummary << "     CR2::Predicted(OF) : " << cr2yields_mllof << " (from OF) " << endl;
  NiceSummary << "     CR2::Ratio         : " << cr2yields_mllsf/cr2yields_mllof << " +/- " << (cr2yields_mllsf/cr2yields_mllof) * sqrt(1.0/cr2yields_mllsf + 1.0/cr2yields_mllof) << endl;
  
  //**************************

    
  TH1F *subtracted = (TH1F*)mllsig->Clone("subtracted");
  TH1F *subtractedControl = (TH1F*)mllscon->Clone("subtractedControl");
  TH1F *baseline = (TH1F*)mllOsig->Clone("baseline");
  TH1F *baselineControl = (TH1F*)mllOscon->Clone("baselineControl");
  for(int i=1;i<=subtracted->GetNbinsX();i++) {
    subtracted->SetBinContent(i,mllsig->GetBinContent(i)-mllOsig->GetBinContent(i));
    subtractedControl->SetBinContent(i,mllscon->GetBinContent(i)-mllOscon->GetBinContent(i));
    baseline->SetBinContent(i,0);
    baselineControl->SetBinContent(i,0);
  }
  
  TH1F *prediction = (TH1F*)mllOsig->Clone("prediction");
  prediction->Add(zlineshape);
  prediction->SetLineColor(TColor::GetColor("#CF35CA"));
  
  TH1F *control_prediction = (TH1F*)mllOscon->Clone("control_prediction");
  control_prediction->SetLineColor(TColor::GetColor("#CF35CA"));
  
  prediction->SetLineColor(TColor::GetColor("#cc0066"));
  mllOsig->SetLineColor(TColor::GetColor("#0000cc"));
  
  if(!doMC) mllOsig->SetLineStyle(2);
  zlineshape->SetLineColor(TColor::GetColor("#006600"));
  zlineshape->SetFillColor(TColor::GetColor("#006600"));
  zlineshape->SetFillStyle(3002); // light dots, not crushing other information
  if(!doMC) {
    mllOsig->SetLineWidth(2);
    prediction->SetLineWidth(2);
    zlineshape->SetLineWidth(2);
  }

  // FIX Y RANGE TO EASE COMPARISON
  mllsig->SetMaximum(ymax);
  float PreviousMinimum=mllsig->GetMinimum();
  mllsig->SetMinimum(0);
  mllsigEE->SetMaximum(ymax);
  mllsigMM->SetMaximum(ymax);
  mllOsig->SetMaximum(ymax);
  mllOscon->SetMaximum(ymax);
  subtracted->SetMaximum(60);
  subtracted->SetMinimum(-30);
  subtractedControl->SetMaximum(65);
  subtractedControl->SetMinimum(-30);


  // 1.- Signal region comparison
  TBox *srbox = new TBox(20,0,70,mllsig->GetMaximum());
  srbox->SetFillStyle(0);
  srbox->SetLineColor(TColor::GetColor("#298A08"));
  srbox->SetLineWidth(3);
    

  stringstream MetHeader;
  MetHeader << "N_{jets}#geq" << snjets.str() << ", E_{T}^{miss}>" << cutat << " GeV";
  stringstream MetHeaderCon;
//   MetHeaderCon << "N_{j}=2, N_{b}>0, 75<MET<150 GeV";
  MetHeaderCon << "N_{jets}=2, 100<E_{T}^{miss}<150 GeV";
  stringstream saveasSig;
  saveasSig << "MetPlots/mll_sig" << cutat << "__" << name;
  
  TLegend* leg;
  
  
  srbox->SetLineColor(TColor::GetColor("#00cc33"));

  
  if ( !doMC ) {
    TPad* rcan = new TPad("rcan","rcan",0,0,1,1);
    rcan->cd();
    mllsig->Draw();
    TGraphErrors *stat3j = MakeErrorGraph(prediction);
    stat3j->Draw("2,same");
    mllOsig->Draw("histo,same");
    zlineshape->Draw("histo,same");
    prediction->Draw("histo,same");
    mllsig->Draw("same");
    DrawPrelim();
    leg = make_legend();
    leg->SetX1(0.52);
    leg->SetHeader(MetHeader.str().c_str());
    leg->AddEntry(mllsig,"Data","PL");
    leg->AddEntry(prediction,"All bg prediction","L");
    leg->AddEntry(mllOsig,"bg without Z","L");
    if(!UseSpecialZprediction) leg->AddEntry(zlineshape,"Z lineshape","L");
    else leg->AddEntry(zlineshape,"bg with Z (JZB)","L");
    leg->AddEntry(stat3j,"stat. uncert.","F");
    leg->AddEntry(srbox,"SR","F");
    leg->Draw();
    srbox->Draw();
    sel->Draw();
    Save_With_Ratio_And_Line( mllsig, prediction, rcan, saveasSig.str() );
    
    //now also add systematic as a nice touch :-) 
    TPad* rcan2 = new TPad("rcan2","rcan2",0,0,1,1);
    
    
    float ConsideredZWidth=20;
    if(Contains((const char*)Restrmasscut,")<10")) ConsideredZWidth=10;
    
    
    TLine *SRline   = new TLine(70,0,70,ymax);
    TLine *ZLowLine = new TLine(91.2-ConsideredZWidth,0,91.2-ConsideredZWidth,ymax);
    TLine *ZHiLine  = new TLine(91.2+ConsideredZWidth,0,91.2+ConsideredZWidth,ymax);
    SRline->SetLineStyle(2);
    ZLowLine->SetLineStyle(2);
    ZHiLine->SetLineStyle(2);
    SRline->SetLineColor(kGray+2);
    ZLowLine->SetLineColor(kGray+2);
    ZHiLine->SetLineColor(kGray+2);
      
    rcan2->cd();
    mllsig->Draw();
    TH1F *SystHisto = (TH1F*)prediction->Clone("SystHisto");
    TGraphErrors *stat3jS = MakeErrorGraphSystematicAndStatistical(mllOsig,zlineshape,prediction,SystHisto);
    stat3jS->Draw("2,same");
    zlineshape->Draw("histo,same");
    prediction->Draw("histo,same");
    mllsig->Draw("same");
    DrawPrelim();
    leg = make_legend();
    leg->SetX1(0.52);
    leg->SetHeader(MetHeader.str().c_str());
    leg->AddEntry(mllsig,"Data","PL");
    leg->AddEntry(prediction,"Total backgrounds","L");
    if(!UseSpecialZprediction) leg->AddEntry(zlineshape,"DY (scaled) ","FL");
    else leg->AddEntry(zlineshape,"DY (JZB)","F");
    leg->AddEntry(stat3jS,"Total uncert.","F");
    leg->Draw();
    sel->Draw();
    SRline->Draw();
    ZLowLine->Draw();
    ZHiLine->Draw();
    

    save_with_ratio_and_sys_band(ConsideredZWidth,mllsig, prediction, rcan2, (saveasSig.str()+"__WithSys"), false, false, "data/pred",SystHisto );
    
    mllsig->GetYaxis()->SetRangeUser(0.1,ymax);
    
    TPad *rcan3 = new TPad("rcan3","rcan3",0,0,1,1);
    rcan3->cd();
    rcan3->SetLogy(1);
    rcan3->cd(); //need to switch back to pad (otherwise it's blank for some reason)
    mllsig->Draw();
    stat3jS->Draw("2,same");
    zlineshape->Draw("histo,same");
    prediction->Draw("histo,same");
    mllsig->Draw("same");
    DrawPrelim();
    sel->Draw();
    save_with_ratio_and_sys_band(ConsideredZWidth, mllsig,   prediction,  rcan3,   (saveasSig.str()+"__WithSys___LOG"), false, false, "data/pred",SystHisto );
    
  } else {
    TPad* rcan = new TPad("rcan","rcan",0,0,1,1);
    rcan->cd();
    mllsig->Draw();
    mcMllsig->Draw("same,hist");
    prediction->Draw("histo,same");
    mllsig->Draw("same");
    DrawPrelim();
    leg = allsamples.allbglegend();
    leg->SetHeader(MetHeader.str().c_str());
    leg->SetX1(0.52);
    leg->AddEntry(prediction,"All bg prediction","L");
    leg->AddEntry(srbox,"SR","F");
    leg->Draw();
    srbox->Draw();
    sel->Draw();
    Save_With_Ratio_And_Line( mllsig, *mcMllsig, rcan, saveasSig.str() );
    
    ExtractScaleFactor(mllsig,mllOsig,mcMllsig,mcMllOsig,prediction,leg,saveasSig.str(),srbox);
  }
  
  // 1b. MC: split ee and mumu
  if ( doMC ) {
    TPad* rcan = new TPad("rcan","rcan",0,0,1,1);
    rcan->cd();
    mllsigEE->Draw();
    mcMllsigEE->Draw("same,hist");
    mllsigEE->Draw("same");
    DrawPrelim();
    leg->Draw();
    srbox->Draw();
    sel->Draw();
    Save_With_Ratio_And_Line( mllsigEE, *mcMllsigEE,rcan->cd(),saveasSig.str()+"_ee" );

    rcan = new TPad("rcan","rcan",0,0,1,1);
    rcan->cd();
    mllsigMM->Draw();
    mcMllsigMM->Draw("histo,same");
    mllsigMM->Draw("same");
    DrawPrelim();
    leg->Draw();
    srbox->Draw();
    sel->Draw();
    Save_With_Ratio_And_Line( mllsigMM, *mcMllsigMM,rcan,saveasSig.str()+"_mm" );
  }

  // 1c. MC: compare of and sf
  if ( doMC ) {
    TH1F* hMcMllsig = CollapseStack( *mcMllsig);
    leg = allsamples.allbglegend("");
    leg->SetHeader(MetHeader.str().c_str());
    // Change "Data" label by hand
    ((TLegendEntry*)leg->GetListOfPrimitives()->At(0))->SetLabel("Same-flavor (MC)");
    TPad* rcan = new TPad("rcan","rcan",0,0,1,1);
    rcan->cd();
    hMcMllsig->SetMaximum(ymax);
    hMcMllsig->Draw("E");
    mcMllOsig->Draw("same,hist");
    hMcMllsig->Draw("same,E");
    DrawMCPrelim();
    leg->SetX1(0.52);
    leg->AddEntry(srbox,"SR","F");
    leg->Draw();
    srbox->Draw();
    sel->Draw();
    Save_With_Ratio_And_Line( hMcMllsig, *mcMllOsig, rcan, saveasSig.str()+"_mconly");
    
  }
  
  // 2.- Signal region comparison - LOG scale
  if ( !doMC ) {
    tcan->cd();
    mllsig->SetMinimum(0.2); // FIX Y RANGE TO EASE COMPARISON
    //mllsig->SetMaximum(mllsig->GetMaximum()*4.0);
    srbox->SetY2(mllsig->GetMaximum());
    tcan->SetLogy(1);
    stringstream saveasSig2;
    saveasSig2 << "MetPlots/mll_sig_ZLINESHAPE_" << cutat << "__" << name;
    
    sel->Draw();
    CompleteSave(tcan,saveasSig2.str());
    tcan->SetLogy(0);
  }
  

  // 3.- Signal region, background subtracted
  if ( !doMC ) {
    tcan->cd();
    for(int i=1;i<=subtracted->GetNbinsX();i++) {
      subtracted->SetBinContent(i,subtracted->GetBinContent(i)-zlineshape->GetBinContent(i));
      subtractedControl->SetBinContent(i,subtractedControl->GetBinContent(i)-zlineshapeControl->GetBinContent(i));
    }

    TGraphErrors *subtrerr = MakeErrorGraph(baseline);
    subtracted->Draw();
    subtrerr->Draw("2,same");
    subtracted->Draw("same");
    DrawPrelim();
    TLegend *DiffLeg = make_legend();
    DiffLeg->SetX1(0.4);
    DiffLeg->SetFillStyle(0);
    DiffLeg->SetHeader(MetHeader.str().c_str());
    DiffLeg->AddEntry(subtracted,"observed - predicted","PL");
    DiffLeg->AddEntry(subtrerr,"stat. uncert","F");
    DiffLeg->AddEntry((TObject*)0,"","");
    DiffLeg->AddEntry((TObject*)0,"","");
    DiffLeg->Draw();
    sel->Draw();
    
    stringstream saveasSigSub;
    saveasSigSub << "MetPlots/mll_sig_SUBTRACTED_" << cutat << "__" << name;
  
    //CompleteSave(tcan,saveasSigSub.str());
      
    // 3a.- Control region, background subtracted
    TGraphErrors *subtrerrControl = MakeErrorGraph(baselineControl);
    subtractedControl->Draw();
    subtrerrControl->Draw("2,same");
    subtractedControl->Draw("same");
    DrawPrelim();
    DiffLeg->SetHeader(MetHeaderCon.str().c_str());
    DiffLeg->Draw();
    saveasSigSub.str("");
    saveasSigSub << "MetPlots/mll_con_SUBTRACTED_" << cutat << "__" << name;
    //CompleteSave(tcan,saveasSigSub.str());
      
      
    // 4.- Signal region, background subtracted, errors added in quadrature
    TGraphErrors *subtrerr2 = (TGraphErrors*)subtrerr->Clone("subtrerr2");
    for(int i=1;i<=subtrerr2->GetN();i++) {
      subtrerr2->SetPoint(i-1,subtracted->GetBinCenter(i),subtracted->GetBinContent(i));
      subtrerr2->SetPointError(i-1,subtrerr2->GetErrorX(i),TMath::Sqrt(subtrerr2->GetErrorY(i)*subtrerr2->GetErrorY(i)+subtracted->GetBinError(i)*subtracted->GetBinError(i)));
    }
    TLine* l = new TLine(subtracted->GetBinLowEdge(1),0.,subtracted->GetBinLowEdge(subtracted->GetNbinsX()-1)+subtracted->GetBinWidth(1),0.);
    l->SetLineWidth(subtracted->GetLineWidth());
    subtracted->Draw();
    subtrerr2->Draw("2,same");
    l->Draw("same");
    subtracted->Draw("same");
    DrawPrelim();
    TLegend *DiffLeg2 = make_legend();
    DiffLeg2->SetX1(0.4);
    DiffLeg2->SetHeader(MetHeader.str().c_str());
    DiffLeg2->SetFillStyle(0);
    DiffLeg2->AddEntry(subtracted,"observed - predicted","PL");
    DiffLeg2->AddEntry(subtrerr2,"stat. uncert","F");
    DiffLeg2->AddEntry((TObject*)0,"","");
    DiffLeg2->AddEntry((TObject*)0,"","");
    DiffLeg2->Draw();
    
    stringstream saveasSigSub2;
    saveasSigSub2 << "MetPlots/mll_sig_SUBTRACTED_quadr_" << cutat << "__" << name;

    sel->Draw();
    CompleteSave(tcan,saveasSigSub2.str());


    //4a.- Control region, background subtracted, errors added in quadrature
    TGraphErrors *subtrerr2Control = (TGraphErrors*)subtrerrControl->Clone("subtrerr2Control");
    for(int i=1;i<=subtrerr2Control->GetN();i++) {
          subtrerr2Control->SetPoint(i-1,subtractedControl->GetBinCenter(i),subtractedControl->GetBinContent(i));
          float width=subtrerr2Control->GetErrorX(i);
          if(i==subtrerr2Control->GetN()) width=subtrerr2Control->GetErrorX(i-1);
          subtrerr2Control->SetPointError(i-1,width,TMath::Sqrt(subtrerr2Control->GetErrorY(i)*subtrerr2Control->GetErrorY(i)+subtractedControl->GetBinError(i)*subtractedControl->GetBinError(i)));
    }
    TLine* lControl = new TLine(subtractedControl->GetBinLowEdge(1),0.,subtractedControl->GetBinLowEdge(subtractedControl->GetNbinsX()-1)+subtractedControl->GetBinWidth(1),0.);
    lControl->SetLineWidth(subtractedControl->GetLineWidth());
    subtractedControl->Draw();
    subtrerr2Control->Draw("2,same");
    lControl->Draw("same");
    subtractedControl->Draw("same");
    DrawPrelim();
    DiffLeg2->SetHeader(MetHeaderCon.str().c_str());
    DiffLeg2->Draw();
      
    saveasSigSub2.str("");
    saveasSigSub2 << "MetPlots/mll_con_SUBTRACTED_quadr_" << cutat << "__" << name;
      
    sel->Draw();
    CompleteSave(tcan,saveasSigSub2.str());
      
    delete DiffLeg;
    delete DiffLeg2;
      
  } // !doMC
    
    
  // 5.- Control region comparison
//   scalefactor = (mllscon->Integral(scaleBinLow,scaleBinHigh)-mllOscon->Integral(scaleBinLow,scaleBinHigh));
//   scalefactor /= zlineshape->Integral(scaleBinLow,scaleBinHigh);
//   zlineshape->Scale(scalefactor);
  control_prediction->Add(zlineshapeControl);

  control_prediction->SetMaximum(ymax); // FIX MAXIMUM TO EASE COMPARISON
  control_prediction->SetMinimum(0);

  TBox *cr1box = new TBox(20,0,70,control_prediction->GetMaximum());
  cr1box->SetFillStyle(0);
  cr1box->SetLineColor(TColor::GetColor("#0404B4"));
  cr1box->SetLineWidth(3);
  
  TBox *cr2box = new TBox(120,0,xmax,control_prediction->GetMaximum());
  cr2box->SetFillStyle(0);
  cr2box->SetLineColor(TColor::GetColor("#0404B4"));
  cr2box->SetLineWidth(3);
  cr2box->SetLineStyle(2);
  
  stringstream saveasCon;
  saveasCon << "MetPlots/mll_con" << cutat << "__" << name;

  TLegend *legc;
  //control_prediction->GetYaxis()->SetRangeUser(0,control_prediction->GetMaximum()*1.3);
  if ( !doMC ) {
    TPad* rcan = new TPad("rcan","rcan",0,0,1,1);
    rcan->cd();
    Color_t control_prediction_color = control_prediction->GetLineColor();
    int LineWidth = control_prediction->GetLineWidth();
    control_prediction->SetLineColor(TColor::GetColor("#FF4000"));
    control_prediction->SetLineWidth(2);

    TH1F *ControlSystHisto = (TH1F*)control_prediction->Clone("SystHisto");
    TGraphErrors *stat3jS = MakeErrorGraphSystematicAndStatistical(mllOscon,zlineshapeControl,control_prediction,ControlSystHisto);
    control_prediction->Draw("hist");
    stat3jS->Draw("2,same");
    zlineshapeControl->Draw("histo,same");
    control_prediction->Draw("histo,same");
    mllscon->Draw("same");
    DrawPrelim();
    legc = make_legend();
    legc->SetX1(0.52);
    legc->SetHeader(MetHeaderCon.str().c_str());
    legc->AddEntry(mllscon,"Data","PL");
    legc->AddEntry(control_prediction,"Total backgrounds","L");
    legc->AddEntry(zlineshapeControl,"DY (scaled)","FL");
    legc->AddEntry(stat3jS,"Total uncert.","F");
    legc->AddEntry(cr1box,"CR1","F");
    legc->AddEntry(cr2box,"CR2","F");
    legc->Draw();
    cr1box->Draw();
    cr2box->Draw();
    sel->Draw();
    
    save_with_ratio_and_sys_band(ConsideredZWidth, mllscon, control_prediction, rcan, saveasCon.str()     , false, false, "data/pred",ControlSystHisto );
    
    control_prediction->SetLineColor(control_prediction_color);
    control_prediction->SetLineWidth(LineWidth);
  }  else {
    control_prediction->SetLineColor(TColor::GetColor("#FF4000"));
    TPad* rcan = new TPad("rcan","rcan",0,0,1,1);
    rcan->cd();
    control_prediction->Draw("hist");
    mcMllscon->Draw("same,hist");
    control_prediction->Draw("histo,same");
    mllscon->Draw("same");
    DrawPrelim();
    legc = allsamples.allbglegend();
    legc->SetX1(0.52);
    legc->SetHeader(MetHeaderCon.str().c_str());
    legc->AddEntry(control_prediction,"All bg","L");
    legc->AddEntry(cr1box,"CR1","F");
    legc->AddEntry(cr2box,"CR2","F");
    legc->Draw();
    cr1box->Draw();
    cr2box->Draw();
    sel->Draw();
    Save_With_Ratio_And_Line( mllscon, *mcMllscon, rcan, saveasCon.str());
  } 
  
  // 6. - Opposite-flavour data/MC comparison
  if ( doMC ) {
    TPad* rcan = new TPad("rcan","rcan",0,0,1,1);    
    rcan->cd();
    mllOsig->SetLineColor(kBlack);
    mllOsig->Draw();
    mcMllOsig->Draw("same,hist");
    mllOsig->Draw("same");
    TLegend *legsdm = allsamples.allbglegend();
    legsdm->SetHeader((MetHeader.str()+", OF").c_str());
    legsdm->SetX1(0.52);
    legsdm->Draw();
    stringstream saveasSigOF;
    saveasSigOF << "MetPlots/mll_sig_of_" << cutat << "__" << name;
    sel->Draw();
    Save_With_Ratio_And_Line( mllOsig, *mcMllOsig, rcan, saveasSigOF.str());

    rcan = new TPad("rcan","rcan",0,0,1,1); 
    rcan->cd();
    mllOscon->SetLineColor(kBlack);
    mllOscon->Draw();
    mcMllOscon->Draw("same,hist");
    mllOscon->Draw("same");
    TLegend *legcdm = allsamples.allbglegend();
    legcdm->SetHeader((MetHeaderCon.str()+", OF").c_str());
    legcdm->SetX1(0.52);
    legcdm->Draw();
    stringstream saveasConOF;
    saveasConOF << "MetPlots/mll_con_of_" << cutat << "__" << name;
    sel->Draw();
    Save_With_Ratio_And_Line( mllOscon, *mcMllOscon, rcan, saveasConOF.str());

    delete legsdm;
    delete legcdm;
  }
  
  // 7. - Opposite flavor data/MC comparison for pt (!)
  if ( doMC ) {
    TPad* rcan = new TPad("rcan","rcan",0,0,1,1);    
    rcan->cd();
    rcan->SetLogy(1);

    ptsig->SetLineColor(kBlack);
    ptsig->Draw();
    mcptsig->Draw("same,hist");
    ptsig->Draw("same");
    TLegend *legsdm = allsamples.allbglegend();
    legsdm->SetHeader((MetHeader.str()+", SF").c_str());
    legsdm->SetX1(0.52);
    legsdm->Draw();
    stringstream saveasSigOF2;
    saveasSigOF2 << "MetPlots/mll_sig_sf_PTdist_" << cutat << "__" << name;
    sel->Draw();
    Save_With_Ratio_And_Line( ptsig, *mcptsig, rcan, saveasSigOF2.str());

    delete legsdm;
  }

  
  // 8. - Opposite flavor data/MC comparison for pt (!)
  if ( doMC ) {
    TPad* rcan = new TPad("rcan","rcan",0,0,1,1);    
    rcan->cd();
    rcan->SetLogy(1);
    
    ptOsig->SetLineColor(kBlack);
    ptOsig->Draw();
    mcptOsig->Draw("same,hist");
    ptOsig->Draw("same");
    TLegend *legsdm = allsamples.allbglegend();
    legsdm->SetHeader((MetHeader.str()+", OF").c_str());
    legsdm->SetX1(0.52);
    legsdm->Draw();
    stringstream saveasSigOF3;
    saveasSigOF3 << "MetPlots/mll_sig_of_PTdist_" << cutat << "__" << name;
    sel->Draw();
    Save_With_Ratio_And_Line( ptOsig, *mcptOsig, rcan, saveasSigOF3.str());

    delete legsdm;
  }

  // 9. - Shape comparison between SR and CR
  if ( !doMC ) { // SF
    TH1F* scaled_conSF = (TH1F*)mllscon->Clone("scaled_conSF");
    scaled_conSF->SetLineColor(kBlue);
    scaled_conSF->Scale(mllsig->Integral()/scaled_conSF->Integral());
    TPad* rcan = new TPad("rcan","rcan",0,0,1,1);    
    rcan->cd();
    mllsig->Draw();
    scaled_conSF->Draw("same,hist");
    mllsig->Draw("same");
    TLegend *leg9 = make_legend("Same-flavor",0.5,0.7,false);
    leg9->SetHeader("Same-flavor");
    leg9->AddEntry(mllsig,"SR","pl");
    leg9->AddEntry(scaled_conSF,"CR (scaled)","l");
    leg9->Draw();
    DrawPrelim();
    stringstream saveas9;
    saveas9 << "MetPlots/mll_compSF_" << cutat << "__" << name;
    sel->Draw();
    Save_With_Ratio_And_Line( mllsig, scaled_conSF, rcan, saveas9.str());
    delete leg9;
  } else {
    TH1F* hMcMllsig    = CollapseStack( *mcMllsig,"hMcMllSig");
    TH1F* scaled_conSF = CollapseStack( *mcMllscon,"scaled_conSF");
    scaled_conSF->SetLineColor(kBlue);
    scaled_conSF->SetFillStyle(0);
    scaled_conSF->Scale(hMcMllsig->Integral()/scaled_conSF->Integral());
    TPad* rcan = new TPad("rcan","rcan",0,0,1,1);    
    rcan->cd();
    hMcMllsig->SetMaximum(ymax);
    hMcMllsig->Draw();
    scaled_conSF->Draw("same,hist");
    hMcMllsig->Draw("same");
    TLegend *leg9 = make_legend("Same-flavor MC",0.5,0.7,false);
    leg9->SetHeader("Same-flavor MC");
    leg9->AddEntry(hMcMllsig,"SF SR","pl");
    leg9->AddEntry(scaled_conSF,"SF CR (scaled)","l");
    leg9->Draw();
    DrawMCPrelim();
    stringstream saveas9;
    saveas9 << "MetPlots/mll_compSF_" << cutat << "__" << name;
    sel->Draw();
    Save_With_Ratio_And_Line( hMcMllsig, scaled_conSF, rcan, saveas9.str());
    delete leg9;
  }
  if ( !doMC ) { // OF
    TH1F* scaled_conOF = (TH1F*)control_prediction->Clone("scaled_conOF");
    scaled_conOF->SetLineColor(kBlue);
    scaled_conOF->Scale(mllOsig->Integral()/scaled_conOF->Integral());
    TPad* rcan = new TPad("rcan","rcan",0,0,1,1);
    rcan->cd();
    mllOsig->SetLineColor(kBlack);
    mllOsig->Draw();
    scaled_conOF->Draw("same,hist");
    mllOsig->Draw("same");
    TLegend *leg9 = make_legend("Opposite-flavor",0.5,0.7,false);
    leg9->AddEntry(mllOsig,"OF SR","pl");
    leg9->AddEntry(scaled_conOF,"OF CR (scaled)","l");
    leg9->Draw();
    DrawPrelim();
    stringstream saveas9;
    saveas9 << "MetPlots/mll_compOF_" << cutat << "__" << name;
    sel->Draw();
    Save_With_Ratio_And_Line( mllOsig, scaled_conOF, rcan, saveas9.str());

    delete leg9;
  } else { // SF MC
    TH1F* hMcMllOsig    = CollapseStack( *mcMllOsig, "hMcMllOsig");
    TH1F* scaled_conOF = CollapseStack( *mcMllOscon, "scaled_conOF");
    scaled_conOF->SetLineColor(kBlue);
    scaled_conOF->SetFillStyle(0);
    scaled_conOF->Scale(hMcMllOsig->Integral()/scaled_conOF->Integral());
    TPad* rcan = new TPad("rcan","rcan",0,0,1,1);    
    rcan->cd();
    hMcMllOsig->SetMaximum(ymax);
    hMcMllOsig->Draw();
    scaled_conOF->Draw("same,hist");
    hMcMllOsig->Draw("same");
    TLegend *leg9 = make_legend("Opposite-flavor MC",0.5,0.7,false);
    leg9->AddEntry(hMcMllOsig,  "OF SR","pl");
    leg9->AddEntry(scaled_conOF,"OF CR (scaled)","l");
    leg9->Draw();
    DrawMCPrelim();
    stringstream saveas9;
    saveas9 << "MetPlots/mll_compOF_" << cutat << "__" << name;
    sel->Draw();
    Save_With_Ratio_And_Line( hMcMllOsig, scaled_conOF, rcan, saveas9.str());
    delete leg9;
  }
  

  // Memory clean-up
  if (doMC) {
    delete mcMllscon;
    delete mcMllOscon;
    delete mcMllsig;
    delete mcMllsigEE;
    delete mcMllsigMM;
    delete mcMllOsig;
  }
  
  delete cr1box;
  delete cr2box;
  delete srbox;
  delete legc;
  delete leg;
  
  delete mllscon;
  delete mllsconEE;
  delete mllsconMM;
  delete mllOscon;
  delete mllOsconEM;
  delete mllOsconME;
  delete mllsig;
  delete mllsigEE;
  delete mllsigMM;
  delete mllOsig;
  delete mllOsigEM;
  delete mllOsigME;
  delete ptsig;
  delete ptOsig;
  delete zlineshape;
  delete Ozlineshape;
  delete zlineshapeControl;
  delete tcan;
  
  NiceSummary.close();
  
}  


void DoMetPlots(string datajzb, string mcjzb) {
  switch_overunderflow(true);
  float metCuts[] = { 100., 150. };
  //float ymax[] = { 180., 170. };
  float ymax[] = { 90., 140. };
  int jetCuts[] = { 3, 2 };
  string leptCuts[] = { "pt1>20&&pt2>20", "pt1>20&&pt2>20&&pfTightHT>100" };
  bool nomc(0),domc(1);
  string backup_basicqualitycut = (const char*) basicqualitycut;
  string backup_essentialcut    = (const char*) essentialcut;
  string backup_basiccut        = (const char*) basiccut;
  string backup_leptoncut       = (const char*) leptoncut;
  
  //zlineshapes(); // Rinout plots
  for ( int i=0; i<2; ++i ) {
    //need to make sure that the above changes actually have some effect. we therefore check all relevant cuts and
    //set the pt condition to 10/10 (yes you read that right). the addition cut (above) will therefore elevate it 
    // to 20,10 or 20,20. otherwise basicqualitycut will impose 20,20 ...
    
    bool isAachen=i;//1=Aachen, 0=not.
    string Sbasicqualitycut = backup_basicqualitycut;
    if(i==1) Sbasicqualitycut = ReplaceAll(Sbasicqualitycut,")<1.4",")<2.4");
    basicqualitycut=TCut(Sbasicqualitycut.c_str());
    
    string Sbasiccut = backup_basiccut;
    if(i==1) Sbasiccut = ReplaceAll(Sbasiccut,")<1.4",")<2.4");
    basiccut=TCut(Sbasiccut.c_str());
    
    string Sessentialcut = backup_essentialcut;
    if(i==1) Sessentialcut = ReplaceAll(Sessentialcut,")<1.4",")<2.4");
    essentialcut=TCut(Sessentialcut.c_str());
    
    string Sleptoncut = backup_leptoncut;
    if(i==1) Sleptoncut = ReplaceAll(Sleptoncut,")<1.4",")<2.4");
    if(i==1) leptoncut=TCut(Sleptoncut.c_str());
    
    cout << "Basic cut : " << (const char*) basiccut << endl;
    cout << "Essential cut : " << (const char*) essentialcut << endl;
  
    ProvideEEOverMMEstimate(cutOSSF&&TCut("pfJetGoodNum40==2")&&TCut(("mll>20&&"+leptCuts[i]).c_str()));
//    ProduceMetPlotsWithCut(isAachen, TCut(("mll>20&&"+leptCuts[i]).c_str()),"",metCuts[i],jetCuts[i],nomc,ymax[i]);
    ProduceMetPlotsWithCut(isAachen, TCut(("mll>20&&"+leptCuts[i]).c_str()),"",metCuts[i],jetCuts[i],domc,ymax[i]);
    
    ProduceMetPlotsWithCut(isAachen, TCut(("mll>20&&pfJetGoodNumBtag30==0&&"+leptCuts[i]).c_str()),"bTagVeto30",metCuts[i], jetCuts[i],nomc,ymax[i]);
    ProduceMetPlotsWithCut(isAachen, TCut(("mll>20&&pfJetGoodNumBtag30>0&&"+leptCuts[i]).c_str()),"AtLeastOneBJet30",metCuts[i],jetCuts[i],nomc,ymax[i]);
    ProduceMetPlotsWithCut(isAachen, TCut(("mll>20&&pfJetGoodNumBtag30==0&&"+leptCuts[i]).c_str()),"bTagVeto30",metCuts[i], jetCuts[i],domc,ymax[i]);
    
    ProduceMetPlotsWithCut(isAachen, TCut(("mll>20&&pfJetGoodNumBtag30>0&&"+leptCuts[i]).c_str()),"AtLeastOneBJet30",metCuts[i], jetCuts[i],domc,ymax[i]);
    
    
    ProduceMetPlotsWithCut(isAachen, TCut(("mll>20&&pfJetGoodNumBtag30>1&&"+leptCuts[i]).c_str()),"AtLeastTwoBJet30",metCuts[i],jetCuts[i],nomc,ymax[i]);
    ProduceMetPlotsWithCut(isAachen, TCut(("mll>20&&pfJetGoodNumBtag30>2&&"+leptCuts[i]).c_str()),"AtLeastThreeBJet30",metCuts[i],jetCuts[i],nomc,ymax[i]);

  }
  basicqualitycut=TCut(backup_basicqualitycut.c_str());
  basiccut       =TCut(backup_basiccut.c_str());
  essentialcut   =TCut(backup_essentialcut.c_str());
  leptoncut      =TCut(backup_leptoncut.c_str());
  switch_overunderflow(false);
}

void LabelHisto(TH1 *MET_ratio,string titlex, string titley) {
  MET_ratio->GetXaxis()->SetTitle(titlex.c_str());
  MET_ratio->GetXaxis()->CenterTitle();
  MET_ratio->GetYaxis()->SetTitle(titley.c_str());
  MET_ratio->GetYaxis()->CenterTitle();
}

TH1F* GetPredictedAndObservedMetShapes(TCut JetCut, string sPositiveCut,string sNegativeCut,string CorrectedMet,string ObservedMet, string JZBPosvar, string JZBNegvar, float MetCut, int is_data, bool isDYonly, bool isAachen) {
  
  //Steps: 
  // 1) Prepare samples and histo definition (with "optimal" binning for MET cut)
  // 2) Fill MET histograms
  // 3) Fill JZB histograms
  // 4) Draw them and store them
  // 5) return predicted MET distribution as is (i.e. not scaled by factor of 2!)
  
  dout << "*************************************" << endl;
//   cout << "** SUMMARY BEFORE STARTING DRAWING **" << endl;
//   cout << "MET variable: " << ObservedMet << endl;
//   cout << "Corr. MET var:" << CorrectedMet << endl;
//   cout << "JZB pos. var: " << JZBPosvar << endl;
//   cout << "JZB neg. var: " << JZBNegvar << endl;
//   cout << "JZB pos cut : " << sPositiveCut << endl;
//   cout << "JZB neg cut : " << sNegativeCut << endl;
  
  if(isAachen) MetPlotsSpace::Zprediction_Uncertainty=0.3;
  
  //Step 1: Prepare samples and histo definition 
  vector<int> SelectedSamples;
  if(is_data==mc&&isDYonly) {
    SelectedSamples=allsamples.FindSample("_em_");
    if(SelectedSamples.size()==0) {
      write_error(__FUNCTION__,"Cannot continue, there seems to be no DY sample without Taus - goodbye!");
      assert(SelectedSamples.size()>0);
    }
  }
  
  float DisplayedBinSize=10.0; // this is the bin size that we use for plotting
  
  float BinWidth=1.0;
  float xmin=0;
  float xmax=150;
  if(isAachen) xmax=160;
  if(MetCut>=xmax) xmax=MetCut+10;
  int nbins=int((xmax-xmin)/BinWidth);
  
  float pt2cut=20;
  if(isAachen)pt2cut=10;
  
  stringstream basiccut;
  basiccut << (const char*) JetCut << "&&" << (const char*) Restrmasscut << "&&" << (const char*) leptoncut << "&&pt1>20&&pt2>" << pt2cut;
  
  stringstream cMET_observed;
  cMET_observed << "(" << basiccut.str() << "&&(" << sPositiveCut << ")&&" << (const char*) cutOSSF << ")";
  stringstream cMET_ttbar_pred;
  cMET_ttbar_pred << "(" << basiccut.str() << "&&(" << sPositiveCut << ")&&" << (const char*) cutOSOF << ")";
  stringstream cMET_osof_pred;
  cMET_osof_pred << "(" << basiccut.str() << "&&(" << sNegativeCut << ")&&" << (const char*) cutOSOF << ")";
  stringstream cMET_ossf_pred;
  cMET_ossf_pred << "(" << basiccut.str() << "&&(" << sNegativeCut << ")&&" << (const char*) cutOSSF << ")";
  
  //Step 2: Fill Met histograms
  float bottommargin=gStyle->GetPadBottomMargin();
  float canvas_height=gStyle->GetCanvasDefH();
  float canvas_width=gStyle->GetCanvasDefW();
  float ratiospace=0.25;// space the ratio should take up (relative to original pad)
  
  float ratiobottommargin=0.3;
  float ratiotopmargin=0.1;
  
  float xstretchfactor=((1-ratiospace)*(1-gStyle->GetPadTopMargin()))/((1)*ratiospace);
  
  TCanvas *main_canvas = new TCanvas("main_canvas","main_canvas",(Int_t)canvas_width,(Int_t)(canvas_height*(1+ratiospace)));
  TPad *mainpad = new TPad("mainpad","mainpad",0,1-(1.0/(1+ratiospace)),1,1);//top (main) pad
  TPad *coverpad = new TPad("coverpad","coverpad",gStyle->GetPadLeftMargin()-0.008,1-(1.0/(1+ratiospace))-0.04,1,1-(1.0/(1+ratiospace))+0.103);//pad covering up the x scale
  TPad *bottompad = new TPad("bottompad", "Ratio Pad",0,0,1,(1-(1-bottommargin)/(1+ratiospace))-0.015); //bottom pad
  
  main_canvas->Range(0,0,1,1);
  main_canvas->SetBorderSize(0);
  main_canvas->SetFrameFillColor(0);
  
  mainpad->Draw();
  mainpad->cd();
  mainpad->SetLogy(1);
  mainpad->Range(0,0,1,1);
  mainpad->SetFillColor(kWhite);
  mainpad->SetBorderSize(0);
  mainpad->SetFrameFillColor(0);

  
  TH1F *MET_observed  = allsamples.Draw("MET_observed",ObservedMet,nbins,xmin,xmax,"MET [GeV]","events",
                                       TCut(cMET_observed.str().c_str()),is_data,PlottingSetup::luminosity,SelectedSamples);
  TH1F *MET_ossf_pred = allsamples.Draw("MET_ossf_pred",CorrectedMet,nbins,xmin,xmax,"MET [GeV]","events",
                                        TCut(cMET_ossf_pred.str().c_str()),is_data,PlottingSetup::luminosity,SelectedSamples);
  TH1F *MET_osof_pred = allsamples.Draw("MET_osof_pred",CorrectedMet,nbins,xmin,xmax,"MET [GeV]","events",
                                        TCut(cMET_osof_pred.str().c_str()),is_data,PlottingSetup::luminosity,SelectedSamples);
  TH1F *MET_ttbar_pred= allsamples.Draw("MET_ttbar_pred",ObservedMet,nbins,xmin,xmax,"MET [GeV]","events",
                                        TCut(cMET_ttbar_pred.str().c_str()),is_data,PlottingSetup::luminosity,SelectedSamples);
  
  
  if((isDYonly && is_data==mc) || is_data==data) {
    TH1F *MET_truth  = allsamples.Draw("MET_truth",ObservedMet,1,MetCut,10000,"MET [GeV]","events",TCut(((string)"met[4]/MetFactor>"+any2string(MetCut)).c_str())&&cutOSSF&&TCut(basiccut.str().c_str()),is_data,PlottingSetup::luminosity,SelectedSamples);
    TH1F *eeMET_truth  = allsamples.Draw("eeMET_truth",ObservedMet,1,MetCut,10000,"MET [GeV]","events",TCut(((string)"id1==0 && met[4]/MetFactor>"+any2string(MetCut)).c_str())&&cutOSSF&&TCut(basiccut.str().c_str()),is_data,PlottingSetup::luminosity,SelectedSamples);
    TH1F *mmMET_truth  = allsamples.Draw("mmMET_truth",ObservedMet,1,MetCut,10000,"MET [GeV]","events",TCut(((string)"id1==1 && met[4]/MetFactor>"+any2string(MetCut)).c_str())&&cutOSSF&&TCut(basiccut.str().c_str()),is_data,PlottingSetup::luminosity,SelectedSamples);
    TH1F *MET_otruth  = allsamples.Draw("MET_otruth",ObservedMet,1,MetCut,10000,"MET [GeV]","events",TCut(((string)"met[4]/MetFactor>"+any2string(MetCut)).c_str())&&cutOSOF&&TCut(basiccut.str().c_str()),is_data,PlottingSetup::luminosity,SelectedSamples);
    if(is_data==mc) write_info(__FUNCTION__,"In Z peak: DY Truth is : "+any2string(MET_truth->Integral()));
    if(is_data==data) {
        ofstream NiceSummary;
        NiceSummary.open("Summary.txt",ios::app);
        write_info(__FUNCTION__,"In Z peak: Observed : " +any2string(MET_truth->Integral()) + "( ee: "+any2string(eeMET_truth->Integral()) + " , mm: "+any2string(mmMET_truth->Integral())+" )");
        write_info(__FUNCTION__,"In Z peak: TTbar est: " +any2string(MET_otruth->Integral()));
        NiceSummary << "  **** Z peak estimation **** " << endl;
        NiceSummary << "Z peak observed : " << any2string(MET_truth->Integral()) <<  "( ee: " <<any2string(eeMET_truth->Integral()) << " , mm: "<<any2string(mmMET_truth->Integral()) << " )" << endl;
        NiceSummary << "Z peak ttbar pred: " << any2string(MET_otruth->Integral()) << endl;
        NiceSummary.close();
    }
    delete MET_truth;
    delete MET_otruth;
    delete eeMET_truth;
    delete mmMET_truth;
  }
    
//  write_info(__FUNCTION__,"Full cut!");
//  cout << (const char*)(TCut(((string)"met[4]/MetFactor>"+any2string(MetCut)).c_str())&&cutOSSF&&TCut(basiccut.str().c_str())) << endl;

  TH1F *MET_predicted=(TH1F*)MET_ossf_pred->Clone("MET_predicted");
  MET_predicted->Add(MET_osof_pred,-1);
  MET_predicted->Add(MET_ttbar_pred);
  MET_predicted->SetLineColor(kRed);
  MET_observed->SetLineColor(kBlack);
  
  TH1F *MET_Z_prediction=(TH1F*)MET_ossf_pred->Clone("MET_Z_prediction");
  MET_Z_prediction->Add(MET_osof_pred,-1);
  MET_Z_prediction->SetLineColor(kBlue);
  
  LabelHisto(MET_observed,"MET (GeV)","events");
  
  //Step 3: Fill JZB histograms
  
  TH1F *JZB_observed  = allsamples.Draw("JZB_observed",JZBPosvar,nbins,xmin,xmax,"JZB [GeV]","events",
                                       TCut(cMET_observed.str().c_str()),is_data,PlottingSetup::luminosity,SelectedSamples);
  TH1F *JZB_ossf_pred = allsamples.Draw("JZB_ossf_pred",JZBNegvar,nbins,xmin,xmax,"JZB [GeV]","events",
                                        TCut(cMET_ossf_pred.str().c_str()),is_data,PlottingSetup::luminosity,SelectedSamples);
  TH1F *JZB_osof_pred = allsamples.Draw("JZB_osof_pred",JZBNegvar,nbins,xmin,xmax,"JZB [GeV]","events",
                                        TCut(cMET_osof_pred.str().c_str()),is_data,PlottingSetup::luminosity,SelectedSamples);
  TH1F *JZB_ttbar_pred= allsamples.Draw("JZB_ttbar_pred",JZBPosvar,nbins,xmin,xmax,"JZB [GeV]","events",
                                        TCut(cMET_ttbar_pred.str().c_str()),is_data,PlottingSetup::luminosity,SelectedSamples);
  
  TH1F *JZB_predicted=(TH1F*)JZB_ossf_pred->Clone("JZB_predicted");
  JZB_predicted->Add(JZB_osof_pred,-1);
  JZB_predicted->Add(JZB_ttbar_pred);
  JZB_predicted->SetLineColor(kRed);
  JZB_observed->SetLineColor(kBlack);
  
  TH1F *JZB_Z_prediction=(TH1F*)JZB_ossf_pred->Clone("JZB_Z_prediction");
  JZB_Z_prediction->Add(JZB_osof_pred,-1);
  MET_Z_prediction->SetLineColor(kBlue);
  
  LabelHisto(JZB_observed,"JZB (GeV)","events");
  
  // Step 4: Draw them and store them
  
  TLegend *legend = new TLegend(0.6,0.6,0.89,0.89);
  
  MET_ttbar_pred->SetLineColor(TColor::GetColor("#005C00"));
  JZB_ttbar_pred->SetLineColor(TColor::GetColor("#005C00"));
  
  legend->SetFillColor(kWhite);
  legend->SetBorderSize(0);
  legend->AddEntry(MET_predicted,"prediction","l");
  legend->AddEntry(MET_observed,"observed","p");
  legend->AddEntry(MET_Z_prediction,"predicted Z","l");
  legend->AddEntry(MET_ttbar_pred,"OF-based prediction","l");
  
  if(is_data==mc) legend->SetHeader("Simulation:");
  if(is_data==mc&&isDYonly) legend->SetHeader("DY #rightarrow ee,#mu#mu only:");
  if(is_data==data) legend->SetHeader("Data:");
  
  stringstream SaveJZBname;
  stringstream SaveMETname;
  if(is_data==data) {
    SaveJZBname << "MetPrediction/JZBdistribution_Data_METCutAt" << MetCut;
    SaveMETname << "MetPrediction/METdistribution_Data_METCutAt" << MetCut;
  }
  if(is_data==mc&&!isDYonly) {
    SaveJZBname << "MetPrediction/JZBdistribution_FullMC_METCutAt" << MetCut;
    SaveMETname << "MetPrediction/METdistribution_FullMC_METCutAt" << MetCut;
  }
  if(is_data==mc&&isDYonly) {
    SaveJZBname << "MetPrediction/JZBdistribution_DYMC_METCutAt" << MetCut;
    SaveMETname << "MetPrediction/METdistribution_DYMC_METCutAt" << MetCut;
  }
  
  dout << "Shape summary (MET>50)  for ";
  if(is_data==data) dout << "data";
  if(is_data==mc&&isDYonly) dout<< "DY ";
  if(is_data==mc&&!isDYonly) dout << " Full MC";
  dout << " : " << endl;
  
  dout << "   observed  : " << MET_observed->Integral(MET_observed->FindBin(50),MET_observed->FindBin(xmax))     << endl;
  dout << "   predicted : " << MET_predicted->Integral(MET_predicted->FindBin(50),MET_predicted->FindBin(xmax))    << endl;
  dout << "     Z pred  : " << MET_Z_prediction->Integral(MET_Z_prediction->FindBin(50),MET_Z_prediction->FindBin(xmax)) << endl;
  dout << "     ttbar   : " << MET_ttbar_pred->Integral(MET_ttbar_pred->FindBin(50),MET_ttbar_pred->FindBin(xmax))   << endl;

  TH1F *ZpredClone = (TH1F*)MET_Z_prediction->Clone("ZpredClone");
  ZpredClone->SetLineColor(kBlue);
  MET_observed->Rebin(int(DisplayedBinSize/BinWidth));
  ZpredClone->Rebin(int(DisplayedBinSize/BinWidth));
  MET_predicted->Rebin(int(DisplayedBinSize/BinWidth));
  MET_ttbar_pred->Rebin(int(DisplayedBinSize/BinWidth));
  
  TH1F *JZBZpredClone = (TH1F*)JZB_Z_prediction->Clone("ZpredClone");
  JZBZpredClone->SetLineColor(kBlue);
  JZB_observed->Rebin(int(DisplayedBinSize/BinWidth));
  JZBZpredClone->Rebin(int(DisplayedBinSize/BinWidth));
  JZB_predicted->Rebin(int(DisplayedBinSize/BinWidth));
  JZB_ttbar_pred->Rebin(int(DisplayedBinSize/BinWidth));
  
  TH1F *JZB_ratio = (TH1F*)JZB_observed->Clone("JZB_ratio");
  JZB_ratio->Divide(JZB_predicted);
  LabelHisto(JZB_ratio,"JZB (GeV)","obs/pred");
  TH1F *MET_ratio = (TH1F*)MET_observed->Clone("MET_ratio");
  MET_ratio->Divide(MET_predicted);
  MET_observed->SetMaximum(5*MET_observed->GetMaximum());
  JZB_observed->SetMaximum(5*JZB_observed->GetMaximum());
  MET_observed->SetMinimum(0.5);
  JZB_observed->SetMinimum(0.5);
  LabelHisto(MET_ratio,"MET (GeV)","obs/pred");
  TBox *sysenvelope = new TBox(xmin,1.0-MetPlotsSpace::Zprediction_Uncertainty,xmax,1.0+MetPlotsSpace::Zprediction_Uncertainty);
  sysenvelope->SetFillColor(TColor::GetColor("#82FA58")); // light green
  sysenvelope->SetLineWidth(0);
  TBox *dsysenvelope = new TBox(xmin,1.0-2*MetPlotsSpace::Zprediction_Uncertainty,xmax,1.0+2*MetPlotsSpace::Zprediction_Uncertainty);
  dsysenvelope->SetFillColor(TColor::GetColor("#F3F781")); // light yellow
  dsysenvelope->SetLineWidth(0);
  
  MET_ratio->GetYaxis()->SetNdivisions(502,false);
  JZB_ratio->GetYaxis()->SetNdivisions(502,false);


  MET_observed->Draw("e1");
  MET_ttbar_pred->Draw("histo,same");
  ZpredClone->Draw("histo,same");
  MET_predicted->Draw("histo,same");
  MET_observed->Draw("e1,same");
  legend->Draw();
  if(is_data==data) DrawPrelim();
  else DrawMCPrelim();
  
  mainpad->Modified();
  main_canvas->cd();
  coverpad->Draw();
  coverpad->cd();
  coverpad->Range(0,0,1,1);
  coverpad->SetFillColor(kWhite);
  coverpad->SetBorderSize(0);
  coverpad->SetFrameFillColor(0);
  coverpad->Modified();
  main_canvas->cd();
  bottompad->SetTopMargin(ratiotopmargin);
  bottompad->SetBottomMargin(ratiobottommargin);
  bottompad->Draw();
  bottompad->cd();
  bottompad->Range(0,0,1,1);
  bottompad->SetFillColor(kWhite);
  
  MET_ratio->GetYaxis()->SetRangeUser(0,2);
  MET_ratio->GetXaxis()->SetLabelSize(xstretchfactor*MET_ratio->GetXaxis()->GetLabelSize());
  MET_ratio->GetYaxis()->SetLabelSize(xstretchfactor*MET_ratio->GetYaxis()->GetLabelSize());
  MET_ratio->GetXaxis()->SetTitleSize(xstretchfactor*gStyle->GetTitleSize());
  MET_ratio->GetYaxis()->SetTitleSize(xstretchfactor*gStyle->GetTitleSize());

  MET_ratio->Draw("e1");
//   dsysenvelope->Draw();
  sysenvelope->Draw();
  MET_ratio->Draw("AXIS,same");
  MET_ratio->Draw("e1,same");
  TLine *metl = new TLine(xmin,1.0,xmax,1.0);
  metl->SetLineColor(kBlue);
  metl->Draw();
  CompleteSave(main_canvas,SaveMETname.str());
  
  //--------------------------------------------------------------------------------------------
  mainpad->cd();
  
  JZB_observed->Draw("e1");
  JZB_ttbar_pred->Draw("histo,same");
  JZBZpredClone->Draw("histo,same");
  JZB_predicted->Draw("histo,same");
  JZB_observed->Draw("e1,same");
  legend->Draw();
  if(is_data==data) DrawPrelim();
  else DrawMCPrelim();
  
  main_canvas->cd();
  coverpad->Draw();
  main_canvas->cd();
  bottompad->Draw();
  bottompad->cd();
  JZB_ratio->GetYaxis()->SetRangeUser(0,2);
  
  JZB_ratio->GetXaxis()->SetLabelSize(xstretchfactor*JZB_ratio->GetXaxis()->GetLabelSize());
  JZB_ratio->GetYaxis()->SetLabelSize(xstretchfactor*JZB_ratio->GetYaxis()->GetLabelSize());
  JZB_ratio->GetXaxis()->SetTitleSize(xstretchfactor*gStyle->GetTitleSize());
  JZB_ratio->GetYaxis()->SetTitleSize(xstretchfactor*gStyle->GetTitleSize());

  JZB_ratio->Draw("e1");
//   dsysenvelope->Draw();
  sysenvelope->Draw();
  JZB_ratio->Draw("AXIS,same");
  JZB_ratio->Draw("e1,same");
  metl->Draw();
  
  CompleteSave(main_canvas,SaveJZBname.str());
  
  //--------------------------------------------------------------------------------------------
  mainpad->cd();
  
  TH1F *SystHisto = (TH1F*)MET_predicted->Clone("SystHisto");
  TGraphErrors *stat3jS = MakeErrorGraphSystematicAndStatistical(MET_ttbar_pred,ZpredClone,MET_predicted,SystHisto);
  MET_predicted->SetLineColor(TColor::GetColor("#cc0066"));
  ZpredClone->SetLineColor(TColor::GetColor("#006600"));
  ZpredClone->SetFillColor(TColor::GetColor("#006600"));
  ZpredClone->SetFillStyle(3002); // light dots, not crushing other information
  ZpredClone->SetLineStyle(2);
  
  TPad *kinpad = new TPad("kinpad","kinpad",0,0,1,1);
  kinpad->SetLogy(1);
  kinpad->cd();
  
  MET_observed->Draw("e1");
  stat3jS->Draw("2,same");
  MET_observed->Draw("e1,same");
  ZpredClone->Draw("histo,same");
  MET_predicted->Draw("histo,same");
  MET_observed->Draw("e1,same");
  
  TLegend *legend2 = make_legend();
  legend2->SetX1(0.52);
  if (isAachen) legend2->SetHeader("N_{j}#geq 2");
  else legend2->SetHeader("N_{j}#geq 3");
  legend2->AddEntry(MET_observed,"Data","PL");
  legend2->AddEntry(MET_predicted,"Total backgrounds","L");
  legend2->AddEntry(ZpredClone,"DY (JZB)","FL");
  legend2->AddEntry(stat3jS,"Total uncert.","F");

  legend2->Draw();
  if(is_data==data) DrawPrelim();
  else DrawMCPrelim();
  
  cout << "About to store syst plot ... " << endl;
  save_with_ratio_and_sys_band(0.0, MET_observed, MET_predicted, kinpad, (SaveMETname.str()+"__WithSys"), false, false, "data/pred",SystHisto );  
  
//  delete main_canvas;
  delete MET_observed;
  delete MET_predicted;
  //do NOT delete MET_Z_prediction (it's the return value)
  delete MET_osof_pred;
  delete MET_ossf_pred;
  delete MET_ttbar_pred;
  
  delete JZB_observed;
  delete JZB_predicted;
  delete JZB_osof_pred;
  delete JZB_ossf_pred;
  delete JZB_Z_prediction;
  delete JZB_ttbar_pred;

  return MET_Z_prediction;
}

float extract_correction(string jzbvariable) {
  int position = (int)jzbvariable.find("[1]");
  if(position==-1) return 0.0;
  string correction=jzbvariable.substr(position+3,jzbvariable.length()-position-3);
  position = (int)correction.find("*");
  if(position==-1) return 0.0;
  correction=correction.substr(0,position);
  float correctionvalue=atof(correction.c_str());
  assert(correctionvalue<1&&correctionvalue>0);
  return correctionvalue;
}
  
float Get_Met_Z_Prediction(TCut JetCut, float MetCut, int isdata, bool isDYonly, bool isAachen=false) {
  dout << "Going to compute Z region prediction for a MET cut at " << MetCut << " GeV" << endl;
  // Steps: 
  // 1) Get peak
  // 2) use the peak and pt correction for sample splitting
  // and for MET distribution shifting
  // 3) compute the estimate for MET>MetCut
  
  // do this for data if isdata==data, otherwise for MC (full closure if isDYonly==false, otherwise use only DY sample)
  
  // Step 0 : If we're dealing with DY, we need to make sure PURW is off!
//   string bkpcutweight = (const char*) cutWeight;
//   if(isdata==mc && isDYonly) cutWeight=TCut("1.0");
  
  // Step 1) Get peak
  float MCPeakNoPtCorr=0,MCPeakErrorNoPtCorr=0,DataPeakNoPtCorr=0,DataPeakErrorNoPtCorr=0,MCSigma=0,DataSigma=0;
  stringstream resultsNoPtCorr;
  stringstream NoPtCorrdatajzb;
  stringstream NoPtCorrmcjzb;
  
  if(isAachen) {
    //need to make sure that none of the typical basic cuts contain problematic selections!
    string Sleptoncut = (const char*) leptoncut;
    if((int)Sleptoncut.find("abs(eta1)<1.4")>-1) {
      write_error(__FUNCTION__,"You're trying to compute the Aachen estimate but are requiring abs(eta1)<1.4 ... please check your config.");
      assert((int)Sleptoncut.find("abs(eta1)<1.4")==-1);
      assert((int)Sleptoncut.find("abs(eta1)<1.4")==-1);
    }
  } else {
    string Sleptoncut = (const char*) leptoncut;
    if((int)Sleptoncut.find("abs(eta1)<2.4")>-1) {
      write_error(__FUNCTION__,"You're trying to compute the ETH estimate but you are using non-central leptons ... please check your config.");
      assert((int)Sleptoncut.find("abs(eta1)<2.4")==-1);
      assert((int)Sleptoncut.find("abs(eta2)<2.4")==-1);
    }
  }
      

  float Ptcorrection=0.0;
  
  if(isdata==data) Ptcorrection=extract_correction(PlottingSetup::jzbvariabledata);
  else Ptcorrection=extract_correction(PlottingSetup::jzbvariablemc);
  
  bool OverFlowStatus=addoverunderflowbins;
  
  find_peaks(MCPeakNoPtCorr,MCPeakErrorNoPtCorr, DataPeakNoPtCorr,DataPeakErrorNoPtCorr,resultsNoPtCorr,true,NoPtCorrdatajzb,NoPtCorrmcjzb,(const char*) JetCut, true);
  
  switch_overunderflow(OverFlowStatus);
  
  float PeakPosition=0.0; 
  string jzbvariable;
  if(isdata==data) {
    PeakPosition=DataPeakNoPtCorr;
    jzbvariable=jzbvariabledata;
    dout << "Found peak in data at " << DataPeakNoPtCorr << " +/- " << DataPeakErrorNoPtCorr << " ; will use this result (" << PeakPosition << ")" << endl;
  } else {
    PeakPosition=MCPeakNoPtCorr;
    jzbvariable=jzbvariablemc;
    dout << "Found peak in mc at " << MCPeakNoPtCorr << " +/- " << MCPeakErrorNoPtCorr << " ; will use this result (" << PeakPosition << ")" << endl;
  }
  
  // Step 2: Use peak for sample splitting and MET shifting
  string CorrectedMet="met[4]/MetFactor-"+any2string(Ptcorrection)+"*pt +"+any2string(abs(1.0*(PeakPosition)));
  if(2*(PeakPosition)<0) CorrectedMet="met[4]/MetFactor-"+any2string(Ptcorrection)+"*pt -"+any2string(abs(1.0*(PeakPosition)));
  
  stringstream sPositiveCut;
  if(PeakPosition>0) sPositiveCut << "((" << jzbvariable << "-" << PeakPosition << ")>0)";
  else sPositiveCut << "(  " << jzbvariable << "+" << abs(PeakPosition) << ")>0)";
  
  stringstream sNegativeCut;
  if(PeakPosition<0) sNegativeCut  << "((" << jzbvariable << "+" << abs(PeakPosition) << ")<0)";
  else sNegativeCut << "((  " << jzbvariable << "-" << abs(PeakPosition) << ")<0)";
  
  string ObservedMet="met[4]/MetFactor";
  
  stringstream JZBPosvar;
  JZBPosvar<<jzbvariable;
  if(PeakPosition>0) JZBPosvar << "-" << PeakPosition;
  else JZBPosvar << "+" << abs(PeakPosition);
  
  stringstream JZBNegvar;
  JZBNegvar<<"-(" << jzbvariable;
  if(PeakPosition>0) JZBNegvar << "-" << PeakPosition << ")";
  else JZBNegvar << "+" << abs(PeakPosition) << ")";
  
  
  // Step 3: Compute estimate
  TH1F *predicted = GetPredictedAndObservedMetShapes(JetCut, sPositiveCut.str(),sNegativeCut.str(),CorrectedMet,ObservedMet,JZBPosvar.str(),JZBNegvar.str(), MetCut, isdata, isDYonly, isAachen);

  float ZregionZestimate=0;
  for(int ibin=1;ibin<=(int)predicted->GetNbinsX();ibin++) {
    if(predicted->GetBinLowEdge(ibin)+predicted->GetBinWidth(ibin)>MetCut) {
      ZregionZestimate+=2*(predicted->GetBinContent(ibin));
    }
  }
  
  dout << "  Z region estimate in MET>" << MetCut << " for this sample: " << ZregionZestimate << endl;
  ofstream NiceSummary;
  NiceSummary.open("Summary.txt",ios::app);
  NiceSummary<< "  Z peak Z estimate in MET>" << MetCut << " for this sample: " << ZregionZestimate << endl;
  NiceSummary.close();

  if(isdata==data) {
    MetPlotsSpace::Zestimate__data=ZregionZestimate;
    MetPlotsSpace::Zestimate__data_stat=2*TMath::Sqrt(ZregionZestimate/2);
    MetPlotsSpace::Zestimate__data_sys=ZregionZestimate*MetPlotsSpace::Zprediction_Uncertainty;
  }
  if(isdata==mc && isDYonly) {
    MetPlotsSpace::Zestimate__dy=ZregionZestimate;
    MetPlotsSpace::Zestimate__dy_stat=2*TMath::Sqrt(ZregionZestimate/2);
    MetPlotsSpace::Zestimate__dy_sys=ZregionZestimate*MetPlotsSpace::Zprediction_Uncertainty;
  }
  if(isdata==mc && !isDYonly) {
    MetPlotsSpace::Zestimate__mc=ZregionZestimate;
    MetPlotsSpace::Zestimate__mc_stat=2*TMath::Sqrt(ZregionZestimate/2);
    MetPlotsSpace::Zestimate__mc_sys=ZregionZestimate*MetPlotsSpace::Zprediction_Uncertainty;
  }
  
  
//   if(isdata==mc && isDYonly) cutWeight=TCut(bkpcutweight.c_str());

  return ZregionZestimate;
}
  
void ProvideEEOverMMEstimate(TCut GeneralCut) {
  TCanvas *eemmcan = new TCanvas("eemmcan","eemmcan");
  TCut completecut = TCut(GeneralCut&&Restrmasscut);
  TH1F *eeh = allsamples.Draw("eeh",  "mll",1,70,120,"m_{ll}","events",completecut&&TCut("id1==0"),data,luminosity);
  TH1F *mmh = allsamples.Draw("mmh",  "mll",1,70,120,"m_{ll}","events",completecut&&TCut("id1==1"),data,luminosity);
  
  float Nee = eeh->Integral();
  float Nmm = mmh->Integral();
  
  dout << "Ratio R(ee/mm) = " << Nee/Nmm << " +/- " << sqrt(1/(Nee) + 1/(Nmm)) << endl;
  
  delete eemmcan;
  delete eeh;
  delete mmh;
}
  
void ExperimentalMetPrediction(bool QuickRun=false, bool isAachen=false) {
  
  switch_overunderflow(true);
  
  bool HighPurityMode=true; // High Purity = |mll-91|<10 GeV ,     else <20
  
  if(QuickRun) {
    HighPurityMode=true;
  }
  
  string restrmasscutbkp=(const char*) PlottingSetup::Restrmasscut;
  
  if(HighPurityMode) PlottingSetup::Restrmasscut=TCut("abs(mll-91)<10");
  else PlottingSetup::Restrmasscut= TCut("abs(mll-91)<20");
  
  dout << "Aachen mode (20/10, 2 jets) ? " << isAachen << endl;
  dout << "High Purity mode? " << HighPurityMode << endl;
  
  string backup_basicqualitycut = (const char*) basicqualitycut;
  string backup_essentialcut    = (const char*) essentialcut;
  string backup_basiccut        = (const char*) basiccut;
  string backup_leptoncut        = (const char*) leptoncut;
  
  if(isAachen) {
    string Sbasicqualitycut = backup_basicqualitycut;
    Sbasicqualitycut = ReplaceAll(Sbasicqualitycut,")<1.4",")<2.4");
    basicqualitycut=TCut(Sbasicqualitycut.c_str());
    
    string Sleptoncut = backup_leptoncut;
    Sleptoncut = ReplaceAll(Sleptoncut,")<1.4",")<2.4");
    leptoncut=TCut(Sleptoncut.c_str());
    
    string Sbasiccut = backup_basiccut;
    Sbasiccut = ReplaceAll(Sbasiccut,")<1.4",")<2.4");
    basiccut=TCut(Sbasiccut.c_str());
    
    string Sessentialcut = backup_essentialcut;
    Sessentialcut = ReplaceAll(Sessentialcut,")<1.4",")<2.4");
    essentialcut=TCut(Sessentialcut.c_str());
    
    cout << "Basic cut : " << (const char*) basiccut << endl;
    cout << "Essential cut : " << (const char*) essentialcut << endl;
  }
  
  
  stringstream snjets;
  if(isAachen) snjets << 2;
  else snjets << 3;
  float maxMET=100;
  if(isAachen) maxMET=150;
    
  TCut nJetsSignal(PlottingSetup::basicqualitycut&&("pfJetGoodNum40>="+snjets.str()).c_str());
  
  
  dout << " ***** TESTING Z PREDICTION ***** " << endl;
  dout << "Notation (you can copy & paste this to evaluate it further)" << endl;
  dout << "Cut;Data;MC;DY;" << endl;
  float DataEstimate = -1;
  DataEstimate = Get_Met_Z_Prediction(Restrmasscut&&nJetsSignal,maxMET, data, false, isAachen);
  float DYEstimate=-1;
  if(!QuickRun) DYEstimate = Get_Met_Z_Prediction(Restrmasscut&&nJetsSignal,maxMET, mc, true, isAachen);
  float MCEstimate=-1;
  if(!QuickRun) MCEstimate = Get_Met_Z_Prediction(Restrmasscut&&nJetsSignal,maxMET, mc, false, isAachen);
  
  dout << "Z prediction (JZB based) " << DataEstimate << endl;
  write_info(__FUNCTION__,"Z prediction (JZB based) "+any2string(DataEstimate));
  if(QuickRun) return;
  dout << maxMET << ";" << DataEstimate << ";" << MCEstimate << ";" << DYEstimate << endl;
  
  float Diff=20.0;
  if(HighPurityMode) Diff=10;
  TCut cut("mll>20&&pt1>20&&pt2>20");
  if (isAachen) cut = TCut("mll>20&&pt1>20&&pt2>20&&pfTightHT>100");
  
  TCanvas *qcan = new TCanvas("qcan","qcan");
  TH1F *zlineshape = allsamples.Draw("zlineshape","mll",int((91+25-18)*5),18,91+25,"m_{ll} (GeV)","events",cutOSSF&&TCut("pfJetGoodNum40==2")&&cut,data,PlottingSetup::luminosity); // bins of 0.2 GeV
  TH1F *Ozlineshape = allsamples.Draw("Ozlineshape","mll",int((91+25-18)*5),18,91+25,"m_{ll} (GeV)","events",cutOSOF&&TCut("pfJetGoodNum40==2")&&cut,data,PlottingSetup::luminosity); // bins of 0.2 GeV
  zlineshape->Add(Ozlineshape,-1);
  delete qcan;
  float a = (zlineshape->Integral(zlineshape->FindBin(20),zlineshape->FindBin(70)));
  float b = (zlineshape->Integral(zlineshape->FindBin(91-Diff),zlineshape->FindBin(91+Diff)));
  float r = a/b;
  float dr= (a/b)*TMath::Sqrt(1/a+1/b);
  
  float SysUncertainty  = TMath::Sqrt(DataEstimate*DataEstimate*dr*dr + r*r*(DataEstimate*MetPlotsSpace::Zprediction_Uncertainty*DataEstimate*MetPlotsSpace::Zprediction_Uncertainty));
  float StatUncertainty = TMath::Sqrt(DataEstimate);
  
  dout << "Z estimate in peak : " << DataEstimate   << " +/- " << DataEstimate*MetPlotsSpace::Zprediction_Uncertainty << " (sys) +/- " << TMath::Sqrt(2*DataEstimate) << " (stat) " << endl;
  dout << "Z ESTIMATE IN SR   : " << DataEstimate*r << " +/- " << SysUncertainty                                      << " (sys) +/- " << StatUncertainty           << " (stat) " << endl;
//   cout << endl;
  dout << "r = " << r << " +/- " << dr << endl;
  
  
  delete Ozlineshape;
  delete zlineshape;
  
  PlottingSetup::Restrmasscut=TCut(restrmasscutbkp.c_str());
  
  basicqualitycut=TCut(backup_basicqualitycut.c_str());
  basiccut       =TCut(backup_basiccut.c_str());
  essentialcut   =TCut(backup_essentialcut.c_str());
  leptoncut      =TCut(backup_leptoncut.c_str());

  switch_overunderflow(false);

}


void ComputeRMuE(string basecut, string name) {
  
  
  dout << "Starting with computation of rmue, but need to reset the essential cut first to make sure no restriction on eta is in there" << endl;
  
  string backup_essentialcut = (const char*) essentialcut;
  string NewEssentialCut = backup_essentialcut;
  NewEssentialCut = ReplaceAll(NewEssentialCut,")<1.4",")<2.4");
  essentialcut = TCut(NewEssentialCut.c_str());
  
  TCut kbase((basecut+" && mll>60 && mll<120").c_str());
  TCut ee("id1==id2&&id1==0");
  TCut mm("id1==id2&&id1==1");
  
  TH1F *Dee   = allsamples.Draw("Dee",  "mll",100,60,120,"m_{ll}","events",kbase&&ee,data,luminosity);
  TH1F *Dmm   = allsamples.Draw("Dmm",  "mll",100,60,120,"m_{ll}","events",kbase&&mm,data,luminosity);
  TH1F *Mee   = allsamples.Draw("Mee",  "mll",100,60,120,"m_{ll}","events",kbase&&ee,mc,luminosity);
  TH1F *Mmm   = allsamples.Draw("Mmm",  "mll",100,60,120,"m_{ll}","events",kbase&&mm,mc,luminosity);
  
  
  float mc_rmue=sqrt(Mmm->Integral()/Mee->Integral());
  float d_mc_rmue   = sqrt(1/(4*Mee->Integral())+Mmm->Integral()/(4*Mee->Integral()*Mee->Integral()));
  float data_rmue=sqrt(Dmm->Integral()/Dee->Integral());
  float d_data_rmue = sqrt(1/(4*Dee->Integral())+Dmm->Integral()/(4*Dee->Integral()*Dee->Integral()));
  
  dout << "Rmu table (" << name << ") : " << endl;
  dout << "   \t n_{\\mu\\mu} \t n_{ee} \t r_{\\mu e} \\pm stat" << endl;
  dout << fixed << "MC \t " << Mmm->Integral() << " \t " << Mee->Integral() << "\t" <<   mc_rmue << " \\pm " <<   d_mc_rmue << endl;
  dout << fixed << "Data\t" << Dmm->Integral() << " \t " << Dee->Integral() << "\t" << data_rmue << " \\pm " << d_data_rmue << endl;
  
  delete Dee;
  delete Dmm;
  delete Mee;
  delete Mmm;
  
  essentialcut   =TCut(backup_essentialcut.c_str());
}


void ComputeRMuE() {
  TCanvas *rmucanvas = new TCanvas("rmucanvas","rmucanvas");
  dout << "Starting with our selection (ETH/Tight) : " << endl;
  ComputeRMuE("abs(eta1)<1.4 && abs(eta2)<1.4 && pt1>20 && pt2>20"," Tight selection");
  dout << "And now the Aachen selection (Aachen/Loose) : " << endl;
  ComputeRMuE("abs(eta1)<2.4 && abs(eta2)<2.4 && pt1>20 && pt2>20"," Loose selection");
  delete rmucanvas;
}  

Revision:	1.52
Committed:	Fri Apr 26 06:44:13 2013 UTC (12 years ago) by buchmann
Content type:	text/plain
Branch:	MAIN
Changes since 1.51:	+161 -97 lines
Log Message:	Added MC yields, made sure everything starts at 20 GeV in the dileptonic mass distribution (no underflow from 15); added cross-checks; added much more information about control regions 1&2, also compatibility for Aachen; adapted code not to use smeared DY for MC yields (in order to be compatible with Aachen)