#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <vector>
#include <string.h>

#include <TROOT.h>
#include <TSystem.h>
#include <TFile.h>
#include <TChain.h>
#include <TString.h>
#include <TH1.h>
#include <THStack.h>
#include <TMinuit.h>
#include <TStyle.h>
#include <TLegend.h>
#include <TCanvas.h>
#include <TPaveText.h>
#include <TColor.h>
#include <TLatex.h>

#include <TRandom.h>
#include <TBenchmark.h>
#include <RooRandom.h>

#include "RooFit.h"
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooExponential.h"
#include "RooChebychev.h"
#include "RooBreitWigner.h"
#include "RooRealVar.h"
#include "RooAbsPdf.h"
#include "RooWorkspace.h"
#include "RooDataSet.h"
#include "RooSimultaneous.h"
#include "RooArgSet.h"
#include "RooCategory.h"
#include "RooFitResult.h"
#include "RooPlot.h"
#include "RooAbsCategoryLValue.h"
#include "RooStats/RooStatsUtils.h"
#include "RooStats/ModelConfig.h"

#include "RooCurve.h"
#include "RooFitResult.h"

using namespace std;
using namespace RooFit;
using namespace RooStats;

#define GeV   1e3

void postfit(TString channel=""){

  TString CHANNEL;
  TString FNAME;
  double MU_VALUE;
  int NBINS;

  double *the_bins;

  if (channel == "tag_el") {
    CHANNEL="tag_el";
    FNAME="lephad_tag_el_mA200tb00_ggA_comb_combined_AllSyst_model.root";
    MU_VALUE = 0.0;
    NBINS = 22;
    the_bins = new double[23];
    double the_bins_[23] = { 0,  50, 70, 80, 90, 95, 100,  105,  110,  115,  120, 125, 130, 135, 140, 150, 160, 170, 180, 200, 250, 320, 420 };
    for (int i=0;i<=NBINS; ++i) { the_bins[i] = the_bins_[i]; }
    
  }
  else if (channel == "tag_mu") {    
    CHANNEL="tag_mu";
    FNAME="lephad_tag_mu_mA200tb00_ggA_comb_combined_AllSyst_model.root";
    MU_VALUE = 0.0;
    NBINS = 22;
    the_bins = new double[23];
    double the_bins_[23] = { 0,  50, 70, 80, 90, 95, 100,  105,  110,  115,  120, 125, 130, 135, 140, 150, 160, 170, 180, 200, 250, 320, 420 };
    for (int i=0;i<=NBINS; ++i) { the_bins[i] = the_bins_[i]; }
  }
  else if (channel == "high_el") {    
    CHANNEL="high_el";
    FNAME="lephad_high_el_mA400tb00_ggA_comb_combined_AllSyst_model.root";
    MU_VALUE = 0.0;
    NBINS = 36;
    the_bins = new double[37];
    double the_bins_[37] =     {0, 160, 172, 180, 186, 194, 200, 206, 212, 218, 224,230,236,242,248,256,264,274,280,288,296,306,314,324,332,348,362,374,392,408,440,468,492,542,630,776,886};
    for (int i=0;i<=NBINS; ++i) { the_bins[i] = the_bins_[i]; }
  }
  else if (channel == "veto_el") {
    CHANNEL= "veto_el";
    FNAME ="lephad_veto_el_mA200tb00_ggA_comb_combined_AllSyst_model.root";
    MU_VALUE = 0.0;
    NBINS = 31;
    the_bins = new double[32];
    double the_bins_[32] = {0,  20, 40, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100,
			    105, 110, 115, 120, 125, 130, 135, 140, 145,
			    150, 160, 170, 180, 200, 225, 250, 300, 400};
    for (int i=0;i<=NBINS; ++i) { the_bins[i] = the_bins_[i]; }
    
  }
  else { // veto_mu is the default value
    CHANNEL= "veto_mu";
    FNAME ="lephad_veto_mu_mA200tb00_ggA_comb_combined_AllSyst_model.root";
    MU_VALUE = 0.0;
    NBINS = 31;
    the_bins = new double[32];
    double the_bins_[32] = {0,  20, 40, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100,
			    105, 110, 115, 120, 125, 130, 135, 140, 145,
			    150, 160, 170, 180, 200, 225, 250, 300, 400};
    for (int i=0;i<=NBINS; ++i) { the_bins[i] = the_bins_[i]; }
  }

//   TString CHANNEL("veto_mu");
//   TString FNAME("lephad_veto_mu_mA200tb00_ggA_comb_combined_AllSyst_model.root");

//   double MU_VALUE = 0.0;
//   const int NBINS = 31;
//   double the_bins[32] = {0,  20, 40, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95,   
// 100,   105,   110,   115,   120,   125,   130,   135,   140,   145,   150,   
// 160,   170,   180,   200,   225,   250,   300, 400};

  /*
  TString CHANNEL("tag_el");
  TString FNAME("lephad_tag_el_mA200tb00_ggA_comb_combined_AllSyst_model.root");

  double MU_VALUE = 0.0;
  const int NBINS = 22;

  double the_bins[23] = { 0,  50, 70, 80, 90, 95, 100,  105,  110,  115,  120, 125, 130, 135, 140, 150, 160, 170, 180, 200, 250, 320, 420 };

   */

  /*
  TString CHANNEL("tag_mu");
  TString FNAME("lephad_tag_mu_mA200tb00_ggA_comb_combined_AllSyst_model.root");

  double MU_VALUE = 0.0;
  const int NBINS = 22;

  double the_bins[23] = { 0,  50, 70, 80, 90, 95, 100,  105,  110,  115,  120, 125, 130, 135, 140, 150, 160, 170, 180, 200, 250, 320, 420 };

   */


  TFile* fin = TFile::Open(FNAME,"read");

  RooWorkspace* ws = (RooWorkspace*) fin->Get("combined");
  RooDataSet* data = (RooDataSet*)   ws->data("obsData");
  RooRealVar* mu   = (RooRealVar*)   ws->var("SigXsecOverSM");
  //  RooRealVar* met  = (RooRealVar*)   ws->var("weightVar");
  // rQCD_veto_mu, rQCD_veto_el, rQCD_tag_mu, rQCD_tag_el, rQCD_high_el
  RooRealVar* rQCD   = (RooRealVar*)   ws->var("rQCD_"+CHANNEL);

  mu ->setVal(MU_VALUE);
  mu ->setConstant(true); 

  double rQCD_prefit = rQCD->getVal();
  
  ModelConfig* mc          = (ModelConfig*)  ws->obj("ModelConfig");
  RooSimultaneous* simPdf  = dynamic_cast<RooSimultaneous*>(mc->GetPdf());
  //RooArgSet* observables   = const_cast<RooArgSet*>(mc->GetObservables());
  //RooArgSet* globalObs     = const_cast<RooArgSet*>(mc->GetGlobalObservables());
  //RooArgSet* nuis          = const_cast<RooArgSet*>(mc->GetNuisanceParameters());
  //RooArgSet* pdfConstraint = const_cast<RooArgSet*>(simPdf->getAllConstraints(*observables, *nuis, false));

  // Fitting
  RooArgSet constrainParams;
  if (mc->GetNuisanceParameters() ) constrainParams.add(*mc->GetNuisanceParameters());
  RooStats::RemoveConstantParameters(&constrainParams);
  //RooStats::SetAllConstant(*nuis);
  
  RooFitResult * fitres   =   mc->GetPdf()->fitTo(*data,InitialHesse(false),Hesse(false),"Migrad",Strategy(0),0,Constrain(constrainParams),Save(true));

  // now that you have done the fit, get the rQCD value
  double rQCD_postfit = rQCD->getVal();
 
  cout << "CHECK: " << rQCD_prefit << " vs " << rQCD_postfit << endl;
  //  RooSimultaneous* simPdf = dynamic_cast<RooSimultaneous*>(mc->GetPdf());
  simPdf = dynamic_cast<RooSimultaneous*>(mc->GetPdf());
  
  //  RooCategory category  = *(dynamic_cast<RooCategory*> (simPdf->indexCat()));
  //  TList* datasets = data->split(category, true);
  TList* datasets = data->split(simPdf->indexCat(), true);
  TIter next(datasets);
  int pdf_index = 0;
  vector<TString> bkgNames;
  cout << "TESTING PREFIT" << endl;
  // Prefit
  TFile* fin_p = TFile::Open(FNAME,"read");
  RooWorkspace* ws_p = (RooWorkspace*) fin_p->Get("combined");
  ModelConfig* mc_p  = (ModelConfig*)  ws_p->obj("ModelConfig");
  RooSimultaneous* simPdf_p  = dynamic_cast<RooSimultaneous*>(mc_p->GetPdf());
  RooRealVar* mu_p   = (RooRealVar*)  ws_p->var("SigXsecOverSM");
  mu_p ->setVal(MU_VALUE);
  mu_p ->setConstant(true); 

  RooRealVar* rQCD_p   = (RooRealVar*)  ws_p->var("rQCD_"+CHANNEL);
  rQCD_p->setVal(rQCD_postfit);
  rQCD_p->setConstant(true);
  // you know need to set the rQCD to the post-fit value for a correct
  // QCD comparison
  

  //
  cout << "ENDING PREFIT" << endl;

  Int_t colort[6];
  bkgNames.push_back("ATLAS_Bkg_LH_qcd_cC");      colort[0]=6;
  bkgNames.push_back("ATLAS_Bkg_LH_Other");    colort[1]=7;
  bkgNames.push_back("ATLAS_Bkg_LH_TT");      colort[2]=5;
  bkgNames.push_back("ATLAS_Bkg_LH_Wlnu");    colort[3]=8;
  bkgNames.push_back("ATLAS_Bkg_LH_Zleplep");    colort[4]=9;
  bkgNames.push_back("ATLAS_Bkg_LH_Ztautau");  colort[5]=10;

  // make all these normal Root histograms
  TH1F * h_dat = new TH1F("h_dat", "dat", NBINS, the_bins);

  // post fit plots

  TH1F * h_all = new TH1F("h_all", "all", NBINS, the_bins);
  TH1F * h_qcd = new TH1F("h_qcd", "qcd", NBINS, the_bins);
  TH1F * h_oth = new TH1F("h_oth", "oth", NBINS, the_bins);
  TH1F * h_top = new TH1F("h_top", "top", NBINS, the_bins);
  TH1F * h_wln = new TH1F("h_wln", "wln", NBINS, the_bins);
  TH1F * h_zll = new TH1F("h_zll", "zll", NBINS, the_bins);
  TH1F * h_ztt = new TH1F("h_ztt", "ztt", NBINS, the_bins);

  TH1F * h_all_up = new TH1F("h_all_up", "all", NBINS, the_bins);
  TH1F * h_all_do = new TH1F("h_all_do", "all", NBINS, the_bins);
  // pre fit plots

  TH1F * h_all_pre = new TH1F("h_all_pre", "all", NBINS, the_bins);
  TH1F * h_qcd_pre = new TH1F("h_qcd_pre", "qcd", NBINS, the_bins);
  TH1F * h_oth_pre = new TH1F("h_oth_pre", "oth", NBINS, the_bins);
  TH1F * h_top_pre = new TH1F("h_top_pre", "top", NBINS, the_bins);
  TH1F * h_wln_pre = new TH1F("h_wln_pre", "wln", NBINS, the_bins);
  TH1F * h_zll_pre = new TH1F("h_zll_pre", "zll", NBINS, the_bins);
  TH1F * h_ztt_pre = new TH1F("h_ztt_pre", "ztt", NBINS, the_bins);

  TH1F * h_all_pre_up = new TH1F("h_all_pre_up", "all", NBINS, the_bins);
  TH1F * h_all_pre_do = new TH1F("h_all_pre_do", "all", NBINS, the_bins);

  RooPlot* xframe[6];

  double background[6]; double background_p[6]; double datan[6];

  TString f_name("yields_"+CHANNEL+"_post.txt");
  TString f_p_name("yields_"+CHANNEL+"_pre.txt");
  TString f_comp_name("yields_"+CHANNEL+"_Comparison.txt");
  ofstream f;     f.open(f_name.Data());
  ofstream f_p;   f_p.open(f_p_name.Data());
  ofstream f_comp;   f_comp.open(f_comp_name.Data());

  f_comp << "Worked on "  << FNAME << endl;
  f_comp << "Assuming mu= "  << MU_VALUE << endl;
  f_comp << "Fixing rQCD in pre-fit to  " << rQCD_postfit << endl;

  double pre_fit_background_in_bins[NBINS+2];
  double post_fit_background_in_bins[NBINS+2];
  double data_total_in_bins[NBINS+2];
  for (int ii=0;ii<NBINS+2; ++ii) {
    pre_fit_background_in_bins[ii] = 0;
    post_fit_background_in_bins[ii] = 0;
    data_total_in_bins[ii] =0;
  }

  for(RooAbsData* ds=(RooAbsData*) next(); ds != 0; ds=(RooAbsData*) next()){

    cout<<"ds->GetName(): "<< ds->GetName() <<"   "<<endl;
    
    TString obsname=ds->GetName();
    if (obsname.Contains("_num_") || obsname.Contains("_den_")) {
      cout << "We won't work on " << obsname << endl;
      continue;
    }
    RooRealVar *x = (RooRealVar*) ws->obj("obs_x_"+obsname); 
    x->setRange(0.,NBINS);
    //
    RooAbsPdf* thisPdf        = simPdf  ->getPdf(ds->GetName());

    RooArgSet* obstmp   = thisPdf   ->getObservables(*mc  ->GetObservables()) ;
    background   [pdf_index] = thisPdf  ->expectedEvents(*obstmp ); 

    datan      [pdf_index] = ds->numEntries();   

    xframe[pdf_index] = x->frame();
    xframe[pdf_index]->SetMinimum(1e-5);
    data          ->plotOn(xframe[pdf_index]);

    // RooCurve* h_data = xframe[pdf_index]->getCurve(); 
    // if (h_data) cout << "yes" << endl;
    // else cout << "NO!!" << endl;
    data->Print("all");
    //    RooAbsRealILValue the_x = data->var("obs_x_lephad_veto_mu");
    TH1F *h_data_1 = (TH1F*) ds->createHistogram("obs_x_lephad_"+CHANNEL, NBINS);
    
    cout << "total data: " << ds->sumEntries() <<endl;
    cout << h_data_1->Integral() << endl;
    double data_total=0;
    for(Int_t binid=0;binid<=NBINS+1;binid++) {
      //f<<"bin:"<<binid<<" "<<h_data->interpolate(binid)<<endl;
      data_total += h_data_1->GetBinContent(binid);
      data_total_in_bins[binid] = h_data_1->GetBinContent(binid);
      h_dat->SetBinContent( binid, h_data_1->GetBinContent(binid)  );
      //      h_dat->SetBinError( binid, h_data_1->GetBinError(binid)  );
      h_dat->SetBinError( binid, sqrt(h_data_1->GetBinContent(binid))  );
    }
    cout << "CHECK: " << data_total << " vs " << ds->sumEntries() << endl;
    thisPdf   ->plotOn(xframe[pdf_index],Range(0.,NBINS),LineColor(1),Normalization(background  [pdf_index],datan[pdf_index]),LineStyle(kDashed));
    
    //f<<(pdf_index==0?"2mu2nu 2011":(pdf_index==1?"2e2nu 2011":(pdf_index==2?"2mu2nu 2012":(pdf_index==3?"2e2nu 2012":0))))<<" "<<endl;    
    RooCurve* h = xframe[pdf_index]->getCurve(); 
    double bkg_total_post_fit=0;
    for(Int_t binid=0;binid<NBINS;binid++) {
      f<<"bin:"<<binid<<" "<<h->interpolate(binid+0.5)<<endl;
      bkg_total_post_fit += h->interpolate(binid+0.5);
      h_all->SetBinContent(binid+1, h->interpolate(binid+0.5));
      //      h_all_up->SetBinContent(binid, h->interpolate(binid+0.5));
      cout << "test: " << binid<< " : " << h->interpolate(binid+0.5) << endl;
    }

    thisPdf   ->plotOn(xframe[pdf_index],Range(0.,NBINS),Normalization(background  [pdf_index],datan[pdf_index]),VisualizeError(*fitres,1),FillColor(kOrange));
    h = xframe[pdf_index]->getCurve(); 
    for(Int_t binid=0;binid<NBINS;binid++) {
      cout << "test error: " << binid+1<< " : " << h->interpolate(binid+0.5) << endl;
      h_all_up->SetBinContent(binid+1, h->interpolate(binid+0.5));
      double diff = (h_all->GetBinContent(binid+1)>0)?h_all->GetBinContent(binid+1)-(fabs(h_all->GetBinContent(binid+1) - h_all_up->GetBinContent(binid+1)))/h_all->GetBinContent(binid+1):0;
      h_all_do->SetBinContent(binid+1, diff);
      cout << "test error: " << binid+1<< " : " << h->interpolate(binid+0.5) << endl;
    }


    //  
    RooAbsPdf* thisPdf_p      = simPdf_p->getPdf(ds->GetName());

    RooArgSet* obstmp_p = thisPdf_p ->getObservables(*mc_p->GetObservables()) ;
    background_p [pdf_index] = thisPdf_p->expectedEvents(*obstmp_p); 

    thisPdf_p ->plotOn(xframe[pdf_index],Range(0.,NBINS),LineColor(1),Normalization(background_p[pdf_index],datan[pdf_index]));
    
    //    f_p<<(pdf_index==0?"2mu2nu 2011":(pdf_index==1?"2e2nu 2011":(pdf_index==2?"2mu2nu 2012":(pdf_index==3?"2e2nu 2012":0))))<<" "<<endl;    
    //    RooCurve* h = xframe[pdf_index]->getCurve(); 
    h = xframe[pdf_index]->getCurve(); 
    Double_t totb=0;
    for(Int_t binid=0;binid<NBINS;binid++){totb+=h->interpolate(binid+0.5); /*f_p<<"bin:"<<binid<<" "<<h->interpolate((binid-1)*30+15+90)<<endl;*/ }
    f_p<<"TotalBkg Prefit: " <<totb << " Postfit: "<< bkg_total_post_fit  <<endl;
    printf("%20s %10.4f \n", obsname.Data(),totb);
    //

    ///break down the background  
    for(UInt_t i = 0; i < bkgNames.size(); i++){
      TString bkgTmpName(thisPdf ->GetName());
      // if(bkgTmpName.Contains("2011"))
      //   bkgTmpName.ReplaceAll("model",bkgNames.at(i)+"_2011");
      // if(bkgTmpName.Contains("2012"))
      //   bkgTmpName.ReplaceAll("model",bkgNames.at(i)+"_2012");
      bkgTmpName = bkgTmpName+"_Hist_alphanominal";
      bkgTmpName = bkgNames.at(i)+"_lephad_"+CHANNEL+"_Hist_alphanominal";
      if (CHANNEL.Contains("veto_")) {
	if ( (bkgNames.at(i)).Contains("_TT") ||
	     (bkgNames.at(i)).Contains("_Other") ) {
	  bkgTmpName = bkgNames.at(i)+"_lephad_"+CHANNEL+"_nominal";
	}
      }
      if (CHANNEL.Contains("tag_el")) {
	if ( (bkgNames.at(i)).Contains("_Zleplep") ||
	     (bkgNames.at(i)).Contains("_Wlnu") ||
	     (bkgNames.at(i)).Contains("_Wlnu") ||
	     (bkgNames.at(i)).Contains("_Ztautau") ||
	     (bkgNames.at(i)).Contains("_TT") ||
	     (bkgNames.at(i)).Contains("_Other") ||
	     (bkgNames.at(i)).Contains("_qcd")
	     )
	bkgTmpName = bkgNames.at(i)+"_lephad_"+CHANNEL+"_nominal";
      }
      if (CHANNEL.Contains("tag_mu")) {
	if ( 
	     (bkgNames.at(i)).Contains("_Wlnu") ||
	     (bkgNames.at(i)).Contains("_Wlnu") ||
	     (bkgNames.at(i)).Contains("_Ztautau") ||
	     (bkgNames.at(i)).Contains("_TT") ||
	     (bkgNames.at(i)).Contains("_Other") ||
	     (bkgNames.at(i)).Contains("_qcd")
	     )
	bkgTmpName = bkgNames.at(i)+"_lephad_"+CHANNEL+"_nominal";
      }
      if (CHANNEL.Contains("high_")) {
	if (  (!(bkgNames.at(i)).Contains("_qcd")) 
	      && (!(bkgNames.at(i)).Contains("_Zleplep")) 
	      )
	  bkgTmpName = bkgNames.at(i)+"_lephad_"+CHANNEL+"_nominal";
      }


      cout << "Found " << bkgTmpName << endl;

      thisPdf    ->plotOn(xframe[pdf_index],Range(0.,NBINS),LineWidth(1),LineStyle(kDashed),LineColor(colort[i]),Normalization(background[pdf_index],datan[pdf_index]),Components(bkgTmpName));
      cout << bkgNames.at(i) << endl;
      f<<bkgNames.at(i)<<endl;
      //      RooCurve* h = xframe[pdf_index]->getCurve(); 
      h = xframe[pdf_index]->getCurve(); 
      double bkg_post_fit=0;
      vector<double> v_post_fit_bins;
      for(Int_t binid=0;binid<NBINS;binid++) {
	bkg_post_fit += h->interpolate(binid+0.5);
	f<<"bin:"<<binid<<" "<<h->interpolate(binid+0.5)<<endl;
	v_post_fit_bins.push_back(h->interpolate(binid+0.5));
	post_fit_background_in_bins[binid+1] += h->interpolate(binid+0.5);
	cout << bkgNames.at(i) << " " << (binid+1)  << " : "
	     << h->interpolate(binid+0.5) 
	     << endl;
      }
      f<<"Total: "<< bkg_post_fit <<endl;
      thisPdf_p  ->plotOn(xframe[pdf_index],Range(0.,NBINS),LineWidth(1),LineColor(colort[i]),Normalization(background_p[pdf_index],datan[pdf_index]),Components(bkgTmpName));
      
      //f_p<<bkgNames.at(i)<<endl;
      //RooCurve* h = xframe[pdf_index]->getCurve(); 
      h = xframe[pdf_index]->getCurve(); 
      Double_t bkg=0;
      vector<double> v_pre_fit_bins;
      for(Int_t binid=0;binid<NBINS;binid++){ 
	bkg+=h->interpolate(binid+0.5);  
	v_pre_fit_bins.push_back(h->interpolate(binid+0.5));
	pre_fit_background_in_bins[binid+1] += h->interpolate(binid+0.5);
/*f_p<<"bin:"<<binid<<" "<<h->interpolate((binid-1)*30+15+90)<<endl;*/ 
      }
      //f_p<<bkg<<endl;
      f_p<<bkgNames.at(i)<< ": " <<bkg<<" Postfit: "<< bkg_post_fit <<endl;
      printf("%25s %10.4f\n", bkgTmpName.Data(),bkg);
      double b_sum_pre=0; double b_sum_post=0;
      for (int jj=1; jj<=NBINS; ++jj) {
	char buffer[100];
	sprintf(buffer, "%25s %2i : %10.3f  %10.3f     %10.3f",  
		(bkgNames.at(i)).Data(), jj, 
		v_pre_fit_bins[jj-1], v_post_fit_bins[jj-1], 
		( (v_pre_fit_bins[jj-1]>0)?v_post_fit_bins[jj-1]/v_pre_fit_bins[jj-1]:0));
	f_comp << buffer <<endl;
	b_sum_pre += v_pre_fit_bins[jj-1];
	b_sum_post += v_post_fit_bins[jj-1];
	// fill the post fit histograms
	if ((bkgNames.at(i)).Contains("_qcd")) {
	  h_qcd->SetBinContent(jj, v_post_fit_bins[jj-1]);
	  h_qcd->SetBinError(jj, 0);
	  h_qcd_pre->SetBinContent(jj, v_pre_fit_bins[jj-1]);
	  h_qcd_pre->SetBinError(jj, 0);
	} else if ((bkgNames.at(i)).Contains("_TT")) {
	  h_top->SetBinContent(jj, v_post_fit_bins[jj-1]);
	  h_top->SetBinError(jj, 0);
	  h_top_pre->SetBinContent(jj, v_pre_fit_bins[jj-1]);
	  h_top_pre->SetBinError(jj, 0);
	} else if ((bkgNames.at(i)).Contains("_Ztautau")) {
	  h_ztt->SetBinContent(jj, v_post_fit_bins[jj-1]);
	  h_ztt->SetBinError(jj, 0);
	  h_ztt_pre->SetBinContent(jj, v_pre_fit_bins[jj-1]);
	  h_ztt_pre->SetBinError(jj, 0);
	} else if ((bkgNames.at(i)).Contains("_Zleplep")) {
	  h_zll->SetBinContent(jj, v_post_fit_bins[jj-1]);
	  h_zll->SetBinError(jj, 0);
	  h_zll_pre->SetBinContent(jj, v_pre_fit_bins[jj-1]);
	  h_zll_pre->SetBinError(jj, 0);
	} else if ((bkgNames.at(i)).Contains("_Wlnu")) {
	  h_wln->SetBinContent(jj, v_post_fit_bins[jj-1]);
	  h_wln->SetBinError(jj, 0);
	  h_wln_pre->SetBinContent(jj, v_pre_fit_bins[jj-1]);
	  h_wln_pre->SetBinError(jj, 0);
	} else if ((bkgNames.at(i)).Contains("_Other")) {
	  h_oth->SetBinContent(jj, v_post_fit_bins[jj-1]);
	  h_oth->SetBinError(jj, 0);
	  h_oth_pre->SetBinContent(jj, v_pre_fit_bins[jj-1]);
	  h_oth_pre->SetBinError(jj, 0);
	}
      }
      char buffer[100];
      sprintf(buffer, "%25s SUM: %10.3f  %10.3f     %10.3f",  
	      (bkgNames.at(i)).Data(), 
	      b_sum_pre, b_sum_post, b_sum_post/b_sum_pre);
      f_comp << buffer <<endl;

      f_comp << "-------------------------------" <<endl;
    }
    pdf_index ++;

    ///////////
    f_comp << "Total sum here:       Prefit      Postfit      Postfit/Prefit"<<endl;
    double sum_pre =0;  double sum_post=0; double sum_data=0;
    for (int jj=1;jj<=NBINS;++jj) {
 	char buffer[100];
	sprintf(buffer, "%25s %2i : %10.3f  %10.3f     %10.3f   %10.3f  %10.3f",  
		"Total", jj, 
		pre_fit_background_in_bins[jj], 
		post_fit_background_in_bins[jj], 
		( (pre_fit_background_in_bins[jj]>0)?post_fit_background_in_bins[jj]/pre_fit_background_in_bins[jj]:0   ),
		data_total_in_bins[jj], 
		( (post_fit_background_in_bins[jj]>0)?data_total_in_bins[jj]/post_fit_background_in_bins[jj]:0) 
		);
	f_comp << buffer <<endl;
	sum_pre += pre_fit_background_in_bins[jj];
	sum_post += post_fit_background_in_bins[jj];
	sum_data += data_total_in_bins[jj];
    }
    char buffer[100];
    sprintf(buffer, "%25s    : %10.3f  %10.3f     %10.3f   %10.3f  %10.3f",
	    "Total SUM", sum_pre, sum_post, sum_post/sum_pre, sum_data,
	    sum_data/sum_post);
    f_comp << buffer <<endl;
  }

  cout<<" Ready to plot "<<endl;
  
  // -----  ATLAS FORMAT -----  
  gROOT->ProcessLine(".L ~/atlasStyle.C");
  gROOT->ProcessLine("setAtlasStyle()");

  // Int_t icol=0;  Int_t font=42;  Double_t tsize=0.05;
  // gStyle->SetOptStat(0);
  // gStyle->SetPalette();
  // gStyle->SetPalette();
  // gStyle->SetFrameBorderMode(icol);
  // gStyle->SetFrameFillColor(icol);
  // gStyle->SetCanvasBorderMode(icol);
  // gStyle->SetPadBorderMode(icol);
  // gStyle->SetPadColor(icol);
  // gStyle->SetStatColor(icol);
  // gStyle->SetPaperSize(20,26);
  // gStyle->SetPadTopMargin(0.05);
  // gStyle->SetPadRightMargin(0.05);
  // gStyle->SetPadBottomMargin(0.16);
  // gStyle->SetPadLeftMargin(0.16);
  // gStyle->SetTitleXOffset(1.3);
  // gStyle->SetTitleYOffset(1.3);
  // gStyle->SetTextFont(font);
  // gStyle->SetTextSize(tsize);
  // gStyle->SetLabelFont(font,"x");
  // gStyle->SetTitleFont(font,"x");
  // gStyle->SetLabelFont(font,"y");
  // gStyle->SetTitleFont(font,"y");
  // gStyle->SetLabelFont(font,"z");
  // gStyle->SetTitleFont(font,"z");
  // gStyle->SetLabelSize(tsize,"x");
  // gStyle->SetTitleSize(tsize,"x");
  // gStyle->SetLabelSize(tsize,"y");
  // gStyle->SetTitleSize(tsize,"y");
  // gStyle->SetLabelSize(tsize,"z");
  // gStyle->SetTitleSize(tsize,"z");
  // gStyle->SetMarkerStyle(20);
  // gStyle->SetMarkerSize(1.2);
  // gStyle->SetHistLineWidth(2);
  // gStyle->SetLineStyleString(2,"[12 12]");
  // gStyle->SetEndErrorSize(0.);
  // gStyle->SetOptTitle(0);
  // gStyle->SetOptStat(1111111);
  // gStyle->SetOptFit(0);
  // gStyle->SetPadTickX(1);
  // gStyle->SetPadTickY(1);

  // make the post fit now
  // ----------------------------------------------------------  
  {
  TCanvas *c = new TCanvas();

  TPad *pad1 = new TPad("pad1","pad1",0,0.25,1,1);
   pad1->SetBottomMargin(0.15);
   pad1->SetGridy();
   pad1->SetGridx();
   pad1->SetTickx();
   pad1->SetTicky();
   pad1->Draw();
   pad1->cd();

  h_dat->Draw("PE");
  h_ztt->Add(h_zll);
  h_ztt->Add(h_wln);
  h_ztt->Add(h_top);
  h_ztt->Add(h_oth);
  h_ztt->Add(h_qcd);
  h_ztt->Print("all");
  h_ztt->SetFillColor(0);
  h_ztt->SetLineWidth(2);
  h_ztt->SetMarkerStyle(20);
  h_ztt->SetMarkerSize(0.2);
  h_dat->SetMarkerColor(1);
  h_dat->SetMarkerStyle(20);
  h_dat->GetYaxis()->SetTitle("Events");
  h_dat->GetXaxis()->SetTitle("MMC mass [GeV]");
  h_dat->Draw("PE");

  h_ztt->Draw("HIST,same");

  // let me compare ztt and all
  cout << "comparison starts" << endl;
  h_ztt->Print("all");
  cout << "middle" << endl;
  h_all->Print("all");
  cout << "comparison ends" << endl;
  h_wln->Add(h_zll);
  h_wln->Add(h_top);
  h_wln->Add(h_oth);
  h_wln->Add(h_qcd);
  h_wln->SetFillColor(6);
  h_wln->Draw("HIST,same");

  h_zll->Add(h_top);
  h_zll->Add(h_oth);
  h_zll->Add(h_qcd);
  h_zll->SetFillColor(9);
  h_zll->Draw("HIST,same");

  h_top->Add(h_oth);
  h_top->Add(h_qcd);
  h_top->SetFillColor(2);
  h_top->Draw("HIST,same");

  h_oth->Add(h_qcd);
  h_oth->SetFillColor(5);
  h_oth->Draw("HIST,same");

  h_qcd->SetFillColor(15);
  h_qcd->Draw("HIST,same");
  //h_ztt->Draw("LE,same");

  // h_all_up->SetLineStyle(2);
  // h_all_do->SetLineStyle(2);

  // h_all_up->Draw("hist,same");
  // h_all_do->Draw("hist,same");

  h_dat->Draw("PE,same");
  TLegend *leg1 = new TLegend(0.15+0.49, 0.55+0.15, 0.45+0.49, 0.77+0.15);
    leg1->SetFillColor(kWhite);
    leg1->SetFillStyle(0);
    leg1->SetBorderSize(0);
    leg1->AddEntry(h_ztt, "Z#rightarrow#tau#tau", "F");
    leg1->AddEntry(h_wln, "W+jets", "F");
    leg1->AddEntry(h_zll, "Z#rightarrow ll", "F");
    leg1->AddEntry(h_top, "t#bar{t}", "F");
    leg1->AddEntry(h_oth, "Diboson,single top", "F");
    leg1->AddEntry(h_qcd, "Multijet", "F");
    leg1->Draw("same");

    TLatex *l = new TLatex;
    l->SetNDC();
    l->SetTextFont(72);
    l->SetTextColor(1);
    l->SetTextSize(0.04);
    l->DrawLatex(0.15,0.96,CHANNEL+ " post-fit");


    pad1->RedrawAxis();
    c->cd();
    
    TPad *pad2 = new TPad("pad2","pad2",0,0,1,0.3);
    pad2->SetTopMargin(0.05);
    pad2->SetBottomMargin(0.15);
    pad2->SetGridy();
    pad2->SetGridx();
    pad2->SetTickx();
    pad2->SetTicky();
    pad2->Draw();
    pad2->cd();


    TH1F *h_ratio = new TH1F(*h_dat);
    h_ratio->Divide(h_ztt);
    h_ratio->GetYaxis()->SetLabelSize(0.1);
    h_ratio->GetYaxis()->SetTitleSize(0.15);
    h_ratio->GetYaxis()->SetTitleOffset(0.4);
    h_ratio->GetXaxis()->SetLabelSize(0.1);
    h_ratio->GetYaxis()->SetRangeUser(0.7, 1.3);
    if (CHANNEL.Contains("veto"))
      h_ratio->GetYaxis()->SetRangeUser(0.9, 1.1);
    h_ratio->GetYaxis()->SetTitle("Data/Prediction");
    h_ratio->GetXaxis()->SetTitle("");

    h_ratio->Draw("P");

    // other ratios
    // TH1F *h_ratio_up = new TH1F(*h_all_up);
    // h_ratio_up->Divide(h_all);
    // h_ratio_up->SetLineColor(2);
    // h_ratio_up->Draw("L,same");

    // TH1F *h_ratio_do = new TH1F(*h_all_do);
    // h_ratio_do->Divide(h_all);
    // h_ratio_do->SetLineColor(4);
    // h_ratio_do->Draw("L,same");
  

    pad2->RedrawAxis();

    c->Print(CHANNEL+"_postfit.png");
  }
    ///////// prefit plots --------------------------------------------

    {
  TCanvas *c1 = new TCanvas();

  TPad *pad1 = new TPad("pad1","pad1",0,0.25,1,1);
   pad1->SetBottomMargin(0.15);
   pad1->SetGridy();
   pad1->SetGridx();
   pad1->SetTickx();
   pad1->SetTicky();
   pad1->Draw();
   pad1->cd();

  h_dat->Draw("PE");
  h_ztt_pre->Add(h_zll_pre);
  h_ztt_pre->Add(h_wln_pre);
  h_ztt_pre->Add(h_top_pre);
  h_ztt_pre->Add(h_oth_pre);
  h_ztt_pre->Add(h_qcd_pre);
  h_ztt_pre->Print("all");
  h_ztt_pre->SetFillColor(0);
  h_ztt_pre->SetLineWidth(2);
  h_ztt_pre->SetMarkerStyle(20);
  h_ztt_pre->SetMarkerSize(0.2);
  h_dat->SetMarkerColor(1);
  h_dat->SetMarkerStyle(20);
  h_dat->GetYaxis()->SetTitle("Events");
  h_dat->GetXaxis()->SetTitle("MMC mass [GeV]");
  h_dat->Draw("PE");
       
  h_ztt_pre->Draw("HIST,same");
  h_ztt_pre->Draw("HIST,same");

  h_wln_pre->Add(h_zll_pre);
  h_wln_pre->Add(h_top_pre);
  h_wln_pre->Add(h_oth_pre);
  h_wln_pre->Add(h_qcd_pre);
  h_wln_pre->SetFillColor(6);
  h_wln_pre->Draw("HIST,same");

  h_zll_pre->Add(h_top_pre);
  h_zll_pre->Add(h_oth_pre);
  h_zll_pre->Add(h_qcd_pre);
  h_zll_pre->SetFillColor(9);
  h_zll_pre->Draw("HIST,same");

  h_top_pre->Add(h_oth_pre);
  h_top_pre->Add(h_qcd_pre);
  h_top_pre->SetFillColor(2);
  h_top_pre->Draw("HIST,same");
       
  h_oth_pre->Add(h_qcd_pre);
  h_oth_pre->SetFillColor(5);
  h_oth_pre->Draw("HIST,same");

  h_qcd_pre->SetFillColor(15);
  h_qcd_pre->Draw("HIST,same");
  //  h_ztt_pre->Draw("LE,same");

  h_dat->Draw("PE,same");
  TLegend *leg1 = new TLegend(0.15+0.49, 0.55+0.15, 0.45+0.49, 0.77+0.15);
    leg1->SetFillColor(kWhite);
    leg1->SetFillStyle(0);
    leg1->SetBorderSize(0);
    leg1->AddEntry(h_ztt, "Z#rightarrow#tau#tau", "F");
    leg1->AddEntry(h_wln, "W+jets", "F");
    leg1->AddEntry(h_zll, "Z#rightarrow ll", "F");
    leg1->AddEntry(h_top, "t#bar{t}", "F");
    leg1->AddEntry(h_oth, "Diboson,single top", "F");
    leg1->AddEntry(h_qcd, "Multijet", "F");
    leg1->Draw("same");

    TLatex *l = new TLatex;
    l->SetNDC();
    l->SetTextFont(72);
    l->SetTextColor(1);
    l->SetTextSize(0.04);
    l->DrawLatex(0.15,0.96,CHANNEL+ " pre-fit");


    pad1->RedrawAxis();
    c1->cd();
    
    TPad *pad2 = new TPad("pad2","pad2",0,0,1,0.3);
    pad2->SetTopMargin(0.05);
    pad2->SetBottomMargin(0.15);
    pad2->SetGridy();
    pad2->SetGridx();
    pad2->SetTickx();
    pad2->SetTicky();
    pad2->Draw();
    pad2->cd();


    TH1F *h_ratio = new TH1F(*h_dat);
    h_ratio->Divide(h_ztt_pre);
    h_ratio->GetYaxis()->SetLabelSize(0.1);
    h_ratio->GetYaxis()->SetTitleSize(0.15);
    h_ratio->GetYaxis()->SetTitleOffset(0.4);
    h_ratio->GetXaxis()->SetLabelSize(0.1);
    h_ratio->GetYaxis()->SetRangeUser(0.7, 1.3);
    h_ratio->GetYaxis()->SetTitle("Data/Prediction");
    h_ratio->GetXaxis()->SetTitle("");

    h_ratio->Draw("P");

   
    pad2->RedrawAxis();

    c1->Print(CHANNEL+"_prefit.png");
    }


}