#include <iomanip>
#include <fstream>
#include <iostream>
#include <sstream>
#include "TH2.h"
#include "TF2.h"
#include "TLine.h"
#include "TFile.h"
#include "TStyle.h"
#include "TLatex.h"
#include "TRandom.h"
#include "TString.h"
#include "TLegend.h"
#include "TRandom.h"
#include "TCanvas.h"
#include "THStack.h"
#include "TPostScript.h"
#include "TGraphErrors.h"

Double_t getGauss( TRandom &gRand, const Double_t mean, const Double_t error );
Double_t getExponential(TRandom &gRand, const Int_t nXsec, const Double_t step, const Double_t tau );
void higgs_combined_limit_cdfd0(Int_t NoSYS=0); // limits with various treatment of systematicNoS: 1=no eff, 2=no back, 3=no sys,4=no correlation between signal and backgrounds; NoSYS=10 for shape study for zh->llbb; NoSYS=11, randon pick two shapes; NoSYS=20; uncorrelated; NoSYS>100 for individual masses       
void higgs_combined_limit_pseudoexperiment_cdfd0(Int_t NoSYS=0 ); // doing pseudo-experiments 

//Int_t Nbins = 0; 
static const bool debug = false; // for debug purpose 
static const Int_t NCHA = 7; // number of channels to combine (lnubb_s,lnubb_d,nunubb,ww)  
static const Int_t NCHAD0=14; // number of channels to combine from D0
static const Int_t NCHAT =NCHA+NCHAD0;  
static const Int_t NCOMM =5; // number of common systematic between CDF and D0 
static const Int_t NACC = 6; // source of systematic for acceptance(lum/btag/lepton/JES/I(F)SR+PDF/Other)
static const Int_t NACCD0 = 19; // source of systematic from D0 
static const Int_t NBKG = 10;  // source of backgrounds systematic from CDF 
static const Int_t NBKGH = 6; // source of backgrounds(HF/mistag/top/nonW/diboson/other) 
static const Int_t NBKGD0 = 23; // source of background systematic from D0 
static const Int_t NBKGHD0 = 8; // bkgd source of histograms for d0  
static const Int_t NMASS = 10; // Higgs mass 110-180 with 10 GeV steps, we can extrapolate in 5 GeV  
static const Int_t HMass[NMASS]={110,115,120,130,140,150,160,170,180,200}; // higgs mass for this 10 points
static Int_t icha[NCHAT]={0}; // the channel to include in combination 
static float scomm[NCOMM]={0.}; // common systematic (sig_lumi, sig_Xsec, bkgd_lumi, bkgd_Xsec_top, bkgd_Xsec_EW)
static float seff[NACC][NCHA]={0.}; // values of acceptance systematic
static float seff_d0[NACCD0][NCHAD0]={0.};   
static float sbkg[NBKG][NMASS][NCHA]={0.}; // values of background systematic 
static float sbkg_d0[NBKGD0][NMASS][NCHAD0]={0.};
static float sjes_nnbb[2][NBKGH]={0.}; // special jes for qcd, top, w+hf,z+hf, diboson for vvbb from CDF   
static const Int_t Nmax[NCHAT] ={6, 6, 4, 10,4,4,6,6,6,6,6,5,5,10,10,10,10,10,10,6,6}; 
static const float sfnnlo[NCHAT]={0.96,0.96,1.29,1.45,0.99,0.99,0.99,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.}; // NNLO Higgs Xsec/NLO:gg->H has 10% additional uncertainties (1.3 is nunubb mass window eff)
static const Int_t Npseudo=1000; // doing psudo-experiment 
TH1F* XDATA[NMASS][NCHAT] = {0}; // the point to data histograms (make sure some of it are 2-d)  
TH1F* XBKG[NBKGH][NMASS][NCHA] = {0}; // the point to background histograms
TH1F* XBKG_d0[NBKGHD0][NMASS][NCHAD0] = {0};  
TH1F* XHIGGS[NMASS][NCHAT] ={0}; // higgs signal 
TH1F* XMC[NCHAT] ={0}; // histograms for pseudo-experiments 
static Double_t lgm[200]={0.}; 
//_____________________________________________________________________________
Double_t getGauss( TRandom &gRand, const Double_t mean, const Double_t error )
{
  Double_t value = -1.;
  do {
    value = gRand.Gaus(mean,error);
  } while ( value <= 0. );

  return value;
}

//_____________________________________________________________________________
Double_t getExponential( TRandom &gRand, const Int_t nXsec, const Double_t step, const Double_t tau )
{
  Double_t value = -1.;
  do {
    value = gRand.Exp(tau);
  } while ( value/step >= nXsec );

  return value;
}

//___________________________________________________________________________________
// Computing the combined limit with the channels defined in higgs_combined_limit.inp 
//-----------------------------------------------------------------------------------
void higgs_combined_limit_cdfd0(Int_t NoSYS)
{
  // Input files: higgs_combined_limit.inp and many various root files needed 
  // Output the plot and limits 
  gStyle->SetOptStat(0);
  gStyle->SetOptTitle(0);
  //
  TString txt[NCHAT] ={"lnubb_s","lnubb_d","nunubb","ww","nunubb_s", "nunubb_d", "llbb","d0evbb_s","d0evbb_d","d0mvbb_s","d0mvbb_d","d0vvbb_s","d0vvbb_d","d0wwee","d0wwem","d0wwmm", "d0whee", "d0whem","d0whmm","d0zhee","d0zhmm"};

  Double_t nfact(0); 
  lgm[0] = 0.; 
  for(Int_t i =1; i<200; ++i){ 
    nfact +=TMath::Log(i); 
    lgm[i] = nfact; 
  }
  char junk[120] =" ";
  FILE * infile = fopen("higgs_combined_limit.inp","r"); 
  for(int i=0; i<NCHA; ++i) {
    fscanf(infile, "%i", &icha[i]);
    if(debug)std::cout<<" icha[ "<< i<<" ]= "<<icha[i]<<std::endl; 
  }

  fgets(junk,60, infile);
  for(int j=0; j<NACC; ++j){
  for(int i=0; i<NCHA; ++i){
    fscanf(infile, "%f", &seff[j][i]); 
    if(debug)std::cout<<" seff[ "<< j<<" ]["<<i<<"] ="<<seff[j][i]<<std::endl; 
    if(NoSYS==20&&j>1)seff[j][i] = -1.*fabs(seff[j][i]); // treat the uncertainties uncorrected except lum, btag 
  }
  fgets(junk,60, infile);
  }

  if(NoSYS==11) seff[NACC-1][NCHA-1] = 0.; // reset other systematic in zh->llbb channel when using different shape 

  for(int j=0; j<NBKG-1; ++j){
  for(int i=0; i<NCHA; ++i) {
    fscanf(infile, "%f", &sbkg[j][0][i]); 
    if(debug)std::cout<< " sbkg["<<j<<"][0]["<<i<<"]="<<sbkg[j][0][i]<<std::endl; 
    if(NoSYS==20)sbkg[j][0][i] = -1.*fabs(sbkg[j][0][i]); 
  }
  fgets(junk,60, infile);
  }
  
  for(int i=0; i<NMASS; ++i) {
    fscanf(infile, "%f", &sbkg[NBKG-1][i][3]); // ww is id 3 
    if(debug)std::cout<<" sbkg["<<NBKG-1<<"]["<<i<<"]["<<3<<"]="<<sbkg[NBKG-1][i][3]<<std::endl; 
  }
  fgets(junk,60, infile); 
  for(int i=0; i<NBKGH; ++i){ 
    fscanf(infile, "%f", &sjes_nnbb[0][i]); 
    if(debug)std::cout<<" sjes_nnbb[0]["<<i<<"]="<<sjes_nnbb[0][i]<<std::endl; 
    if(NoSYS==20) sjes_nnbb[0][i] = -1.*fabs(sjes_nnbb[0][i]); 
  }
  fgets(junk,60, infile); 
  for(int i=0; i<NBKGH; ++i){ 
    fscanf(infile, "%f", &sjes_nnbb[1][i]); 
    if(debug)std::cout<<" sjes_nnbb[1]["<<i<<"]="<<sjes_nnbb[1][i]<<std::endl; 
    if(NoSYS==20) sjes_nnbb[1][i] = -1.*fabs(sjes_nnbb[1][i]); 
  }

  for(int j=0; j<NBKG-1; ++j){
    for(int i=0; i<NCHA; ++i){ 
      for(int k=1; k<NMASS; ++k){ 
	sbkg[j][k][i] = sbkg[j][0][i]; 
      }
    }
  }
  // read in d0 inputs 
  infile = fopen("dumpD0_data.inp","r"); 
  for(int i=0; i<NCHAD0; ++i) {
    fscanf(infile, "%i", &icha[NCHA+i]);
    if(debug)std::cout<<" icha[ "<< NCHA+i<<" ]= "<<icha[NCHA+i]<<std::endl; 
  }
  fgets(junk,60, infile);
  for(int j=0; j<NCOMM; ++j){
    fscanf(infile, "%f", &scomm[j]); 
    if(debug)std::cout<<" scomm[ "<< j<<"] ="<<scomm[j]<<std::endl; 
  fgets(junk,110, infile);
  }

  for(int j=0; j<NACCD0; ++j){
  for(int i=0; i<NCHAD0; ++i){
    fscanf(infile, "%f", &seff_d0[j][i]); 
    if(debug)std::cout<<" seff_d0[ "<< j<<" ]["<<i<<"] ="<<seff_d0[j][i]<<std::endl; 
  }
  fgets(junk,60, infile);
  }

  for(int j=0; j<NBKGD0; ++j){
  for(int i=0; i<NCHAD0; ++i) {
    fscanf(infile, "%f", &sbkg_d0[j][0][i]); 
    if(debug)std::cout<< " sbkg_d0["<<j<<"][0]["<<i<<"]="<<sbkg_d0[j][0][i]<<std::endl; 
  }
  fgets(junk,60, infile);
  }
  
  for(int j=0; j<NBKGD0; ++j){
    for(int i=0; i<NCHAD0; ++i){ 
      for(int k=1; k<NMASS; ++k){ 
	sbkg_d0[j][k][i] = sbkg_d0[j][0][i]; 
      }
    }
  }

  // correct the tevatron common systematic for cdf 
  for(int i =0; i<NCHA; ++i){ 
    seff[0][i] = sqrt(pow(seff[0][i],2)-pow(scomm[0],2)); 
  }

  for(int i=0; i<NCHA; ++i){ // top xsec  
    for(int k=0; k<NMASS; ++k){
      if(fabs(sbkg[3][k][i])>0.){
	sbkg[3][k][i] = sbkg[3][k][i]/fabs(sbkg[3][k][i])*sqrt(pow(sbkg[3][k][i],2)-pow(scomm[3],2)-pow(scomm[0],2));
      }
    }
  }

  for(int i=0; i<NCHA; ++i){ // EW xsec 
    for(int k=0; k<NMASS; ++k){
      if(fabs(sbkg[4][k][i])>0.){
	sbkg[4][k][i] = sbkg[4][k][i]/fabs(sbkg[4][k][i])*sqrt(pow(sbkg[4][k][i],2)-pow(scomm[4],2)-pow(scomm[0],2));
      }
    }
  }

  TFile* f[2];
  f[0] = new TFile("shape_snntag_1fb.root","read"); 
  f[1] = new TFile("shape_secvtx_1fb.root","read");
  TFile* f3 = new TFile("shape_nunubb.root","read");   
  TFile* f4 = new TFile("shape_hww.root","read"); 
  TFile* f5 = new TFile("cdfvvbb973.root","read"); 
  TFile* f6 = new TFile("ZHllbb_1fb.root","read"); 

  TFile* f1_d0[6]; 
  f1_d0[0] = new TFile("fv_whevbb_1tag.root","read"); 
  f1_d0[1] = new TFile("fv_whevbb_2tag.root","read"); 
  f1_d0[2] = new TFile("fv_whmvbb_1tag.root","read"); 
  f1_d0[3] = new TFile("fv_whmvbb_2tag.root","read"); 
  f1_d0[4] = new TFile("fv_zhvvbb_1Tag.root","read"); 
  f1_d0[5] = new TFile("fv_zhvvbb_2Tag.root","read"); 

  TFile* f2_d0[3];   
  f2_d0[0] = new TFile("fv_hww_ee.root","read"); 
  f2_d0[1] = new TFile("fv_hww_emu.root","read"); 
  f2_d0[2] = new TFile("fv_hww_mumu_1fb.root","read"); 

  TFile* f3_d0[3]; 
  f3_d0[0] = new TFile("fv_www_ee.root","read"); 
  f3_d0[1] = new TFile("fv_www_emu.root","read"); 
  f3_d0[2] = new TFile("fv_www_mumu.root","read"); 

  TFile* f4_d0[2]; 
  f4_d0[0] = new TFile("fv_zheebb_2tag.root","read"); 
  f4_d0[1] = new TFile("fv_zhmmbb_2tag.root","read"); 


  for(int i=0; i<2; ++i){ 
    XDATA[0][i]  = (TH1F*)f[i]->Get("M2j_data"); 
    XDATA[0][i]->Rebin(4); 
    XHIGGS[0][i] = (TH1F*)f[i]->Get("M2j_WH110_CTEQ5L_ISRon_FSRon"); 
    XHIGGS[1][i] = (TH1F*)f[i]->Get("M2j_WH115_CTEQ5L_ISRon_FSRon"); 
    XHIGGS[2][i] = (TH1F*)f[i]->Get("M2j_WH120_CTEQ5L_ISRon_FSRon"); 
    XHIGGS[3][i] = (TH1F*)f[i]->Get("M2j_WH130_CTEQ5L_ISRon_FSRon"); 
    XHIGGS[4][i] = (TH1F*)f[i]->Get("M2j_WH140_CTEQ5L_ISRon_FSRon"); 
    XHIGGS[5][i] = (TH1F*)f[i]->Get("M2j_WH150_CTEQ5L_ISRon_FSRon"); 
    XBKG[0][0][i] = (TH1F*)f[i]->Get("M2j_mistag"); 
    XBKG[1][0][i] = (TH1F*)f[i]->Get("M2j_nonW"); 
    XBKG[2][0][i] = (TH1F*)f[i]->Get("M2j_HF");     
    XBKG[3][0][i] = (TH1F*)f[i]->Get("M2j_Top"); 
    XBKG[4][0][i] = (TH1F*)f[i]->Get("M2j_diboson");
    for(int j = 0; j<6; ++j){ 
      XHIGGS[j][i] ->Rebin(4); 
      if(j<5)XBKG[j][0][i] ->Rebin(4);       
    }


    for(int j =1; j<NMASS; ++j){
      XDATA[j][i] = XDATA[0][i]; 
      for(int k=0; k<NBKGH; ++k){ 
	XBKG[k][j][i] = XBKG[k][0][i];
      }
    }
  }
  // what about ZH->nunubb ?   
    XDATA[0][2]  = (TH1F*)f3->Get("Dijet_data"); 
    XHIGGS[0][2] = (TH1F*)f3->Get("Dijet_mh110"); 
    XHIGGS[2][2] = (TH1F*)f3->Get("Dijet_mh120");
    XHIGGS[1][2] = (TH1F*)XHIGGS[0][2]->Clone(); 
    XHIGGS[1][2] ->Add(XHIGGS[2][2]); 
    XHIGGS[1][2] ->Scale(0.5); 
    XHIGGS[3][2] = (TH1F*)f3->Get("Dijet_mh130"); 
    
    std::cout<<" CDF ZH->nunubb mass "<< HMass[0]<< "signal "<<XHIGGS[0][2]->Integral()<<std::endl; 
    std::cout<<" CDF ZH->nunubb mass "<< HMass[1]<< "signal "<<XHIGGS[1][2]->Integral()<<std::endl; 
    std::cout<<" CDF ZH->nunubb mass "<< HMass[2]<< "signal "<<XHIGGS[2][2]->Integral()<<std::endl; 
    std::cout<<" CDF ZH->nunubb mass "<< HMass[3]<< "signal "<<XHIGGS[3][2]->Integral()<<std::endl; 

    XBKG[0][0][2] = (TH1F*)f3->Get("Dijet_mistag"); 
    XBKG[1][0][2] = (TH1F*)f3->Get("Dijet_qcd"); 
    XBKG[2][0][2] = (TH1F*)f3->Get("Dijet_hf");     
    XBKG[3][0][2] = (TH1F*)f3->Get("Dijet_top"); 
    XBKG[4][0][2] = (TH1F*)f3->Get("Dijet_diboson");

    for(int j =1; j<NMASS; ++j){
      XDATA[j][2] = XDATA[0][2]; 
      for(int k=0; k<NBKGH; ++k){ 
	XBKG[k][j][2] = XBKG[k][0][2];
      }
    }
  // what about ww->llnunu 
  for(int i = 0; i<10; ++i){ 
    TString a ="data_";
    TString a1 ="hww_";
    TString a2 ="ww_";
    TString a3 ="bg_";
    int i0 = HMass[i]; 
    a +=i0; 
    a1 +=i0; 
    a2 +=i0; 
    a3 +=i0;
    XDATA[i][3] = (TH1F*) f4->Get(a); 
    XHIGGS[i][3] = (TH1F*) f4->Get(a1); 
    XBKG[4][i][3] = (TH1F*) f4->Get(a2);
    XBKG[5][i][3] = (TH1F*) f4->Get(a3);
    std::cout<<" Hww Hmass="<<i0<<" data="<<XDATA[i][3]->Integral()<<" XHIGGS="<< XHIGGS[i][3]->Integral()<<" WW="<<XBKG[4][i][3]->Integral()<<" Other="<<XBKG[5][i][3]->Integral()<<std::endl; 
  }
  // new zh->nunubb single tag analysis 
    XDATA[0][4]  = (TH1F*)f5->Get("mjj_DATA_1Tag_CENT");
    XDATA[0][4] ->Rebin(2); 
    XHIGGS[0][4] = (TH1F*)f5->Get("mjj_zh110_1Tag_CENT"); 
    XHIGGS[0][4] ->Add((TH1F*)f5->Get("mjj_wh110_1Tag_CENT"));
    XHIGGS[1][4] = (TH1F*)f5->Get("mjj_zh115_1Tag_CENT"); 
    XHIGGS[1][4] ->Add((TH1F*)f5->Get("mjj_wh115_1Tag_CENT"));
    XHIGGS[2][4] = (TH1F*)f5->Get("mjj_zh120_1Tag_CENT"); 
    XHIGGS[2][4] ->Add((TH1F*)f5->Get("mjj_wh120_1Tag_CENT"));
    XHIGGS[3][4] = (TH1F*)f5->Get("mjj_zh130_1Tag_CENT"); 
    XHIGGS[3][4] ->Add((TH1F*)f5->Get("mjj_wh130_1Tag_CENT"));
    for(int i=0; i<4; i++){ 
      XHIGGS[i][4]->Rebin(2); 
    }
    XBKG[0][0][4] = (TH1F*)f5->Get("mjj_MTAG_1Tag_CENT"); // mistag 
    XBKG[1][0][4] = (TH1F*)f5->Get("mjj_QCD1_1Tag_CENT"); // QCD
    XBKG[1][0][4] ->Add((TH1F*)f5->Get("mjj_QCD2_1Tag_CENT")); 
    XBKG[2][0][4] = (TH1F*)f5->Get("mjj_Wev_1Tag_CENT"); // W+HF 
    XBKG[2][0][4] ->Add((TH1F*)f5->Get("mjj_Wmv_1Tag_CENT")); 
    XBKG[2][0][4] ->Add((TH1F*)f5->Get("mjj_Wtauv_1Tag_CENT")); 
    XBKG[3][0][4] = (TH1F*)f5->Get("mjj_ttbar_1Tag_CENT"); //top  
    XBKG[3][0][4] ->Add((TH1F*)f5->Get("mjj_tops_1Tag_CENT")); 
    XBKG[3][0][4] ->Add((TH1F*)f5->Get("mjj_topt_1Tag_CENT")); 
    XBKG[4][0][4] = (TH1F*)f5->Get("mjj_ZZ_1Tag_CENT"); //diboson  
    XBKG[4][0][4] ->Add((TH1F*)f5->Get("mjj_WZ_1Tag_CENT")); 
    XBKG[4][0][4] ->Add((TH1F*)f5->Get("mjj_WW_1Tag_CENT")); 
    XBKG[5][0][4] = (TH1F*)f5->Get("mjj_Zmm_1Tag_CENT"); //Z+HF  
    XBKG[5][0][4] ->Add((TH1F*)f5->Get("mjj_Ztautau_1Tag_CENT")); 
    XBKG[5][0][4] ->Add((TH1F*)f5->Get("mjj_Zvv_1Tag_CENT")); 
    for(int i=0; i<6; i++){ 
      XBKG[i][0][4]->Rebin(2); 
    }
    
    for(int j =1; j<NMASS; ++j){
      XDATA[j][4] = XDATA[0][4]; 
      for(int k=0; k<NBKGH; ++k){ 
	XBKG[k][j][4] = XBKG[k][0][4];
      }
    }

  // new zh->nunubb double tag analysis 
    XDATA[0][5]  = (TH1F*)f5->Get("mjj_DATA_2Tag_CENT"); 
    XDATA[0][5] ->Rebin(2); 
    XHIGGS[0][5] = (TH1F*)f5->Get("mjj_zh110_2Tag_CENT"); 
    XHIGGS[0][5] ->Add((TH1F*)f5->Get("mjj_wh110_2Tag_CENT"));
    XHIGGS[1][5] = (TH1F*)f5->Get("mjj_zh115_2Tag_CENT"); 
    XHIGGS[1][5] ->Add((TH1F*)f5->Get("mjj_wh115_2Tag_CENT"));
    XHIGGS[2][5] = (TH1F*)f5->Get("mjj_zh120_2Tag_CENT"); 
    XHIGGS[2][5] ->Add((TH1F*)f5->Get("mjj_wh120_2Tag_CENT"));
    XHIGGS[3][5] = (TH1F*)f5->Get("mjj_zh130_2Tag_CENT"); 
    XHIGGS[3][5] ->Add((TH1F*)f5->Get("mjj_wh130_2Tag_CENT"));
    for(int i =0; i<4; ++i){ 
      XHIGGS[i][5] ->Rebin(2); 
    }
    XBKG[0][0][5] = (TH1F*)f5->Get("mjj_MTAG_2Tag_CENT"); // mistag 
    XBKG[1][0][5] = (TH1F*)f5->Get("mjj_QCD1_2Tag_CENT"); // QCD
    XBKG[1][0][5] ->Add((TH1F*)f5->Get("mjj_QCD2_2Tag_CENT")); 
    XBKG[2][0][5] = (TH1F*)f5->Get("mjj_Wev_2Tag_CENT"); // W+HF 
    XBKG[2][0][5] ->Add((TH1F*)f5->Get("mjj_Wmv_2Tag_CENT")); 
    XBKG[2][0][5] ->Add((TH1F*)f5->Get("mjj_Wtauv_2Tag_CENT")); 
    XBKG[3][0][5] = (TH1F*)f5->Get("mjj_ttbar_2Tag_CENT"); //top  
    XBKG[3][0][5] ->Add((TH1F*)f5->Get("mjj_tops_2Tag_CENT")); 
    XBKG[3][0][5] ->Add((TH1F*)f5->Get("mjj_topt_2Tag_CENT")); 
    XBKG[4][0][5] = (TH1F*)f5->Get("mjj_ZZ_2Tag_CENT"); //diboson  
    XBKG[4][0][5] ->Add((TH1F*)f5->Get("mjj_WZ_2Tag_CENT")); 
    XBKG[4][0][5] ->Add((TH1F*)f5->Get("mjj_WW_2Tag_CENT")); 
    XBKG[5][0][5] = (TH1F*)f5->Get("mjj_Zmm_2Tag_CENT"); //Z+HF  
    XBKG[5][0][5] ->Add((TH1F*)f5->Get("mjj_Ztautau_2Tag_CENT")); 
    XBKG[5][0][5] ->Add((TH1F*)f5->Get("mjj_Zvv_2Tag_CENT")); 
    for(int i =0; i<6; ++i){ 
      XBKG[i][0][5] ->Rebin(2); 
    }
    for(int j =1; j<NMASS; ++j){
      XDATA[j][5] = XDATA[0][5]; 
      for(int k=0; k<NBKGH; ++k){ 
	XBKG[k][j][5] = XBKG[k][0][5];
      }
    }
    // get zh->llbb  (these are 2-d, but treat as 1-d for now) 
    XDATA[0][6]  = (TH1F*)f6->Get("DATATot");
    XHIGGS[0][6] = (TH1F*)f6->Get("ZH110Tot"); //zh110 
    XHIGGS[1][6] = (TH1F*)f6->Get("ZH115Tot"); //zh115 
    XHIGGS[2][6] = (TH1F*)f6->Get("ZH120Tot"); //zh120
    XHIGGS[3][6] = (TH1F*)f6->Get("ZH130Tot"); //zh130
    XHIGGS[4][6] = (TH1F*)f6->Get("ZH140Tot"); //zh140
    XHIGGS[5][6] = (TH1F*)f6->Get("ZH150Tot"); //zh150

    XBKG[0][0][6] = (TH1F*)f6->Get("ZmisTot"); // mistag 
    XBKG[1][0][6] = (TH1F*)f6->Get("FakeTot"); // QCD or fake 
    XBKG[2][0][6] = (TH1F*)f6->Get("ZbbTot"); // Z+bb 
    if(NoSYS==10) XBKG[2][0][6] = (TH1F*)f6->Get("ZbbshapeTot"); // Z+bb 
    XBKG[3][0][6] = (TH1F*)f6->Get("TTTot"); //top  
    XBKG[4][0][6] = (TH1F*)f6->Get("ZZTot"); //diboson  
    XBKG[4][0][6] ->Add((TH1F*)f6->Get("ZWTot")); 
    XBKG[5][0][6] = (TH1F*)f6->Get("ZccTot"); //Z+cc  
    if(NoSYS==10) XBKG[5][0][6] = (TH1F*)f6->Get("ZccshapeTot"); //Z+cc  
    for(int j =1; j<NMASS; ++j){
      XDATA[j][6] = XDATA[0][6]; 
      for(int k=0; k<NBKGH; ++k){ 
	XBKG[k][j][6] = XBKG[k][0][6];
      }
    }
   
    TH2F* XZHshape[2]; // shape systematic for zh->llbb pseudo experiment for generating different shape  
    XZHshape[0] = (TH2F*)f6->Get("ZbbshapeTot"); // Z+bb
    XZHshape[1] =  (TH2F*)f6->Get("ZccshapeTot"); //Z+cc  


  // get histograms for d0 
  for(int i=0; i<2; ++i){ //whevbb
    
    XDATA[0][i+NCHA]  = (TH1F*)f1_d0[i]->Get("Data Final Variable (mH)");
    XDATA[0][i+NCHA] ->Rebin(5); 
    XHIGGS[0][i+NCHA] = (TH1F*)f1_d0[i]->Get("Signal Interpolate 105 + 5");
    XHIGGS[1][i+NCHA] = (TH1F*)f1_d0[i]->Get("Signal Interpolate 115 + 0");  
    XHIGGS[2][i+NCHA] = (TH1F*)f1_d0[i]->Get("Signal Interpolate 115 + 5"); 
    XHIGGS[3][i+NCHA] = (TH1F*)f1_d0[i]->Get("Signal Interpolate 125 + 5");
    XHIGGS[4][i+NCHA] = (TH1F*)f1_d0[i]->Get("Signal Interpolate 135 + 5");
    XHIGGS[5][i+NCHA] = (TH1F*)f1_d0[i]->Get("Signal Interpolate 145 + 5");
    for(int j =0; j<5; ++j)XHIGGS[j][i+NCHA]->Rebin(5); 
    XBKG_d0[0][0][i] = (TH1F*)f1_d0[i]->Get("WZ-Bkgd Final Variable (mH)"); 
    XBKG_d0[1][0][i] = (TH1F*)f1_d0[i]->Get("Wbb-Bkgd Final Variable (mH)"); 
    XBKG_d0[2][0][i] = (TH1F*)f1_d0[i]->Get("Wjj-Bkgd Final Variable (mH)"); 
    XBKG_d0[3][0][i] = (TH1F*)f1_d0[i]->Get("Zjj-Bkgd Final Variable (mH)"); 
    XBKG_d0[4][0][i] = (TH1F*)f1_d0[i]->Get("tt-Bkgd Final Variable (mH)"); 
    XBKG_d0[5][0][i] = (TH1F*)f1_d0[i]->Get("single top-Bkgd Final Variable (mH)"); 
    XBKG_d0[6][0][i] = (TH1F*)f1_d0[i]->Get("qcd-Bkgd Final Variable (mH)"); 
    for(int j=0; j<7; j++)XBKG_d0[j][0][i]->Rebin(5);

    for(int j =1; j<Nmax[i+NCHA]; ++j){
      XDATA[j][i+NCHA] = XDATA[0][i+NCHA]; 
      for(int k=0; k<NBKGHD0; ++k){ 
	XBKG_d0[k][j][i] = XBKG_d0[k][0][i];
      }
    }
  }

  for(int i=2; i<4; ++i){ // wmvbb
    
    XDATA[0][i+NCHA]  = (TH1F*)f1_d0[i]->Get("Data Final Variable (mH)");
    XDATA[0][i+NCHA] ->Rebin(5); 
    XHIGGS[0][i+NCHA] = (TH1F*)f1_d0[i]->Get("hwlv Interpolate 105 + 5"); 
    XHIGGS[1][i+NCHA] = (TH1F*)f1_d0[i]->Get("hwlv Interpolate 115 + 0"); 
    XHIGGS[2][i+NCHA] = (TH1F*)f1_d0[i]->Get("hwlv Interpolate 115 + 5"); 
    XHIGGS[3][i+NCHA] = (TH1F*)f1_d0[i]->Get("hwlv Interpolate 125 + 5");
    XHIGGS[4][i+NCHA] = (TH1F*)f1_d0[i]->Get("hwlv Interpolate 135 + 5");
    XHIGGS[5][i+NCHA] = (TH1F*)f1_d0[i]->Get("hwlv Interpolate 145 + 5");
    for(int j =0; j<5; ++j)XHIGGS[j][i+NCHA]->Rebin(5); 
    XBKG_d0[0][0][i] = (TH1F*)f1_d0[i]->Get("WZ-Bkgd Final Variable (mH)"); 
    XBKG_d0[1][0][i] = (TH1F*)f1_d0[i]->Get("Wbb-Bkgd Final Variable (mH)"); 
    XBKG_d0[2][0][i] = (TH1F*)f1_d0[i]->Get("Wjj-Bkgd Final Variable (mH)"); 
    XBKG_d0[3][0][i] = (TH1F*)f1_d0[i]->Get("Zmm-Bkgd Final Variable (mH)");
    XBKG_d0[3][0][i] ->Add((TH1F*)f1_d0[i]->Get("Wtv-Bkgd Final Variable (mH)"));
    XBKG_d0[4][0][i] = (TH1F*)f1_d0[i]->Get("tt-Bkgd Final Variable (mH)"); 
    XBKG_d0[5][0][i] = (TH1F*)f1_d0[i]->Get("st-Bkgd Final Variable (mH)"); 
    XBKG_d0[6][0][i] = (TH1F*)f1_d0[i]->Get("qcd-Bkgd Final Variable (mH)"); 
    for(int j=0; j<7; j++)XBKG_d0[j][0][i]->Rebin(5);

    for(int j =1; j<Nmax[i+NCHA]; ++j){
      XDATA[j][i+NCHA] = XDATA[0][i+NCHA]; 
      for(int k=0; k<NBKGHD0; ++k){ 
	XBKG_d0[k][j][i] = XBKG_d0[k][0][i];
      }
    }
  }


  for(int i=4; i<6; ++i){ // zhvvbb
    XDATA[0][i+NCHA]  = (TH1F*)f1_d0[i]->Get("Data Final Variable (mH)");
    XDATA[0][i+NCHA] ->Rebin(5); 
    XHIGGS[0][i+NCHA] = (TH1F*)f1_d0[i]->Get("Signal Interpolate 105 + 5"); 
    XHIGGS[1][i+NCHA] = (TH1F*)f1_d0[i]->Get("Signal Interpolate 115 + 0"); 
    XHIGGS[2][i+NCHA] = (TH1F*)f1_d0[i]->Get("Signal Interpolate 115 + 5");
    XHIGGS[3][i+NCHA] = (TH1F*)f1_d0[i]->Get("Signal Interpolate 125 + 5");
    XHIGGS[4][i+NCHA] = (TH1F*)f1_d0[i]->Get("Signal Interpolate 135 + 5");
    for(int j =0; j<4; ++j)XHIGGS[j][i+NCHA]->Rebin(5); 

    XBKG_d0[0][0][i] = (TH1F*)f1_d0[i]->Get("WZ-Bkgd Final Variable (mH)"); 
    XBKG_d0[0][0][i] ->Add((TH1F*)f1_d0[i]->Get("ZZ-Bkgd Final Variable (mH)"));
    XBKG_d0[1][0][i] = (TH1F*)f1_d0[i]->Get("Wbb-Bkgd Final Variable (mH)");
    XBKG_d0[1][0][i] ->Add((TH1F*)f1_d0[i]->Get("Zbb-Bkgd Final Variable (mH)"));
    XBKG_d0[2][0][i] = (TH1F*)f1_d0[i]->Get("Wjj-Bkgd Final Variable (mH)"); 
    XBKG_d0[3][0][i] = (TH1F*)f1_d0[i]->Get("Zjj-Bkgd Final Variable (mH)"); 
    XBKG_d0[4][0][i] = (TH1F*)f1_d0[i]->Get("tt-Bkgd Final Variable (mH)"); 
    XBKG_d0[5][0][i] = (TH1F*)f1_d0[i]->Get("qtb-Bkgd Final Variable (mH)"); 
    XBKG_d0[6][0][i] = (TH1F*)f1_d0[i]->Get("qcd-Bkgd Final Variable (mH)"); 
    for(int j=0; j<7; j++)XBKG_d0[j][0][i]->Rebin(5);

    for(int j =1; j<Nmax[i+NCHA]; ++j){
      XDATA[j][i+NCHA] = XDATA[0][i+NCHA]; 
      for(int k=0; k<NBKGHD0; ++k){ 
	XBKG_d0[k][j][i] = XBKG_d0[k][0][i];
      }
    }
  }

  // what about ww->llnunu from d0
  TString tit[NMASS]={" ", " ","120 + 0","120 + 10","140 + 0", "140 + 10","160 + 0","160 + 10","180 + 0","180 + 20"}; 
  TString tits[NMASS]={" ", " ","120","130","140","150","160","170","180","200"}; 
  for(int j=2; j<NMASS; ++j){  //hww_ee
    XDATA[j][NCHA+6] = (TH1F*) f2_d0[0]->Get("Data Final Variable (dphi) - "+tits[j]); 
    if(HMass[j] !=200){
    XHIGGS[j][NCHA+6] = (TH1F*) f2_d0[0]->Get("sig Interpolate "+tit[j]);
    XBKG_d0[0][j][6] = (TH1F*) f2_d0[0]->Get("wz Interpolate "+tit[j]);
    XBKG_d0[1][j][6] = (TH1F*) f2_d0[0]->Get("zz Interpolate "+tit[j]);
    XBKG_d0[2][j][6] = (TH1F*) f2_d0[0]->Get("ww Interpolate "+tit[j]);
    XBKG_d0[3][j][6] = (TH1F*) f2_d0[0]->Get("tt Interpolate "+tit[j]);
    XBKG_d0[4][j][6] = (TH1F*) f2_d0[0]->Get("wjet Interpolate "+tit[j]);
    XBKG_d0[5][j][6] = (TH1F*) f2_d0[0]->Get("ztt Interpolate "+tit[j]);
    XBKG_d0[6][j][6] = (TH1F*) f2_d0[0]->Get("zee Interpolate "+tit[j]);
    XBKG_d0[7][j][6] = (TH1F*) f2_d0[0]->Get("qcd Interpolate "+tit[j]);
    }
    else{ 
    XHIGGS[j][NCHA+6] = (TH1F*) f2_d0[0]->Get("Signal Final Variable (200 mH)"); 
    XBKG_d0[0][j][6] = (TH1F*) f2_d0[0]->Get("WZ-Bkgd Final Variable (dphi) - 200"); 
    XBKG_d0[1][j][6] = (TH1F*) f2_d0[0]->Get("ZZ-Bkgd Final Variable (dphi) - 200"); 
    XBKG_d0[2][j][6] = (TH1F*) f2_d0[0]->Get("WW-Bkgd Final Variable (dphi) - 200"); 
    XBKG_d0[3][j][6] = (TH1F*) f2_d0[0]->Get("tt-Bkgd Final Variable (dphi) - 200");
    XBKG_d0[4][j][6] = (TH1F*) f2_d0[0]->Get("Wjet-Bkgd Final Variable (dphi) - 200"); 
    XBKG_d0[5][j][6] = (TH1F*) f2_d0[0]->Get("ztt-Bkgd Final Variable (dphi) - 200"); 
    XBKG_d0[6][j][6] = (TH1F*) f2_d0[0]->Get("zee-Bkgd Final Variable (dphi) - 200"); 
    XBKG_d0[7][j][6] = (TH1F*) f2_d0[0]->Get("qcd-Bkgd Final Variable (dphi) - 200"); 
    }
  }

  for(int j=2; j<NMASS; ++j){  //hww_emu
    XDATA[j][NCHA+7] = (TH1F*) f2_d0[1]->Get("Data Final Variable (dphi) - "+tits[j]); 
    if(HMass[j] !=200){
    XHIGGS[j][NCHA+7] = (TH1F*) f2_d0[1]->Get("sig Interpolate "+tit[j]);
    XBKG_d0[0][j][7] = (TH1F*) f2_d0[1]->Get("wz Interpolate "+tit[j]);
    XBKG_d0[1][j][7] = (TH1F*) f2_d0[1]->Get("zz Interpolate "+tit[j]);
    XBKG_d0[2][j][7] = (TH1F*) f2_d0[1]->Get("ww Interpolate "+tit[j]);
    XBKG_d0[3][j][7] = (TH1F*) f2_d0[1]->Get("tt Interpolate "+tit[j]);
    XBKG_d0[4][j][7] = (TH1F*) f2_d0[1]->Get("wjet Interpolate "+tit[j]);
    XBKG_d0[5][j][7] = (TH1F*) f2_d0[1]->Get("ztt Interpolate "+tit[j]);
    XBKG_d0[6][j][7] = (TH1F*) f2_d0[1]->Get("zmm Interpolate "+tit[j]);
    XBKG_d0[7][j][7] = (TH1F*) f2_d0[1]->Get("qcd Interpolate "+tit[j]);
    }
    else{
    XHIGGS[j][NCHA+7] = (TH1F*) f2_d0[1]->Get("Signal Final Variable (200 mH)");
    XBKG_d0[0][j][7] = (TH1F*) f2_d0[1]->Get("WZ-Bkgd Final Variable (dphi) - 200"); 
    XBKG_d0[1][j][7] = (TH1F*) f2_d0[1]->Get("ZZ-Bkgd Final Variable (dphi) - 200"); 
    XBKG_d0[2][j][7] = (TH1F*) f2_d0[1]->Get("WW-Bkgd Final Variable (dphi) - 200"); 
    XBKG_d0[3][j][7] = (TH1F*) f2_d0[1]->Get("tt-Bkgd Final Variable (dphi) - 200"); 
    XBKG_d0[4][j][7] = (TH1F*) f2_d0[1]->Get("Wjet-Bkgd Final Variable (dphi) - 200"); 
    XBKG_d0[5][j][7] = (TH1F*) f2_d0[1]->Get("ztt-Bkgd Final Variable (dphi) - 200"); 
    XBKG_d0[6][j][7] = (TH1F*) f2_d0[1]->Get("zmm-Bkgd Final Variable (dphi) - 200"); 
    XBKG_d0[7][j][7] = (TH1F*) f2_d0[1]->Get("qcd-Bkgd Final Variable (dphi) - 200"); 
    }      
  }

  for(int j=2; j<NMASS; ++j){  //hww_mumu
    XDATA[j][NCHA+8] = (TH1F*) f2_d0[2]->Get("Data Final Variable (dphi)"); 
    XHIGGS[j][NCHA+8] = (TH1F*) f2_d0[2]->Get("Signal Interpolate "+tit[j]);
    if(j==(NMASS-1))XHIGGS[j][NCHA+8] = (TH1F*) f2_d0[2]->Get("Signal Final Variable (200 dphi)");
    XBKG_d0[0][j][8] = (TH1F*) f2_d0[2]->Get("WZ-Bkgd Final Variable (dphi)"); 
    XBKG_d0[1][j][8] = (TH1F*) f2_d0[2]->Get("Zbb-Bkgd Final Variable (dphi)"); 
    XBKG_d0[2][j][8] = (TH1F*) f2_d0[2]->Get("Zjj-Bkgd Final Variable (dphi)"); 
    XBKG_d0[3][j][8] = (TH1F*) f2_d0[2]->Get("Wtn-Bkgd Final Variable (dphi)"); 
    XBKG_d0[4][j][8] = (TH1F*) f2_d0[2]->Get("tt-Bkgd Final Variable (dphi)"); 
    XBKG_d0[5][j][8] = (TH1F*) f2_d0[2]->Get("qcd-Bkgd Final Variable (dphi)"); 
  }

  // 
  // what about wh->www->llnunu from d0
  TString titx[NMASS]={"100 + 10", "115 + 0","115 + 5","115 + 15","135 + 5", "135 + 15","155 + 5","155 + 15","175 + 5","175 + 25"}; 
  TString titxs[NMASS]={"110", "115","120","130","140","150","160","170","180","200"}; 
  for(int i =0; i<3; ++i){ // www_ee, www_emu, www_mumu 
  for(int j=0; j<NMASS; ++j){ 
    XDATA[j][NCHA+9+i] = (TH1F*) f3_d0[i]->Get("Data Final Variable (dphi) - "+titxs[j]); 
    XHIGGS[j][NCHA+9+i] = (TH1F*) f3_d0[i]->Get("sig Interpolate "+titx[j]);
    XHIGGS[j][NCHA+9+i] ->Rebin(4); 
    XBKG_d0[0][j][9+i] = (TH1F*) f3_d0[i]->Get("WZ Interpolate "+titx[j]);
    XBKG_d0[1][j][9+i] = (TH1F*) f3_d0[i]->Get("ZZ Interpolate "+titx[j]);
    XBKG_d0[2][j][9+i] = (TH1F*) f3_d0[i]->Get("flip Interpolate "+titx[j]);
    XBKG_d0[3][j][9+i] = (TH1F*) f3_d0[i]->Get("qcd Interpolate "+titx[j]);
  }
  }
  //
  // what about zh->llbb from d0
  TString tity[NMASS]={"105 + 5", "115 + 0","115 + 5","125 + 5","135 + 5", "145 + 5"," "," "," "," "}; 
  for(int i =0; i<2; ++i){ // eebb, mmbb
    XDATA[0][NCHA+12+i] = (TH1F*) f4_d0[i]->Get("Data Final Variable (mH)");
    XDATA[0][NCHA+12+i] ->Rebin(2); 
    for(int j = 0;  j<Nmax[NCHA+12+i]; ++j){ 
    XHIGGS[j][NCHA+12+i] = (TH1F*) f4_d0[i]->Get("hwlv Interpolate "+tity[j]);
    XHIGGS[j][NCHA+12+i] ->Rebin(2);
    }
    XBKG_d0[0][0][12+i] = (TH1F*) f4_d0[i]->Get("WZ-Bkgd Final Variable (mH)");
    XBKG_d0[1][0][12+i] = (TH1F*) f4_d0[i]->Get("ZZ-Bkgd Final Variable (mH)");
    XBKG_d0[2][0][12+i] = (TH1F*) f4_d0[i]->Get("Zbb-Bkgd Final Variable (mH)");
    XBKG_d0[3][0][12+i] = (TH1F*) f4_d0[i]->Get("Zjj-Bkgd Final Variable (mH)");
    XBKG_d0[4][0][12+i]  = (TH1F*) f4_d0[i]->Get("tt-Bkgd Final Variable (mH)");
    XBKG_d0[5][0][12+i]  = (TH1F*) f4_d0[i]->Get("qcd-Bkgd Final Variable (mH)");
    for(int k=0; k<6; k++)XBKG_d0[k][0][12+i]->Rebin(2);
 
    for(int j =1; j<Nmax[i+NCHA+12]; ++j){
      XDATA[j][i+NCHA+12] = XDATA[0][i+NCHA+12]; 
      for(int k=0; k<NBKGHD0-2; ++k){ 
	XBKG_d0[k][j][12+i] = XBKG_d0[k][0][12+i];
      }
    }
  }

  const Int_t nXsecR = 600;
  static Int_t nIter =3000000; //1000000; // 100k or 500k 
  if(NoSYS==3) nIter = nXsecR; 
  const float step =0.5; 
  // random number
  TRandom gRand;
  gRand.SetSeed();

  // cross section step size
  Int_t NMASSX(0); 
  Int_t nchannel(0); 
  for(Int_t i=0; i<NCHAT; ++i){ 
    if(icha[i]>0){
      nchannel++; 
      if(NMASSX<Nmax[i])NMASSX = Nmax[i]; 
    }
  }
   
  // combined chi^{2} and likelihood histograms
  
  TH1F *hchi[NMASS];
  TH1F *hlh [NMASS];
  Int_t nmbin = NMASS*2; 
  TH1F *hl = new TH1F("Limit","Limit", nmbin, 105., 105+NMASSX*10);
  for ( Int_t i = 0 ; i < NMASSX ; i++ ) {
      TString nchi("hchi_");
      int j = HMass[i]; 
      nchi += (Int_t)j;

      TString nlh("hlh_");
      nlh += (Int_t)j;

      hchi[i] = new TH1F(nchi,nchi,nXsecR,0.,step*nXsecR);
      hchi[i]->SetXTitle("#sigma Ratio");
      hchi[i]->SetYTitle("#chi^{2}");

      hlh[i]  = new TH1F(nlh,nlh,nXsecR,0.,step*nXsecR);
      hlh[i]->SetXTitle("#sigma Ratio");
      hlh[i]->SetYTitle("likelihood");
  }

  std::cout <<" Channels Combined="<<nchannel<< " " <<" nIter "<<nIter<<std::endl; 
  if(NoSYS==1||NoSYS==2||NoSYS==3) std::cout<<" With no systematic included "<<" "<<NoSYS<<std::endl; 
  // 95% C.L. upper limit
  Double_t fit[NMASS][3] = {0.}; // central value, +sigma, -sigma
  Double_t limit[NMASS] = {0.};

  // start mass loop
  for ( Int_t im = 0 ; im < NMASSX ; im++ ) {
    if(NoSYS>100&&HMass[im]!=NoSYS)continue;
    // start sampling loop
    Double_t lh[nXsecR] = {0.};
    Double_t nc[nXsecR] =  {0.};

    for ( Int_t Iter = 0 ; Iter < nIter ; Iter++ ) {
      // get fluctuation: for eff
      if(Iter%100000==0)std::cout<<" Processed so far Iter "<<Iter<<std::endl; 
      float gcomm[NCOMM]; // common systematic between cdf and d0 
      float gcommb_lum; // allow to turn off common systematic between signal and background 

      float geff_d0[NACCD0][NCHAD0]; // d0 acceptances of systematic 
      float gbkg_d0[NBKGD0][NCHAD0]; // d0 backgrounds of systematic 
      float geff[NACC][NCHA]; //cdf acceptance of systematic 
      float geffb[NACC][NCHA]; // possible to break down the signal and backgrounds 
      float gbkg[NBKG][NCHA]; // cdf backgrounds of systematic 
      float gjes_nnbb[2][NBKGH]; // for nunubb jes on the backgrounds  
      float gtagc; // tagging c probability at cdf 
      // modeling efficiencies 
      for(int i = 0; i<NCOMM; ++i){ 
	gcomm[i] = getGauss(gRand, 1., fabs(scomm[i]) );
      }
      if(NoSYS==4){
	gcommb_lum = getGauss(gRand, 1., fabs(scomm[0]) );
      }
      else{ 
	gcommb_lum = gcomm[0]; 
      }
      
      for(int i =0; i<NACCD0; ++i){  // for D0 treat the signal and background correlated 
	float val; 
	do{
	  float a = gRand.Gaus(0.,1.0);
	  float ab = a; 
	  if(NoSYS==4) ab =gRand.Gaus(0.,1.0);
	  val =1.0; 
	  for(int j =0; j<NCHAD0; ++j){ 
	      geff_d0[i][j] =(1.0+a*seff_d0[i][j]);
	      gbkg_d0[i][j] =(1.0+ab*sbkg_d0[i][im][j]); 
	      if(geff_d0[i][j]<=0.)val = geff_d0[i][j];
	      if(gbkg_d0[i][j]<=0.)val = gbkg_d0[i][j]; 
	  }
	}while(val<=0.); 
      }
	
      for(int i =NACCD0; i<NBKGD0; ++i){ 
	float val; 
	do{
	  float a = gRand.Gaus(0.,1.0);
	  val =1.0; 
	  for(int j =0; j<NCHAD0; ++j){ 
	      gbkg_d0[i][j] =(1.0+a*sbkg_d0[i][im][j]);
	      if(gbkg_d0[i][j]<=0.)val = gbkg_d0[i][j]; 
	  }
	}while(val<=0.); 
      }

      // modeling efficiency >0 for CDF 
      for(int i =0; i<NACC; ++i){  // for CDF treat the signal correlated 
	float val; 
	do{
	  float a = gRand.Gaus(0.,1.0);
	  float ab = a; 
	  if(NoSYS==4)ab=gRand.Gaus(0.,1.0);
	  val =1.0; 
	  for(int j =0; j<NCHA; ++j){ 
	    if(seff[i][j]>0){  // some problems here 
		geff[i][j] =(1.0+a*seff[i][j]);
		geffb[i][j] =(1.0+ab*seff[i][j]);
	      }
	      else{
		geff[i][j] = getGauss(gRand, 1., fabs(seff[i][j]));
		geffb[i][j] = geff[i][j]; 
		if(NoSYS==4) geffb[i][j] = getGauss(gRand, 1., fabs(seff[i][j]));
	      }
	      if(geff[i][j]<=0.)val = geff[i][j];
	      if(geffb[i][j]<=0.)val = geffb[i][j];
	  }
	  if(i==1){ 
	    gtagc = 1.0+0.16*ab; 
	    if(gtagc<=0.)val = gtagc; 
	  }
	  if(i==3){ 
	    for(int j=0; j<5; ++j){ // jes uncertainties for nunubb 
	      if(sjes_nnbb[0][j]>0){
		gjes_nnbb[0][j] = (1.0 +a*sjes_nnbb[0][j]); 
	      }
	      else{ 
		gjes_nnbb[0][j] = getGauss(gRand, 1.0, fabs(sjes_nnbb[0][j])); 
	      }
              if(sjes_nnbb[1][j]>0){
		gjes_nnbb[1][j] = (1.0 +a*sjes_nnbb[1][j]);
	      }
	      else{
		gjes_nnbb[1][j] = getGauss(gRand, 1.0, fabs(sjes_nnbb[1][j])); 
	      }
	      if(gjes_nnbb[0][j]<=0.)val = gjes_nnbb[0][j];
	      if(gjes_nnbb[1][j]<=0.)val = gjes_nnbb[1][j]; 
	    }
	  }
	}while(val<=0.); 
      }
	
      // modeling background >0 
      // for CDF treat the background systematic and including some of systematic from acceptance too
      // For example some of uncertainties are partial only in nunubb and zh->llbb 
      for(int i =0; i<NBKG; ++i){  
	float val; 
	do{
	  float a = gRand.Gaus(0.,1.0);
	  val =1.0; 
	  for(int j =0; j<NCHA; ++j){ 
	      if(sbkg[i][im][j]>0){
		gbkg[i][j] =(1.0+a*sbkg[i][im][j]);
	      }
	      else{
		gbkg[i][j] = getGauss(gRand, 1., fabs(sbkg[i][im][j]));
	      }
	      if(gbkg[i][j]<=0.)val = gbkg[i][j];
	  }
	}while(val<=0.); 
      }

      // additional smear of ww xsec 
      float xww = getGauss(gRand,1.0, 0.1); //assign additional 10% 
      if(NoSYS==3||NoSYS==1)xww =1.; 
      
      // shape sampling  between two with prob = 0,1
      float xsampling = gRand.Uniform(); 
      // get cross section Ratio from 0 to nXsecR*step  
      Double_t tXsec = gRand.Uniform()*nXsecR*step; 
      Int_t tBin = (Int_t)(tXsec/step);
      if(NoSYS==3){
	tXsec =(Iter+0.5)*step; 
	tBin = Iter; 
      }
      nc[tBin]++;

      // likelihood calculation
      Double_t loglhS = 0.;
      Double_t efftot = 0.; 
      for( Int_t ich=0; ich<NCHAT; ++ich){ 
	if(icha[ich]>0&&im<Nmax[ich]&&(ich<(NCHAT-8)||(ich>=(NCHAT-8)&&ich<(NCHAT-5)&&im>1)||ich>=(NCHAT-5))){ // there is no 110 for D0 ww analysis 
	  if(ich<NCHA-1){ // cdf channels 
	  // eff 
	  Double_t acc(1.); 
	  for( Int_t i=0; i<NACC; ++i) {
	    acc *=geff[i][ich]; 
	  }
	  acc *=gcomm[0]; // commom tev luminosity 
	  acc *=sfnnlo[ich]; // correct the nnlo factor for Higgs normalization 
	  if(txt[ich]=="ww")acc *=xww; // assign additional 10% uncertainties 
	  float xb[NBKGH]={0.}; // background normalization
	  if(ich<4){ //lvbb_s, lvbb_d, nnbb, ww 
	  for( Int_t j=0; j<NBKGH-1; ++j){ 	   
	    if(fabs(sbkg[j][im][ich])>0.)xb[j] =gbkg[j][ich]; 
	    if(j==3) xb[j] *=gcomm[3]*gcommb_lum*geffb[0][ich]; //common top xsec 
	    if(j==4) xb[j] *=gcomm[4]*gcommb_lum*geffb[0][ich]; // common diboson 
	  }
	  if(fabs(sbkg[NBKG-1][im][ich])>0.)xb[NBKGH-1] = gbkg[NBKG-1][ich]; // special for ww 
	  }
	  else if(ich<6){ // nunubb_s, nunubb_d

	  for( Int_t j=0; j<NBKGH; ++j){ 	   
	    if(j<2&&fabs(sbkg[j][im][ich])>0.){
	      xb[j] =gbkg[j][ich]; 
	    }
	    else if(j<(NBKGH-1)&&(fabs(sbkg[j][im][ich])>0.||fabs(sbkg[j+3][im][ich])>0.)){ // correlated or uncorrected present 
	      xb[j]=gbkg[j][ich]*gbkg[j+3][ich]; //w+hf
	    }
	    else if(j==(NBKGH-1)&&(fabs(sbkg[2][im][ich])>0.||fabs(sbkg[j+3][im][ich])>0.)){
	      xb[j] = gbkg[2][ich]*gbkg[j+3][ich]; //z+hf
	    }
	    if(j==3) xb[j] *=gcomm[3]; //common top xsec 
	    if(j==4) xb[j] *=gcomm[4]; // common diboson
	    float aj = gjes_nnbb[0][j]; 
	    if(ich==5)aj=gjes_nnbb[1][j]; 
	    if(j>1)xb[j] *=aj*gcommb_lum*geffb[0][ich]*geffb[5][ich]*geffb[1][ich]*geffb[2][ich]*geffb[4][ich]; //no addition sys for mistag/qcd
	  }

	  }
	  // No SYS applied here 
	    if(NoSYS==3||NoSYS==1)acc =sfnnlo[ich]; 
	    if(NoSYS==3||NoSYS==2){
	    for(Int_t j =0; j<NBKGH; ++j){
	      if(xb[j]>0)xb[j] = 1.;
	    }
	    }

	  for(Int_t k =1; k<XDATA[im][ich]->GetNbinsX()+1; ++k){ 

	    Double_t mc = 0.; 
	    mc += tXsec*acc*XHIGGS[im][ich]->GetBinContent(k); 
	    efftot +=acc*XHIGGS[im][ich]->GetBinContent(k);

	    //background 
	    for( Int_t j=0; j<NBKGH; ++j){ 
	      if(xb[j]>0.)mc +=xb[j]*XBKG[j][im][ich]->GetBinContent(k); // xb[j]==0, the hist is not defined
	    }
	    Int_t nx = (Int_t) XDATA[im][ich]->GetBinContent(k); 
	    if(mc>0) loglhS +=XDATA[im][ich]->GetBinContent(k)*TMath::Log(mc)-mc-lgm[nx]; 
	  } //loop over the histogram bins 
	  } // for cdf channels
	  else if(ich==(NCHA-1)){
            //*************************************************** 
	    // for special treat of 2-d llbb fit from cdf 
	    //***************************************************
	    // eff 
	    Double_t acc(1.); 
	    for( Int_t i=0; i<NACC; ++i) {
	      acc *=geff[i][ich]; 
	    }
	    acc *=gcomm[0]; // commom tev luminosity 
	    acc *=sfnnlo[ich]; // correct the nnlo factor for Higgs normalization 
	    if(txt[ich]=="ww")acc *=xww; // assign additional 10% uncertainties 
	    float xb[NBKGH]={0.}; // background normalization
	    for( Int_t j=0; j<NBKGH; ++j){ 	   
	      
	    if(j<(NBKGH-1)&&fabs(sbkg[j][im][ich])>0.){
	      xb[j] =gbkg[j][ich];
	      if(j>1)xb[j] *=geffb[1][ich]; // btag
	    }
	    else if(j==(NBKGH-1)&&fabs(sbkg[2][im][ich])>0.){ 
	      xb[j] =gbkg[2][ich]*gtagc; //ctag
	    }
	    if(j==3) xb[j] *=gcomm[3]; //common top xsec 
	    if(j==4) xb[j] *=gcomm[4]; // common diboson 
	    if(j>1)xb[j] *=gcommb_lum*geffb[0][ich]*geffb[2][ich]*geffb[3][ich]*geffb[4][ich]; // add additional systematic:lum,lep id, jes, isr/fsr+pdf   
	    }
	    // No SYS applied here 
	    if(NoSYS==3||NoSYS==1)acc =sfnnlo[ich]; 
	    if(NoSYS==3||NoSYS==2){
	      for(Int_t j =0; j<NBKGH; ++j){
		if(xb[j]>0)xb[j] = 1.;
	      }
	    }
	    // convert 1-d to 2d; 
	    TH2F* xd2d =(TH2F*)XDATA[im][ich]; 
	    TH2F* xh2d =(TH2F*)XHIGGS[im][ich]; 
	    TH2F* xb2d[NBKGH]; 
	    for(Int_t i=0; i<NBKGH; ++i) xb2d[i] =(TH2F*)XBKG[i][im][ich]; 

	    for(Int_t k =1; k<xd2d->GetNbinsX()+1; ++k){ 
	      for(Int_t k1=1; k1<xd2d->GetNbinsY()+1; ++k1){  
		
		Double_t mc = 0.; 
		mc += tXsec*acc*xh2d->GetBinContent(k,k1); 
		efftot +=acc*xh2d->GetBinContent(k,k1);

		//background 
		for( Int_t j=0; j<NBKGH; ++j){ 
		  if(xb[j]>0.){
		    float a1 = xb2d[j]->GetBinContent(k,k1); 
		    if(NoSYS==11){ 
		      if(j==2)a1 =a1*xsampling +(1-xsampling)*XZHshape[0]->GetBinContent(k,k1); 
		      if(j==5)a1 =a1*xsampling+(1-xsampling)*XZHshape[1]->GetBinContent(k,k1); 
		    }
		    mc +=xb[j]*a1; // xb[j]==0, the hist is not defined
		  }
		}
		Int_t nx = (Int_t) xd2d->GetBinContent(k,k1); 
		if(mc>0) loglhS +=xd2d->GetBinContent(k,k1)*TMath::Log(mc)-mc-lgm[nx]; 
	      } //loop over the histogram bins x 
	    } //
	  }
	  //*****************************************************************************
          // liklihood for D0 channels 
          //*****************************************************************************
	  else{ // for d0 channels 
	  // eff 

	  Double_t acc(1.); 
	  
	  for( Int_t i=0; i<NACCD0; ++i) {
	    acc *=geff_d0[i][ich-NCHA];
	      }
	  acc *=gcomm[0]; // common luminosity 
	  if(txt[ich].Contains("d0ww"))acc *=xww; // assign additional 10% uncertainties 

	  Double_t acc_bkg(1.0); 
	  for( Int_t i=0; i<NACCD0; ++i) {
	    acc_bkg *=gbkg_d0[i][ich-NCHA];
	      }
	  acc_bkg *=gcommb_lum; // common luminosity 
	  
	  float xb[NBKGHD0]={0.}; // background normalization 
	  // signal id
	  float eids=geff_d0[8][ich-NCHA]*geff_d0[9][ich-NCHA]*geff_d0[10][ich-NCHA]; 
	  float mids=geff_d0[12][ich-NCHA]*geff_d0[13][ich-NCHA];
	  float bids = geff_d0[5][ich-NCHA]*geff_d0[6][ich-NCHA]; 
	  // background id
	  float eid=gbkg_d0[8][ich-NCHA]*gbkg_d0[9][ich-NCHA]*gbkg_d0[10][ich-NCHA]; 
	  float mid=gbkg_d0[12][ich-NCHA]*gbkg_d0[13][ich-NCHA];
          float bid =gbkg_d0[5][ich-NCHA]*gbkg_d0[6][ich-NCHA]; 	   
	  // dileptons extra id (ee,mm), double btags; ignore for now by reset to 1.  
	  eids =1.; 
	  mids =1.; 
	  bids =1.; 
	  eid = 1.; 
	  mid = 1.; 
	  bid = 1.; 
	  if(txt[ich].Contains("d0wwee")){
	    acc *=eids;
	    for( Int_t j=0; j<NBKGHD0; ++j){ 	   
	      if(j <3){
		xb[j] = eid*acc_bkg*gcomm[4]; //*Xsec_EM
	      }
	      else if( j==3){
		xb[j] = eid*acc_bkg*gcomm[3]; // *Xsec_top
	      }
	      else if(j<NBKGHD0-1){
		xb[j] = eid*acc_bkg*gcomm[4]; //*Xsec_EM 		
	      }
	      else{
		xb[j] = eid*acc_bkg*gbkg_d0[22][ich-NCHA]; //*Xsec_qcd 	
	      }
	    }	    
	  }
	  else if(txt[ich].Contains("d0wwem")){
	    for( Int_t j=0; j<NBKGHD0; ++j){ 	   
	      if(j <3){
		xb[j] = acc_bkg*gcomm[4]; //*Xsec_EM
	      }
	      else if( j==3){
		xb[j] = acc_bkg*gcomm[3]; // *Xsec_top
	      }
	      else if(j<NBKGHD0-1){ 
		xb[j] = acc_bkg*gcomm[4]; //*Xsec_EM 	
	      }
	      else{
		xb[j] = acc_bkg*gbkg_d0[22][ich-NCHA]; //*Xsec_qcd 
	      }
	    }	    
	  }
	  else if(txt[ich].Contains("d0wwmm")){
	    acc *=mids;
	    for( Int_t j=0; j<NBKGHD0-2; ++j){ 	   
	      if(j==0||j==3){
		xb[j] = mid*acc_bkg*gcomm[4]; //*Xsec_EM
	      }
	      else if(j==1){ 
		xb[j] = mid*acc_bkg*gbkg_d0[20][ich-NCHA]; //*Xsec_HF
	      }
	      else if(j==2){ 
		xb[j] = mid*acc_bkg*gbkg_d0[21][ich-NCHA]; //*Xsec_LF
	      }
              else if(j==4){ 
		xb[j] = mid*acc_bkg*gcomm[3]; // *Xsec_top
	      }
	      else{
		xb[j] = mid*acc_bkg*gbkg_d0[22][ich-NCHA]; //*Xsec_qcd 	
	      }
	    }	    
	  }
	  else if(txt[ich].Contains("d0wh")){ // wh->www same sign 
	    for( Int_t j=0; j<4; ++j){ 	   
	      if(j<2){
		xb[j] = acc_bkg*gcomm[4]; //*Xsec_EM
	      }
	      else{
		xb[j] = acc_bkg*gbkg_d0[22][ich-NCHA]; //*Xsec_qcd 	
	      }
	    }	    
	  }
	  else if((txt[ich]=="d0zh")){ // d0zheebb or d0zhmmbb 2tags
	    for( Int_t j=0; j<NBKGHD0-2; ++j){ 	   
	      if(j <2){
		xb[j] = acc_bkg*gcomm[4]; //eid*Lum*trigger*Xsec_EM
	      }
	      else if(j<3){
		xb[j] = acc_bkg*gbkg_d0[20][ich-NCHA]; //*Xsec_HF
	      }
	      else if( j<4){
		xb[j] = acc_bkg*gbkg_d0[19][ich-NCHA]; //*Xsec_bTag_LF
	      }
	      else if( j<5){
		xb[j] = acc_bkg*gcomm[3]; // *Xsec_top
	      }
	      else{ 
		xb[j] = acc_bkg*gbkg_d0[22][ich-NCHA]; //*Xsec_qcd 	
	      }	      	    
	    }
	  }  
	  else if((txt[ich]=="d0evbb_s")||(txt[ich]=="d0mvbb_s")||(txt[ich]=="d0vvbb_s")){ // single tag
	    for( Int_t j=0; j<NBKGHD0-1; ++j){ 	   
	      if(j <1){
		xb[j] = acc_bkg*gcomm[4]; //eid*Lum*trigger*Xsec_EM
	      }
	      else if(j<2){
		xb[j] = acc_bkg*gbkg_d0[20][ich-NCHA]; //*Xsec_HF
	      }
	      else if( j<4){
		xb[j] = acc_bkg*gbkg_d0[19][ich-NCHA]; //*Xsec_bTag_LF
	      }
	      else if( j<6){
		xb[j] = acc_bkg*gcomm[3]; // *Xsec_top
	      }
	      else{ 
		xb[j] = acc_bkg*gbkg_d0[22][ich-NCHA]; //*Xsec_qcd 	
	      }	      	    
	    }
	  }  
	  else{ // double tags 
	    acc *=bids; 
	    for( Int_t j=0; j<NBKGHD0-1; ++j){ 	   
	      if(j <1){
		xb[j] = bid*acc_bkg*gcomm[4]; //*Xsec_EM
	      }
	      else if( j<2){ 
		xb[j] = bid*acc_bkg*gbkg_d0[20][ich-NCHA]; //*Xsec_HF
	      }
	      else if( j<4){
		xb[j] = bid*acc_bkg*gbkg_d0[19][ich-NCHA]; //*Xsec_bTag_LF
	      }
	      else if( j<6){
		xb[j] = bid*acc_bkg*gcomm[3]; // *Xsec_top
	      }
	      else{ 
		xb[j] = bid*acc_bkg*gbkg_d0[22][ich-NCHA]; //*Xsec_qcd 	
	      }	      	    
	    }
	  }  
	  //		
	  // No SYS applied here 
	    if(NoSYS==3||NoSYS==1)acc =sfnnlo[ich]; 
	    if(NoSYS==3||NoSYS==2){ 
	    for(Int_t j =0; j<NBKGHD0; ++j){
	      if(xb[j]>0)xb[j] = 1.;
	    }
	    }
	  for(Int_t k =1; k<XDATA[im][ich]->GetNbinsX()+1; ++k){ 

	    //if(txt[ich].Contains("d0wwmm")&&XDATA[im][ich]->GetBinCenter(k)>2.5)continue; // requiring dphi<2.5 for ww->mm channel to remove Z->mumu events 
	    Double_t mc = 0.; 
	    mc += tXsec*acc*XHIGGS[im][ich]->GetBinContent(k); 
	    efftot +=acc*XHIGGS[im][ich]->GetBinContent(k);

	    //background 
	    Int_t nx = NBKGHD0-1; 
	    if(txt[ich].Contains("d0wwmm")){
	      nx = NBKGHD0-2; 
	    }
	    else if(txt[ich].Contains("d0ww")){
	      nx = NBKGHD0;
	    }
	    if(txt[ich].Contains("d0wh"))nx = NBKGHD0-4; 
	    if(txt[ich].Contains("d0zh"))nx = NBKGHD0-2;
	    for( Int_t j=0; j<nx; ++j){
	      if(xb[j]>0.)mc +=xb[j]*XBKG_d0[j][im][ich-NCHA]->GetBinContent(k); // xb[j]==0, the hist is not defined
	    }
	    
	    nx = (Int_t) XDATA[im][ich]->GetBinContent(k); 
	    if(mc>0)loglhS +=XDATA[im][ich]->GetBinContent(k)*TMath::Log(mc)-mc-lgm[nx]; 
	    
	  } //loop over the histogram bins
	  }
	} // channels in consideration 	
      } // loop all the channels
      lh[tBin] +=efftot*TMath::Exp(loglhS);
    }   // end of sampling loop

    // average of likelihood
    Double_t avelh [nXsecR] = {0.};
    Double_t sumlh(0.); 
    for ( Int_t i = 0 ; i < nXsecR ; i++ ) {
      if ( nc[i] > 0. ) {
	avelh[i] = lh[i]/nc[i];
	sumlh +=avelh[i]; 
      }
    }
    // normalize
    Double_t maxlh(0.);
    Double_t avechi[nXsecR] = {0.};
    Int_t XsecMax(0); 
    for ( Int_t Xsec = 0 ; Xsec < nXsecR ; Xsec++ ) {
      avelh[Xsec] /= sumlh;
      if ( avelh[Xsec] > 0. ) avechi[Xsec] = -2.*TMath::Log(avelh[Xsec]);
      if ( maxlh  < avelh[Xsec] ) {
	fit[im][0] = (Xsec+0.5)*step;
	maxlh  = avelh[Xsec];
	XsecMax = Xsec; 
      }

      hchi[im]->Fill(Xsec*step+step/2,avechi[Xsec]);
      hlh [im]->Fill(Xsec*step+step/2,avelh [Xsec]);
    }
    // compute the +- sigma
    Int_t di(0); 
    for(int i =XsecMax; i<nXsecR; ++i) {
      if((avechi[i]-avechi[XsecMax])>1.0){ 	
	di = i;
	break; 
      }
    }
    if(di==0) di = nXsecR; 
    fit[im][1] = (di-XsecMax)*step; 
    di=0;
    for(int i =XsecMax; i>-1; --i) { 
      if((avechi[i]-avechi[XsecMax])>1.0){ 
	di = i;
	break; 
      }
    }
    fit[im][2] = (XsecMax-di)*step; 

    // 95% C.L. upper limit
    Double_t A = 0.;
    for ( Int_t Xsec = 0 ; Xsec < nXsecR ; Xsec++ ) {
      A += avelh[Xsec];
      if ( A > 0.95 ) {
	limit[im] = Xsec*step+step/2;
	break;
      }
    }

    std::cout<<" Mass="<< HMass[im]
	     <<" XsecR="<< fit  [im][0]
	     <<"+"<< fit  [im][1]
	     <<"-"<< fit  [im][2]	
	     << " 95% CL limit ="<<limit[im]
	      << std::endl;
  }

  // set 95% C.L. upper limit line for likelihood
 
  for(Int_t i=0; i<NMASSX; ++i){ 
    Int_t ib = hl->GetXaxis()->FindBin(HMass[i]); 
    hl->SetBinContent(ib, limit[i]); 
    hl->SetBinError(ib,0.); 
  }

  TLine line[NMASS];
  for ( Int_t i = 0 ; i < NMASSX ; i++ ) {
      line[i].SetX1(limit[i]);
      line[i].SetY1(0.);
      line[i].SetX2(limit[i]);
      line[i].SetY2(hlh[i]->GetMaximum()*0.2);
      line[i].SetLineColor(2);
  }

  TLatex pPreliminary(0.20,0.92,"CDF D0 Run II Preliminary (1 fb^{-1})");
  pPreliminary.SetNDC(true);
  pPreliminary.SetTextSize(0.05);

  // draw histograms
  TString mytxt("_"); 
  if(nchannel==NCHAT){ 
    mytxt += "all"; }
  else{
    for(Int_t i=0; i<NCHAT; ++i){ 
      if(icha[i]>0){
	mytxt +="1";
      }
      else{
	mytxt +="0";
      }
    }
  }


  TString psFile("R_combined_limit_cdfd0_new");
  if(NoSYS==3)psFile +="_NoSys";
  if(NoSYS==1)psFile +="_NoSysEff";
  if(NoSYS==2)psFile +="_NoSysBkg";
  if(NoSYS==4)psFile +="_NoCor"; // no correlation between signal and backgrounds
  if(NoSYS==10)psFile +="_shape"; // ZH->llbb with different zh shape 
  if(NoSYS==11)psFile +="_combshape"; // ZH->llbb with different zh shape with alternative combined shape 
  if(NoSYS==20)psFile +="_uncorr"; // treat most of systematic uncorrelated, except lum, btag, top ew xsec   
  if(NoSYS>100){
    psFile +="_mass";
    psFile +=(Int_t)NoSYS; 
  }
  psFile +=mytxt; 
  psFile += ".ps";

  TCanvas cx("cx","cx",640,640);
  TPostScript ps(psFile);

  TLatex nLimit[NMASS];
  for ( Int_t i = 0 ; i < NMASSX ; i++ ) {
      Char_t cLimit[200];
      Int_t Mass = HMass[i];
      sprintf(cLimit,"%3i GeV/c^{2} Limit = %3.1f pb",Mass,limit[i]);

      nLimit[i].SetNDC(true);
      nLimit[i].SetTextSize(0.03);
      nLimit[i].SetText(0.48,0.80-i*0.05,cLimit);
  }

    cx.Clear();
    Int_t nx0 =(NMASSX+1)/2; 
    cx.Divide(nx0,2);
    for ( Int_t i = 0 ; i < NMASSX ; i++ ) {
      cx.cd(i+1);
      hchi[i]->SetLineColor(1);
      hchi[i]->Draw();
    }
    cx.Update();
 

    cx.Clear();
    cx.Divide(nx0,2);
    for ( Int_t i = 0 ; i < NMASSX ; i++ ) {
      cx.cd(i+1);
      hlh[i]->SetLineColor(1);
      float ax =hlh[i]->GetXaxis()->GetXmax(); 
      if(ax>(2.0*limit[i]))ax = 2.0*limit[i]; 
      hlh[i]->GetXaxis()->SetRangeUser(0.,ax); 
      hlh[i]->Draw();
      line[i].Draw();
    }
    cx.Update();
    
    cx.Clear(); 
    cx.SetLogy();
    hl->SetMaximum(1000.); 
    hl->SetMinimum(1.); 
    hl->SetLineColor(2); 
    hl->SetLineWidth(4); 
    hl->Draw("AC");
    cx.Update(); 

  ps.Close();
  cx.Close();

  // output data file
  TString dataFile("R_combined_limit_cdfd0_new");
  if(NoSYS==3)dataFile +="_NoSys"; 
  if(NoSYS==1)dataFile +="_NoSysEff"; 
  if(NoSYS==2)dataFile +="_NoSysBKG"; 
  if(NoSYS==4)dataFile +="_NoCor"; // no correlation between signal and backgrounds
  if(NoSYS==10)dataFile +="_shape"; // ZH->llbb with different zh shape 
  if(NoSYS==11)dataFile +="_combshape"; // ZH->llbb with different zh shape with alternative combined shape 
  if(NoSYS==20)dataFile +="_uncorr"; // treat most of systematic uncorrelated, except lum, btag, top ew xsec   
  if(NoSYS>100){
    dataFile +="_mass";
    dataFile +=(Int_t)NoSYS;
  }

  dataFile += mytxt; 
  dataFile += ".dat";
  ofstream fout(dataFile);
  fout <<" Channels Combined="<<nchannel<<std::endl; 
  if(NoSYS==1||NoSYS==2||NoSYS==3)fout<<" with no systematic included "<< NoSYS<<std::endl; 

  for ( Int_t i = 0 ; i < NMASSX ; i++ ) {

    fout <<" Mass="<< HMass[i]
	 <<" XsecR="<< fit  [i][0]
	 <<"+"<< fit  [i][1]
	 <<"-"<< fit  [i][2]	
	 << " 95% CL limit="<<limit[i]
	 << std::endl;
  }
  fout.close();

  // output root file
  TString rootFile("R_combined_limit_cdfd0_new");
  if(NoSYS==3)rootFile +="_NoSys"; 
  if(NoSYS==1)rootFile +="_NoSysEff"; 
  if(NoSYS==2)rootFile +="_NoSysBkg"; 
  if(NoSYS==4)rootFile +="_NoCor"; // no correlation between signal and backgrounds
  if(NoSYS==10)rootFile +="_shape"; // ZH->llbb with different zh shape 
  if(NoSYS==11)rootFile +="_combshape"; // ZH->llbb with different zh shape with alternative combined shape 
  if(NoSYS==20)rootFile +="_uncorr"; // treat most of systematic uncorrelated, except lum, btag, top ew xsec   
  if(NoSYS>100){
    rootFile +="_mass";
    rootFile +=(Int_t)NoSYS;
  }
  rootFile += mytxt; 
  rootFile += ".root";
  TFile froot(rootFile,"recreate");
  for ( Int_t i = 0 ; i < NMASSX ; i++ ) {
      hchi[i]->Write();
      hlh [i]->Write();
  }
  hl->Write(); 

  froot.Close();

  // delete histograms
  for ( Int_t i = 0 ; i < NMASSX ; i++ ) {
      hchi[i]->Delete();
      hlh [i] ->Delete();
  }
  hl->Delete(); 

}   

//_________________________________________________________________________________
// Making pseudo-experiments for given channels defined in higgs_combined_limit.inp
//---------------------------------------------------------------------------------

void higgs_combined_limit_pseudoexperiment_cdfd0(Int_t NoSYS )
{
  // Input files: higgs_combined_limit.inp and many various root files needed 
  // Output the plot and limits
  gROOT->Time();
  gStyle->SetOptStat(0);
  gStyle->SetOptTitle(0);
  //
  TString txt[NCHAT] ={"lnubb_s","lnubb_d","nunubb","ww","nunubb_s", "nunubb_d", "llbb","d0evbb_s","d0evbb_d","d0mvbb_s","d0mvbb_d","d0vvbb_s","d0vvbb_d","d0wwee","d0wwem","d0wwmm","d0whee","d0whem","d0whmm","d0zhee","d0zhmm"};

  Double_t nfact(0); 
  lgm[0] = 0.; 
  for(Int_t i =1; i<200; ++i){ 
    nfact +=TMath::Log(i); 
    lgm[i] = nfact; 
  }
  char junk[120] =" ";
  FILE * infile = fopen("higgs_combined_limit.inp","r"); 
  for(int i=0; i<NCHA; ++i) {
    fscanf(infile, "%i", &icha[i]);
    if(debug)std::cout<<" icha[ "<< i<<" ]= "<<icha[i]<<std::endl; 
  }

  fgets(junk,60, infile);
  for(int j=0; j<NACC; ++j){
  for(int i=0; i<NCHA; ++i){
    fscanf(infile, "%f", &seff[j][i]); 
    if(debug)std::cout<<" seff[ "<< j<<" ]["<<i<<"] ="<<seff[j][i]<<std::endl; 
  }
  fgets(junk,60, infile);
  }

  for(int j=0; j<NBKG-1; ++j){
  for(int i=0; i<NCHA; ++i) {
    fscanf(infile, "%f", &sbkg[j][0][i]); 
    if(debug)std::cout<< " sbkg["<<j<<"][0]["<<i<<"]="<<sbkg[j][0][i]<<std::endl; 
  }
  fgets(junk,60, infile);
  }
  
  for(int i=0; i<NMASS; ++i) {
    fscanf(infile, "%f", &sbkg[NBKG-1][i][3]); // ww chid=3
    if(debug)std::cout<<" sbkg["<<NBKG-1<<"]["<<i<<"]["<<3<<"]="<<sbkg[NBKG-1][i][3]<<std::endl; 
  }
  fgets(junk,60, infile); 
  for(int i=0; i<NBKGH; ++i){ 
    fscanf(infile, "%f", &sjes_nnbb[0][i]); 
    if(debug)std::cout<<" sjes_nnbb[0]["<<i<<"]="<<sjes_nnbb[0][i]<<std::endl; 
  }

  fgets(junk,60, infile); 
  for(int i=0; i<NBKGH; ++i){ 
    fscanf(infile, "%f", &sjes_nnbb[1][i]); 
    if(debug)std::cout<<" sjes_nnbb[1]["<<i<<"]="<<sjes_nnbb[1][i]<<std::endl; 
  }

  for(int j=0; j<NBKG-1; ++j){
    for(int i=0; i<NCHA; ++i){ 
      for(int k=1; k<NMASS; ++k){ 
	sbkg[j][k][i] = sbkg[j][0][i]; 
      }
    }
  }
  // read in d0 inputs 
  infile = fopen("dumpD0_data.inp","r"); 
  for(int i=0; i<NCHAD0; ++i) {
    fscanf(infile, "%i", &icha[NCHA+i]);
    if(debug)std::cout<<" icha[ "<< NCHA+i<<" ]= "<<icha[NCHA+i]<<std::endl; 
  }
  fgets(junk,60, infile);
  for(int j=0; j<NCOMM; ++j){
    fscanf(infile, "%f", &scomm[j]); 
    if(debug)std::cout<<" scomm[ "<< j<<"] ="<<scomm[j]<<std::endl; 
  fgets(junk,110, infile);
  }

  for(int j=0; j<NACCD0; ++j){
  for(int i=0; i<NCHAD0; ++i){
    fscanf(infile, "%f", &seff_d0[j][i]); 
    if(debug)std::cout<<" seff_d0[ "<< j<<" ]["<<i<<"] ="<<seff_d0[j][i]<<std::endl; 
  }
  fgets(junk,60, infile);
  }

  for(int j=0; j<NBKGD0; ++j){
  for(int i=0; i<NCHAD0; ++i) {
    fscanf(infile, "%f", &sbkg_d0[j][0][i]); 
    if(debug)std::cout<< " sbkg_d0["<<j<<"][0]["<<i<<"]="<<sbkg_d0[j][0][i]<<std::endl; 
  }
  fgets(junk,60, infile);
  }
  
  for(int j=0; j<NBKGD0; ++j){
    for(int i=0; i<NCHAD0; ++i){ 
      for(int k=1; k<NMASS; ++k){ 
	sbkg_d0[j][k][i] = sbkg_d0[j][0][i]; 
      }
    }
  }

  // correct the tevatron common systematic for cdf 
  for(int i =0; i<NCHA; ++i){ 
    seff[0][i] = sqrt(pow(seff[0][i],2)-pow(scomm[0],2)); 
  }

  for(int i=0; i<NCHA; ++i){ // top xsec  
    for(int k=0; k<NMASS; ++k){
      if(fabs(sbkg[3][k][i])>0.){
	sbkg[3][k][i] = sbkg[3][k][i]/fabs(sbkg[3][k][i])*sqrt(pow(sbkg[3][k][i],2)-pow(scomm[3],2)-pow(scomm[0],2));
      }
    }
  }

  for(int i=0; i<NCHA; ++i){ // EW xsec 
    for(int k=0; k<NMASS; ++k){
      if(fabs(sbkg[4][k][i])>0.){
	sbkg[4][k][i] = sbkg[4][k][i]/fabs(sbkg[4][k][i])*sqrt(pow(sbkg[4][k][i],2)-pow(scomm[4],2)-pow(scomm[0],2));
      }
    }
  }

  TFile* f[2];
  f[0] = new TFile("shape_snntag_1fb.root","read"); 
  f[1] = new TFile("shape_secvtx_1fb.root","read");
  TFile* f3 = new TFile("shape_nunubb.root","read");   
  TFile* f4 = new TFile("shape_hww.root","read"); 
  TFile* f5 = new TFile("cdfvvbb973.root","read"); 
  TFile* f6 = new TFile("ZHllbb_1fb.root","read"); 

  TFile* f1_d0[6]; 
  f1_d0[0] = new TFile("fv_whevbb_1tag.root","read"); 
  f1_d0[1] = new TFile("fv_whevbb_2tag.root","read"); 
  f1_d0[2] = new TFile("fv_whmvbb_1tag.root","read"); 
  f1_d0[3] = new TFile("fv_whmvbb_2tag.root","read"); 
  f1_d0[4] = new TFile("fv_zhvvbb_1Tag.root","read"); 
  f1_d0[5] = new TFile("fv_zhvvbb_2Tag.root","read"); 

  TFile* f2_d0[3];   
  f2_d0[0] = new TFile("fv_hww_ee.root","read"); 
  f2_d0[1] = new TFile("fv_hww_emu.root","read"); 
  f2_d0[2] = new TFile("fv_hww_mumu_1fb.root","read"); 

  TFile* f3_d0[3];   
  f3_d0[0] = new TFile("fv_www_ee.root","read"); 
  f3_d0[1] = new TFile("fv_www_emu.root","read"); 
  f3_d0[2] = new TFile("fv_www_mumu.root","read"); 

  TFile* f4_d0[2]; 
  f4_d0[0] = new TFile("fv_zheebb_2tag.root","read"); 
  f4_d0[1] = new TFile("fv_zhmmbb_2tag.root","read"); 

  for(int i=0; i<2; ++i){ 
    XDATA[0][i]  = (TH1F*)f[i]->Get("M2j_data"); 
    XDATA[0][i] ->Rebin(4); 
    XHIGGS[0][i] = (TH1F*)f[i]->Get("M2j_WH110_CTEQ5L_ISRon_FSRon"); 
    XHIGGS[1][i] = (TH1F*)f[i]->Get("M2j_WH115_CTEQ5L_ISRon_FSRon"); 
    XHIGGS[2][i] = (TH1F*)f[i]->Get("M2j_WH120_CTEQ5L_ISRon_FSRon"); 
    XHIGGS[3][i] = (TH1F*)f[i]->Get("M2j_WH130_CTEQ5L_ISRon_FSRon"); 
    XHIGGS[4][i] = (TH1F*)f[i]->Get("M2j_WH140_CTEQ5L_ISRon_FSRon"); 
    XHIGGS[5][i] = (TH1F*)f[i]->Get("M2j_WH150_CTEQ5L_ISRon_FSRon"); 
    XBKG[0][0][i] = (TH1F*)f[i]->Get("M2j_mistag"); 
    XBKG[1][0][i] = (TH1F*)f[i]->Get("M2j_nonW"); 
    XBKG[2][0][i] = (TH1F*)f[i]->Get("M2j_HF");     
    XBKG[3][0][i] = (TH1F*)f[i]->Get("M2j_Top"); 
    XBKG[4][0][i] = (TH1F*)f[i]->Get("M2j_diboson");
    for(int j = 0; j<6; ++j){ 
      XHIGGS[j][i] ->Rebin(4); 
      if(j<5)XBKG[j][0][i] ->Rebin(4);       
    }

    for(int j =1; j<NMASS; ++j){
      XDATA[j][i] = XDATA[0][i]; 
      for(int k=0; k<NBKGH; ++k){ 
	XBKG[k][j][i] = XBKG[k][0][i];
      }
    }
  }
  // what about ZH->nunubb ?   
    XDATA[0][2]  = (TH1F*)f3->Get("Dijet_data"); 
    XHIGGS[0][2] = (TH1F*)f3->Get("Dijet_mh110"); 
    XHIGGS[2][2] = (TH1F*)f3->Get("Dijet_mh120"); 
    XHIGGS[3][2] = (TH1F*)f3->Get("Dijet_mh130"); 
    XHIGGS[1][2] =(TH1F*)XHIGGS[0][2]->Clone(); 
    XHIGGS[1][2]->Add(XHIGGS[2][2]); 
    XHIGGS[1][2]->Scale(0.5); 
    XBKG[0][0][2] = (TH1F*)f3->Get("Dijet_mistag"); 
    XBKG[1][0][2] = (TH1F*)f3->Get("Dijet_qcd"); 
    XBKG[2][0][2] = (TH1F*)f3->Get("Dijet_hf");     
    XBKG[3][0][2] = (TH1F*)f3->Get("Dijet_top"); 
    XBKG[4][0][2] = (TH1F*)f3->Get("Dijet_diboson");

    for(int j =1; j<NMASS; ++j){
      XDATA[j][2] = XDATA[0][2]; 
      for(int k=0; k<NBKGH; ++k){ 
	XBKG[k][j][2] = XBKG[k][0][2];
      }
    }
  // what about ww->llnunu 
  for(int i = 0; i<10; ++i){ 
    TString a ="data_";
    TString a1 ="hww_";
    TString a2 ="ww_";
    TString a3 ="bg_";
    int i0 = HMass[i]; 
    a +=i0; 
    a1 +=i0; 
    a2 +=i0; 
    a3 +=i0;
    XDATA[i][3] = (TH1F*) f4->Get(a); 
    XHIGGS[i][3] = (TH1F*) f4->Get(a1); 
    XBKG[4][i][3] = (TH1F*) f4->Get(a2);
    XBKG[5][i][3] = (TH1F*) f4->Get(a3);
    std::cout<<" Hww Hmass="<<i0<<" data="<<XDATA[i][3]->Integral()<<" XHIGGS="<< XHIGGS[i][3]->Integral()<<" WW="<<XBKG[4][i][3]->Integral()<<" Other="<<XBKG[5][i][3]->Integral()<<std::endl; 
  }
  // new zh->nunubb single tag analysis 
    XDATA[0][4]  = (TH1F*)f5->Get("mjj_DATA_1Tag_CENT"); 
    XDATA[0][4] -> Rebin(2); 
    XHIGGS[0][4] = (TH1F*)f5->Get("mjj_zh110_1Tag_CENT"); 
    XHIGGS[0][4] ->Add((TH1F*)f5->Get("mjj_wh110_1Tag_CENT"));
    XHIGGS[1][4] = (TH1F*)f5->Get("mjj_zh115_1Tag_CENT"); 
    XHIGGS[1][4] ->Add((TH1F*)f5->Get("mjj_wh115_1Tag_CENT"));

    XHIGGS[2][4] = (TH1F*)f5->Get("mjj_zh120_1Tag_CENT"); 
    XHIGGS[2][4] ->Add((TH1F*)f5->Get("mjj_wh120_1Tag_CENT"));
    XHIGGS[3][4] = (TH1F*)f5->Get("mjj_zh130_1Tag_CENT"); 
    XHIGGS[3][4] ->Add((TH1F*)f5->Get("mjj_wh130_1Tag_CENT"));
    for(int i =0; i<4; ++i){ 
      XHIGGS[i][4]->Rebin(2); 
    }
    XBKG[0][0][4] = (TH1F*)f5->Get("mjj_MTAG_1Tag_CENT"); // mistag 
    XBKG[1][0][4] = (TH1F*)f5->Get("mjj_QCD1_1Tag_CENT"); // QCD
    XBKG[1][0][4] ->Add((TH1F*)f5->Get("mjj_QCD2_1Tag_CENT")); 
    XBKG[2][0][4] = (TH1F*)f5->Get("mjj_Wev_1Tag_CENT"); // W+HF 
    XBKG[2][0][4] ->Add((TH1F*)f5->Get("mjj_Wmv_1Tag_CENT")); 
    XBKG[2][0][4] ->Add((TH1F*)f5->Get("mjj_Wtauv_1Tag_CENT")); 
    XBKG[3][0][4] = (TH1F*)f5->Get("mjj_ttbar_1Tag_CENT"); //top  
    XBKG[3][0][4] ->Add((TH1F*)f5->Get("mjj_tops_1Tag_CENT")); 
    XBKG[3][0][4] ->Add((TH1F*)f5->Get("mjj_topt_1Tag_CENT")); 
    XBKG[4][0][4] = (TH1F*)f5->Get("mjj_ZZ_1Tag_CENT"); //diboson  
    XBKG[4][0][4] ->Add((TH1F*)f5->Get("mjj_WZ_1Tag_CENT")); 
    XBKG[4][0][4] ->Add((TH1F*)f5->Get("mjj_WW_1Tag_CENT")); 
    XBKG[5][0][4] = (TH1F*)f5->Get("mjj_Zmm_1Tag_CENT"); //Z+HF  
    XBKG[5][0][4] ->Add((TH1F*)f5->Get("mjj_Ztautau_1Tag_CENT")); 
    XBKG[5][0][4] ->Add((TH1F*)f5->Get("mjj_Zvv_1Tag_CENT")); 
    for(int i=0; i<6; ++i){ 
      XBKG[i][0][4]->Rebin(2); 
    }

    for(int j =1; j<NMASS; ++j){
      XDATA[j][4] = XDATA[0][4]; 
      for(int k=0; k<NBKGH; ++k){ 
	XBKG[k][j][4] = XBKG[k][0][4];
      }
    }

  // new zh->nunubb double tag analysis 
    XDATA[0][5]  = (TH1F*)f5->Get("mjj_DATA_2Tag_CENT");
    XDATA[0][5] ->Rebin(2); 
    XHIGGS[0][5] = (TH1F*)f5->Get("mjj_zh110_2Tag_CENT"); 
    XHIGGS[0][5] ->Add((TH1F*)f5->Get("mjj_wh110_2Tag_CENT"));
    XHIGGS[1][5] = (TH1F*)f5->Get("mjj_zh115_2Tag_CENT"); 
    XHIGGS[1][5] ->Add((TH1F*)f5->Get("mjj_wh115_2Tag_CENT"));
    XHIGGS[2][5] = (TH1F*)f5->Get("mjj_zh120_2Tag_CENT"); 
    XHIGGS[2][5] ->Add((TH1F*)f5->Get("mjj_wh120_2Tag_CENT"));
    XHIGGS[3][5] = (TH1F*)f5->Get("mjj_zh130_2Tag_CENT"); 
    XHIGGS[3][5] ->Add((TH1F*)f5->Get("mjj_wh130_2Tag_CENT"));
    for(int i =0; i<4; ++i){ 
      XHIGGS[i][5] ->Rebin(2); 
    }
    XBKG[0][0][5] = (TH1F*)f5->Get("mjj_MTAG_2Tag_CENT"); // mistag 
    XBKG[1][0][5] = (TH1F*)f5->Get("mjj_QCD1_2Tag_CENT"); // QCD
    XBKG[1][0][5] ->Add((TH1F*)f5->Get("mjj_QCD2_2Tag_CENT")); 
    XBKG[2][0][5] = (TH1F*)f5->Get("mjj_Wev_2Tag_CENT"); // W+HF 
    XBKG[2][0][5] ->Add((TH1F*)f5->Get("mjj_Wmv_2Tag_CENT")); 
    XBKG[2][0][5] ->Add((TH1F*)f5->Get("mjj_Wtauv_2Tag_CENT")); 
    XBKG[3][0][5] = (TH1F*)f5->Get("mjj_ttbar_2Tag_CENT"); //top  
    XBKG[3][0][5] ->Add((TH1F*)f5->Get("mjj_tops_2Tag_CENT")); 
    XBKG[3][0][5] ->Add((TH1F*)f5->Get("mjj_topt_2Tag_CENT")); 
    XBKG[4][0][5] = (TH1F*)f5->Get("mjj_ZZ_2Tag_CENT"); //diboson  
    XBKG[4][0][5] ->Add((TH1F*)f5->Get("mjj_WZ_2Tag_CENT")); 
    XBKG[4][0][5] ->Add((TH1F*)f5->Get("mjj_WW_2Tag_CENT")); 
    XBKG[5][0][5] = (TH1F*)f5->Get("mjj_Zmm_2Tag_CENT"); //Z+HF  
    XBKG[5][0][5] ->Add((TH1F*)f5->Get("mjj_Ztautau_2Tag_CENT")); 
    XBKG[5][0][5] ->Add((TH1F*)f5->Get("mjj_Zvv_2Tag_CENT")); 
    for(int i=0; i<6; ++i){ 
      XBKG[i][0][5] ->Rebin(2); 
    }
    for(int j =1; j<NMASS; ++j){
      XDATA[j][5] = XDATA[0][5]; 
      for(int k=0; k<NBKGH; ++k){ 
	XBKG[k][j][5] = XBKG[k][0][5];
      }
    }
    // get zh->llbb  (these are 2-d, but treat as 1-d for now) 
    XDATA[0][6]  = (TH1F*)f6->Get("DATATot");
    XHIGGS[0][6] = (TH1F*)f6->Get("ZH110Tot"); //zh110
    XHIGGS[1][6] = (TH1F*)f6->Get("ZH115Tot"); //zh110
    XHIGGS[2][6] = (TH1F*)f6->Get("ZH120Tot"); //zh120
    XHIGGS[3][6] = (TH1F*)f6->Get("ZH130Tot"); //zh130
    XHIGGS[4][6] = (TH1F*)f6->Get("ZH140Tot"); //zh140
    XHIGGS[5][6] = (TH1F*)f6->Get("ZH150Tot"); //zh150
    XBKG[0][0][6] = (TH1F*)f6->Get("ZmisTot"); // mistag 
    XBKG[1][0][6] = (TH1F*)f6->Get("FakeTot"); // QCD or fake 
    XBKG[2][0][6] = (TH1F*)f6->Get("ZbbTot"); // Z+bb 
    XBKG[3][0][6] = (TH1F*)f6->Get("TTTot"); //top  
    XBKG[4][0][6] = (TH1F*)f6->Get("ZZTot"); //diboson  
    XBKG[4][0][6] ->Add((TH1F*)f6->Get("ZWTot")); 
    XBKG[5][0][6] = (TH1F*)f6->Get("ZccTot"); //Z+cc  

    for(int j =1; j<NMASS; ++j){
      XDATA[j][6] = XDATA[0][6]; 
      for(int k=0; k<NBKGH; ++k){ 
	XBKG[k][j][6] = XBKG[k][0][6];
      }
    }
    TH2F* XZHshape[2]; // shape systematic for zh->llbb pseudo experiment for generating different shape  
    XZHshape[0] = (TH2F*)f6->Get("ZbbshapeTot"); // Z+bb
    XZHshape[1] =  (TH2F*)f6->Get("ZccshapeTot"); //Z+cc  

  // get histograms for d0 
  for(int i=0; i<2; ++i){ //whevbb
    
    XDATA[0][i+NCHA]  = (TH1F*)f1_d0[i]->Get("Data Final Variable (mH)");
    XDATA[0][i+NCHA] ->Rebin(5); 
    XHIGGS[0][i+NCHA] = (TH1F*)f1_d0[i]->Get("Signal Interpolate 105 + 5"); 
    XHIGGS[1][i+NCHA] = (TH1F*)f1_d0[i]->Get("Signal Interpolate 115 + 0"); 
    XHIGGS[2][i+NCHA] = (TH1F*)f1_d0[i]->Get("Signal Interpolate 115 + 5"); 
    XHIGGS[3][i+NCHA] = (TH1F*)f1_d0[i]->Get("Signal Interpolate 125 + 5");
    XHIGGS[4][i+NCHA] = (TH1F*)f1_d0[i]->Get("Signal Interpolate 135 + 5");
    XHIGGS[5][i+NCHA] = (TH1F*)f1_d0[i]->Get("Signal Interpolate 145 + 5");
    for(int j =0; j<6; ++j)XHIGGS[j][i+NCHA]->Rebin(5); 
    XBKG_d0[0][0][i] = (TH1F*)f1_d0[i]->Get("WZ-Bkgd Final Variable (mH)"); 
    XBKG_d0[1][0][i] = (TH1F*)f1_d0[i]->Get("Wbb-Bkgd Final Variable (mH)"); 
    XBKG_d0[2][0][i] = (TH1F*)f1_d0[i]->Get("Wjj-Bkgd Final Variable (mH)"); 
    XBKG_d0[3][0][i] = (TH1F*)f1_d0[i]->Get("Zjj-Bkgd Final Variable (mH)"); 
    XBKG_d0[4][0][i] = (TH1F*)f1_d0[i]->Get("tt-Bkgd Final Variable (mH)"); 
    XBKG_d0[5][0][i] = (TH1F*)f1_d0[i]->Get("single top-Bkgd Final Variable (mH)"); 
    XBKG_d0[6][0][i] = (TH1F*)f1_d0[i]->Get("qcd-Bkgd Final Variable (mH)"); 
    for(int j=0; j<7; j++)XBKG_d0[j][0][i]->Rebin(5);

    for(int j =1; j<Nmax[i+NCHA]; ++j){
      XDATA[j][i+NCHA] = XDATA[0][i+NCHA]; 
      for(int k=0; k<NBKGHD0; ++k){ 
	XBKG_d0[k][j][i] = XBKG_d0[k][0][i];
      }
    }
  }

  for(int i=2; i<4; ++i){ // wmvbb
    
    XDATA[0][i+NCHA]  = (TH1F*)f1_d0[i]->Get("Data Final Variable (mH)");
    XDATA[0][i+NCHA] ->Rebin(5); 
    XHIGGS[0][i+NCHA] = (TH1F*)f1_d0[i]->Get("hwlv Interpolate 105 + 5"); 
    XHIGGS[1][i+NCHA] = (TH1F*)f1_d0[i]->Get("hwlv Interpolate 115 + 0"); 
    XHIGGS[2][i+NCHA] = (TH1F*)f1_d0[i]->Get("hwlv Interpolate 115 + 5"); 
    XHIGGS[3][i+NCHA] = (TH1F*)f1_d0[i]->Get("hwlv Interpolate 125 + 5");
    XHIGGS[4][i+NCHA] = (TH1F*)f1_d0[i]->Get("hwlv Interpolate 135 + 5");
    XHIGGS[5][i+NCHA] = (TH1F*)f1_d0[i]->Get("hwlv Interpolate 145 + 5");
    for(int j =0; j<6; ++j)XHIGGS[j][i+NCHA]->Rebin(5); 
    XBKG_d0[0][0][i] = (TH1F*)f1_d0[i]->Get("WZ-Bkgd Final Variable (mH)"); 
    XBKG_d0[1][0][i] = (TH1F*)f1_d0[i]->Get("Wbb-Bkgd Final Variable (mH)"); 
    XBKG_d0[2][0][i] = (TH1F*)f1_d0[i]->Get("Wjj-Bkgd Final Variable (mH)"); 
    XBKG_d0[3][0][i] = (TH1F*)f1_d0[i]->Get("Zmm-Bkgd Final Variable (mH)");
    XBKG_d0[3][0][i] ->Add((TH1F*)f1_d0[i]->Get("Wtv-Bkgd Final Variable (mH)"));
    XBKG_d0[4][0][i] = (TH1F*)f1_d0[i]->Get("tt-Bkgd Final Variable (mH)"); 
    XBKG_d0[5][0][i] = (TH1F*)f1_d0[i]->Get("st-Bkgd Final Variable (mH)"); 
    XBKG_d0[6][0][i] = (TH1F*)f1_d0[i]->Get("qcd-Bkgd Final Variable (mH)"); 
    for(int j=0; j<7; j++)XBKG_d0[j][0][i]->Rebin(5);

    for(int j =1; j<Nmax[i+NCHA]; ++j){
      XDATA[j][i+NCHA] = XDATA[0][i+NCHA]; 
      for(int k=0; k<NBKGHD0; ++k){ 
	XBKG_d0[k][j][i] = XBKG_d0[k][0][i];
      }
    }
  }


  for(int i=4; i<6; ++i){ // zhvvbb
    XDATA[0][i+NCHA]  = (TH1F*)f1_d0[i]->Get("Data Final Variable (mH)");
    XDATA[0][i+NCHA] ->Rebin(5); 
    XHIGGS[0][i+NCHA] = (TH1F*)f1_d0[i]->Get("Signal Interpolate 105 + 5"); 
    XHIGGS[1][i+NCHA] = (TH1F*)f1_d0[i]->Get("Signal Interpolate 115 + 0");
    XHIGGS[2][i+NCHA] = (TH1F*)f1_d0[i]->Get("Signal Interpolate 115 + 5");
    XHIGGS[3][i+NCHA] = (TH1F*)f1_d0[i]->Get("Signal Interpolate 125 + 5");
    XHIGGS[4][i+NCHA] = (TH1F*)f1_d0[i]->Get("Signal Interpolate 135 + 5");
    for(int j =0; j<5; ++j)XHIGGS[j][i+NCHA]->Rebin(5); 

    XBKG_d0[0][0][i] = (TH1F*)f1_d0[i]->Get("WZ-Bkgd Final Variable (mH)"); 
    XBKG_d0[0][0][i] ->Add((TH1F*)f1_d0[i]->Get("ZZ-Bkgd Final Variable (mH)"));
    XBKG_d0[1][0][i] = (TH1F*)f1_d0[i]->Get("Wbb-Bkgd Final Variable (mH)");
    XBKG_d0[1][0][i] ->Add((TH1F*)f1_d0[i]->Get("Zbb-Bkgd Final Variable (mH)"));
    XBKG_d0[2][0][i] = (TH1F*)f1_d0[i]->Get("Wjj-Bkgd Final Variable (mH)"); 
    XBKG_d0[3][0][i] = (TH1F*)f1_d0[i]->Get("Zjj-Bkgd Final Variable (mH)"); 
    XBKG_d0[4][0][i] = (TH1F*)f1_d0[i]->Get("tt-Bkgd Final Variable (mH)"); 
    XBKG_d0[5][0][i] = (TH1F*)f1_d0[i]->Get("qtb-Bkgd Final Variable (mH)"); 
    XBKG_d0[6][0][i] = (TH1F*)f1_d0[i]->Get("qcd-Bkgd Final Variable (mH)"); 
    for(int j=0; j<7; j++)XBKG_d0[j][0][i]->Rebin(5);

    for(int j =1; j<Nmax[i+NCHA]; ++j){
      XDATA[j][i+NCHA] = XDATA[0][i+NCHA]; 
      for(int k=0; k<NBKGHD0; ++k){ 
	XBKG_d0[k][j][i] = XBKG_d0[k][0][i];
      }
    }
  }

  // what about ww->llnunu from d0
  TString tit[NMASS]={" ", " ","120 + 0","120 + 10","140 + 0", "140 + 10","160 + 0","160 + 10","180 + 0","180 + 20"}; 
  TString tits[NMASS]={" ", " ","120","130","140","150","160","170","180","200"}; 
  for(int j=2; j<NMASS; ++j){  //hww_ee
    XDATA[j][NCHA+6] = (TH1F*) f2_d0[0]->Get("Data Final Variable (dphi) - "+tits[j]); 
    if(HMass[j] !=200){
    XHIGGS[j][NCHA+6] = (TH1F*) f2_d0[0]->Get("sig Interpolate "+tit[j]);
    XBKG_d0[0][j][6] = (TH1F*) f2_d0[0]->Get("wz Interpolate "+tit[j]);
    XBKG_d0[1][j][6] = (TH1F*) f2_d0[0]->Get("zz Interpolate "+tit[j]);
    XBKG_d0[2][j][6] = (TH1F*) f2_d0[0]->Get("ww Interpolate "+tit[j]);
    XBKG_d0[3][j][6] = (TH1F*) f2_d0[0]->Get("tt Interpolate "+tit[j]);
    XBKG_d0[4][j][6] = (TH1F*) f2_d0[0]->Get("wjet Interpolate "+tit[j]);
    XBKG_d0[5][j][6] = (TH1F*) f2_d0[0]->Get("ztt Interpolate "+tit[j]);
    XBKG_d0[6][j][6] = (TH1F*) f2_d0[0]->Get("zee Interpolate "+tit[j]);
    XBKG_d0[7][j][6] = (TH1F*) f2_d0[0]->Get("qcd Interpolate "+tit[j]);
    }
    else{ 
    XHIGGS[j][NCHA+6] = (TH1F*) f2_d0[0]->Get("Signal Final Variable (200 mH)"); 
    XBKG_d0[0][j][6] = (TH1F*) f2_d0[0]->Get("WZ-Bkgd Final Variable (dphi) - 200"); 
    XBKG_d0[1][j][6] = (TH1F*) f2_d0[0]->Get("ZZ-Bkgd Final Variable (dphi) - 200"); 
    XBKG_d0[2][j][6] = (TH1F*) f2_d0[0]->Get("WW-Bkgd Final Variable (dphi) - 200"); 
    XBKG_d0[3][j][6] = (TH1F*) f2_d0[0]->Get("tt-Bkgd Final Variable (dphi) - 200");
    XBKG_d0[4][j][6] = (TH1F*) f2_d0[0]->Get("Wjet-Bkgd Final Variable (dphi) - 200"); 
    XBKG_d0[5][j][6] = (TH1F*) f2_d0[0]->Get("ztt-Bkgd Final Variable (dphi) - 200"); 
    XBKG_d0[6][j][6] = (TH1F*) f2_d0[0]->Get("zee-Bkgd Final Variable (dphi) - 200"); 
    XBKG_d0[7][j][6] = (TH1F*) f2_d0[0]->Get("qcd-Bkgd Final Variable (dphi) - 200"); 
    }
  }

  for(int j=2; j<NMASS; ++j){  //hww_emu
    XDATA[j][NCHA+7] = (TH1F*) f2_d0[1]->Get("Data Final Variable (dphi) - "+tits[j]); 
    if(HMass[j] !=200){
    XHIGGS[j][NCHA+7] = (TH1F*) f2_d0[1]->Get("sig Interpolate "+tit[j]);
    XBKG_d0[0][j][7] = (TH1F*) f2_d0[1]->Get("wz Interpolate "+tit[j]);
    XBKG_d0[1][j][7] = (TH1F*) f2_d0[1]->Get("zz Interpolate "+tit[j]);
    XBKG_d0[2][j][7] = (TH1F*) f2_d0[1]->Get("ww Interpolate "+tit[j]);
    XBKG_d0[3][j][7] = (TH1F*) f2_d0[1]->Get("tt Interpolate "+tit[j]);
    XBKG_d0[4][j][7] = (TH1F*) f2_d0[1]->Get("wjet Interpolate "+tit[j]);
    XBKG_d0[5][j][7] = (TH1F*) f2_d0[1]->Get("ztt Interpolate "+tit[j]);
    XBKG_d0[6][j][7] = (TH1F*) f2_d0[1]->Get("zmm Interpolate "+tit[j]);
    XBKG_d0[7][j][7] = (TH1F*) f2_d0[1]->Get("qcd Interpolate "+tit[j]);
    }
    else{
    XHIGGS[j][NCHA+7] = (TH1F*) f2_d0[1]->Get("Signal Final Variable (200 mH)");
    XBKG_d0[0][j][7] = (TH1F*) f2_d0[1]->Get("WZ-Bkgd Final Variable (dphi) - 200"); 
    XBKG_d0[1][j][7] = (TH1F*) f2_d0[1]->Get("ZZ-Bkgd Final Variable (dphi) - 200"); 
    XBKG_d0[2][j][7] = (TH1F*) f2_d0[1]->Get("WW-Bkgd Final Variable (dphi) - 200"); 
    XBKG_d0[3][j][7] = (TH1F*) f2_d0[1]->Get("tt-Bkgd Final Variable (dphi) - 200"); 
    XBKG_d0[4][j][7] = (TH1F*) f2_d0[1]->Get("Wjet-Bkgd Final Variable (dphi) - 200"); 
    XBKG_d0[5][j][7] = (TH1F*) f2_d0[1]->Get("ztt-Bkgd Final Variable (dphi) - 200"); 
    XBKG_d0[6][j][7] = (TH1F*) f2_d0[1]->Get("zmm-Bkgd Final Variable (dphi) - 200"); 
    XBKG_d0[7][j][7] = (TH1F*) f2_d0[1]->Get("qcd-Bkgd Final Variable (dphi) - 200"); 
    }      
  }

  for(int j=2; j<NMASS; ++j){  //hww_mumu
    XDATA[j][NCHA+8] = (TH1F*) f2_d0[2]->Get("Data Final Variable (dphi)"); 
    XHIGGS[j][NCHA+8] = (TH1F*) f2_d0[2]->Get("Signal Interpolate "+tit[j]);
    if(j==(NMASS-1))XHIGGS[j][NCHA+8] = (TH1F*) f2_d0[2]->Get("Signal Final Variable (200 dphi)");
    XBKG_d0[0][j][8] = (TH1F*) f2_d0[2]->Get("WZ-Bkgd Final Variable (dphi)"); 
    XBKG_d0[1][j][8] = (TH1F*) f2_d0[2]->Get("Zbb-Bkgd Final Variable (dphi)"); 
    XBKG_d0[2][j][8] = (TH1F*) f2_d0[2]->Get("Zjj-Bkgd Final Variable (dphi)"); 
    XBKG_d0[3][j][8] = (TH1F*) f2_d0[2]->Get("Wtn-Bkgd Final Variable (dphi)"); 
    XBKG_d0[4][j][8] = (TH1F*) f2_d0[2]->Get("tt-Bkgd Final Variable (dphi)"); 
    XBKG_d0[5][j][8] = (TH1F*) f2_d0[2]->Get("qcd-Bkgd Final Variable (dphi)"); 
  }
  //
  // what about wh->www->llnunu from d0
  TString titx[NMASS]={"100 + 10", "115 + 0", "115 + 5","115 + 15","135 + 5", "135 + 15","155 + 5","155 + 15","175 + 5","175 + 25"}; 
  TString titxs[NMASS]={"110", "115","120","130","140","150","160","170","180","200"}; 
  for(int i =0; i<3; ++i){ // www_ee, www_emu, www_mumu 
  for(int j=0; j<NMASS; ++j){ 
    XDATA[j][NCHA+9+i] = (TH1F*) f3_d0[i]->Get("Data Final Variable (dphi) - "+titxs[j]); 
    XHIGGS[j][NCHA+9+i] = (TH1F*) f3_d0[i]->Get("sig Interpolate "+titx[j]);
    XHIGGS[j][NCHA+9+i] -> Rebin(4); 
    XBKG_d0[0][j][9+i] = (TH1F*) f3_d0[i]->Get("WZ Interpolate "+titx[j]);
    XBKG_d0[1][j][9+i] = (TH1F*) f3_d0[i]->Get("ZZ Interpolate "+titx[j]);
    XBKG_d0[2][j][9+i] = (TH1F*) f3_d0[i]->Get("flip Interpolate "+titx[j]);
    XBKG_d0[3][j][9+i] = (TH1F*) f3_d0[i]->Get("qcd Interpolate "+titx[j]);
  }
  }

  // what about zh->llbb from d0
  TString tity[NMASS]={"105 + 5", "115 + 0","115 + 5","125 + 5","135 + 5", "145 + 5"," "," "," "," "}; 
  for(int i =0; i<2; ++i){ // eebb, mmbb
    XDATA[0][NCHA+12+i] = (TH1F*) f4_d0[i]->Get("Data Final Variable (mH)");
    XDATA[0][NCHA+12+i] ->Rebin(2); 
    for(int j = 0;  j<Nmax[NCHA+12+i]; ++j){ 
    XHIGGS[j][NCHA+12+i] = (TH1F*) f4_d0[i]->Get("hwlv Interpolate "+tity[j]);
    XHIGGS[j][NCHA+12+i] ->Rebin(2);
    }
    XBKG_d0[0][0][12+i] = (TH1F*) f4_d0[i]->Get("WZ-Bkgd Final Variable (mH)");
    XBKG_d0[1][0][12+i] = (TH1F*) f4_d0[i]->Get("ZZ-Bkgd Final Variable (mH)");
    XBKG_d0[2][0][12+i] = (TH1F*) f4_d0[i]->Get("Zbb-Bkgd Final Variable (mH)");
    XBKG_d0[3][0][12+i] = (TH1F*) f4_d0[i]->Get("Zjj-Bkgd Final Variable (mH)");
    XBKG_d0[4][0][12+i]  = (TH1F*) f4_d0[i]->Get("tt-Bkgd Final Variable (mH)");
    XBKG_d0[5][0][12+i]  = (TH1F*) f4_d0[i]->Get("qcd-Bkgd Final Variable (mH)");
    for(int k=0; k<6; k++)XBKG_d0[k][0][12+i]->Rebin(2);
 
    for(int j =1; j<Nmax[i+NCHA+12]; ++j){
      XDATA[j][i+NCHA+12] = XDATA[0][i+NCHA+12]; 
      for(int k=0; k<NBKGHD0-2; ++k){ 
	XBKG_d0[k][j][12+i] = XBKG_d0[k][0][12+i];
      }
    }
  }

  //
  const Int_t nXsecR = 300; // 300
  static Int_t nIter = 20000; // 100k or 500k 
  if(NoSYS==3) nIter =nXsecR; 
  const float step =1.; 
  // random number
  TRandom gRand;
  gRand.SetSeed();

  // cross section step size
  Int_t NMASSX(0); 
  Int_t nchannel(0); 
  for(Int_t i=0; i<NCHAT; ++i){ 
    if(icha[i]>0){
      nchannel++; 
      if(NMASSX<Nmax[i])NMASSX = Nmax[i]; 
    }
  }
   
  // make copy of data: (starting 120 GeV) 
  for(Int_t i =0; i<NCHAT; ++i){ 
    if(XDATA[2][i]){ 
      if(XDATA[2][i]->GetDimension()<2){ 
	XMC[i] = (TH1F*)XDATA[2][i]->Clone(); 
      }
      else{ 
	TH2F* x2d = (TH2F*)XDATA[2][i];
	XMC[i] = (TH1F*)x2d->Clone(); 
      }
    }
  }
  // combined chi^{2} and likelihood histograms

  TH1F *hchi[NMASS];
  TH1F *hlh [NMASS];
  TH1F *hl[NMASS]; 
  TH1F *hevt[NCHAT][NMASS]; 
  for ( Int_t i = 0 ; i < NMASSX ; i++ ) {
      TString nchi("hchi_");
      int j = HMass[i]; 
      nchi += (Int_t)j;

      TString nlh("hlh_");
      nlh += (Int_t)j;

      TString nl("hlimit_");
      nl += (Int_t)j;

      hchi[i] = new TH1F(nchi,nchi,nXsecR,0.,step*nXsecR);
      hchi[i]->SetXTitle("#sigma Ratio");
      hchi[i]->SetYTitle("#chi^{2}");

      hlh[i]  = new TH1F(nlh,nlh,nXsecR,0.,step*nXsecR);
      hlh[i]->SetXTitle("#sigma Ratio");
      hlh[i]->SetYTitle("likelihood");

      hl[i]  = new TH1F(nl,nl,nXsecR,0.,step*nXsecR);
      hl[i]->SetXTitle("#sigma Ratio");
      hl[i]->SetYTitle("Pseudoexperiments");
      for(Int_t ix=0; ix<NCHAT; ++ix){ 
	TString ne("nevt_"); 
	ne +=(Int_t)j;
	ne +="_"; 
	ne +=txt[ix];
	hevt[ix][i] = new TH1F(ne,ne, 500,0., 500.);
      }
  }

  std::cout  <<" Channels Combined="<<nchannel<< " " <<std::endl; 
  if(NoSYS==1||NoSYS==2||NoSYS==3) std::cout<<" With no systematic included "<<" "<<NoSYS<<std::endl; 
  // 95% C.L. upper limit

  for(Int_t ie =0; ie<Npseudo; ++ie){ // starting the pseudo-experiments here
   
  Double_t fit[NMASS][3] = {0.}; // central value, +sigma, -sigma
  Double_t limit[NMASS] = {0.};

  // start mass loop
  for ( Int_t im = 0 ; im < NMASSX ; im++ ) {
    if(NoSYS>100&&HMass[im]!=NoSYS)continue;
    hchi[im]->Reset(); 
    hlh[im]->Reset(); 
    // make data: 
      float gcomm[NCOMM]; // common systematic between cdf and d0
      float gcommb_lum; // try to uncorrelated the signal and background for luminosity 
      float geff_d0[NACCD0][NCHAD0]; // d0 acceptances of systematic 
      float gbkg_d0[NBKGD0][NCHAD0]; // d0 backgrounds of systematic 
      float geff[NACC][NCHA]; //cdf acceptance of systematic 
      float geffb[NACC][NCHA]; // possible to break down the signal and backgrounds 
      float gbkg[NBKG][NCHA]; // cdf backgrounds of systematic 
      float gjes_nnbb[2][NBKGH]; // for nunubb jes on the backgrounds  
      float gtagc; // tagging c probability at cdf 
      // modeling efficiencies 
      for(int i = 0; i<NCOMM; ++i){ 
	gcomm[i] = getGauss(gRand, 1., fabs(scomm[i]) );
      }
      if(NoSYS==4){ 
	gcommb_lum = getGauss(gRand, 1., fabs(scomm[0]) ); 
      }
      else{
	gcommb_lum = gcomm[0]; 
      }

      for(int i =0; i<NACCD0; ++i){  // for D0 treat the signal and background correlated 
	float val; 
	do{
	  float a = gRand.Gaus(0.,1.0);
	  float ab = a; 
	  if(NoSYS==4)ab = gRand.Gaus(0.,1.0);
	  val =1.0; 
	  for(int j =0; j<NCHAD0; ++j){ 
	      geff_d0[i][j] =(1.0+a*seff_d0[i][j]);
	      gbkg_d0[i][j] =(1.0+ab*sbkg_d0[i][im][j]); 
	      if(geff_d0[i][j]<=0.)val = geff_d0[i][j];
	      if(gbkg_d0[i][j]<=0.)val = gbkg_d0[i][j]; 
	  }
	}while(val<=0.); 
      }
	
      for(int i =NACCD0; i<NBKGD0; ++i){ 
	float val; 
	do{
	  float a = gRand.Gaus(0.,1.0);
	  val =1.0; 
	  for(int j =0; j<NCHAD0; ++j){ 
	      gbkg_d0[i][j] =(1.0+a*sbkg_d0[i][im][j]);
	      if(gbkg_d0[i][j]<=0.)val = gbkg_d0[i][j]; 
	  }
	}while(val<=0.); 
      }

      // modeling efficiency >0 for CDF 
      for(int i =0; i<NACC; ++i){  // for CDF treat the signal correlated 
	float val; 
	do{
	  float a = gRand.Gaus(0.,1.0);
	  float ab = a; 
	  if(NoSYS==4)ab =gRand.Gaus(0.,1.0); 
	  val =1.0; 
	  for(int j =0; j<NCHA; ++j){ 
	      if(seff[i][j]>0){
		geff[i][j] =(1.0+a*seff[i][j]);
		geffb[i][j] =(1.0+ab*seff[i][j]);
	      }
	      else{
		geff[i][j] = getGauss(gRand, 1., fabs(seff[i][j]));
		geffb[i][j] = geff[i][j]; 
		if(NoSYS==4)geffb[i][j] = getGauss(gRand, 1., fabs(seff[i][j]));
	      }
	      if(geff[i][j]<=0.)val = geff[i][j];
	      if(geffb[i][j]<=0.)val = geffb[i][j];
	  }
	  if(i==1){ 
	    gtagc = 1.0+0.16*ab; 
	    if(gtagc<=0.)val = gtagc; 
	  }
	  if(i==3){ 
	    for(int j=0; j<5; ++j){ // jes uncertainties for nunubb
	      if(sjes_nnbb[0][j]>0){
		gjes_nnbb[0][j] = (1.0 +a*sjes_nnbb[0][j]); 
	      }
	      else{ 
		gjes_nnbb[0][j] = getGauss(gRand, 1.0, fabs(sjes_nnbb[0][j])); 
	      }
              if(sjes_nnbb[1][j]>0){
		gjes_nnbb[1][j] = (1.0 +a*sjes_nnbb[1][j]);
	      }
	      else{
		gjes_nnbb[1][j] = getGauss(gRand, 1.0, fabs(sjes_nnbb[1][j])); 
	      } 
	      //gjes_nnbb[0][j] = (1.0 +a*sjes_nnbb[0][j]);
              //gjes_nnbb[1][j] = (1.0 +a*sjes_nnbb[1][j]); 
	      if(gjes_nnbb[0][j]<=0.)val = gjes_nnbb[0][j];
	      if(gjes_nnbb[1][j]<=0.)val = gjes_nnbb[1][j];
	    }
	  }
	}while(val<=0.); 
      }
	
      // modeling background >0 
      // for CDF treat the background systematic and including some of systematic from acceptance too
      // For example some of uncertainties are partial only in nunubb and zh->llbb 
      for(int i =0; i<NBKG; ++i){  
	float val; 
	do{
	  float a = gRand.Gaus(0.,1.0);
	  val =1.0; 
	  for(int j =0; j<NCHA; ++j){ 
	      if(sbkg[i][im][j]>0){
		gbkg[i][j] =(1.0+a*sbkg[i][im][j]);
	      }
	      else{
		gbkg[i][j] = getGauss(gRand, 1., fabs(sbkg[i][im][j]));
	      }
	      if(gbkg[i][j]<=0.)val = gbkg[i][j];
	  }
	}while(val<=0.); 
      }

    //
      for( Int_t ich=0; ich<NCHAT; ++ich){ 
	if(icha[ich]>0&&im<Nmax[ich]&&(ich<(NCHAT-8)||(ich>=(NCHAT-8)&&ich<(NCHAT-5)&&im>1)||ich>=(NCHAT-5)) ){// there is no 110 for D0 ww analysis 
 
	  XMC[ich]->Reset();
	  if(ich<NCHA-1){ // cdf channels 
	  // eff 
	  float xb[NBKGH]={0.}; // background normalization
	  if(ich<4){ //lvbb_s, lvbb_d, nnbb, ww 
	  for( Int_t j=0; j<NBKGH-1; ++j){ 	   
	    if(fabs(sbkg[j][im][ich])>0.)xb[j] =gbkg[j][ich]; 
	    if(j==3) xb[j] *=gcomm[3]*gcommb_lum*geffb[0][ich]; //common top xsec 
	    if(j==4) xb[j] *=gcomm[4]*gcommb_lum*geffb[0][ich]; // common diboson 
	  }
	  if(fabs(sbkg[NBKG-1][im][ich])>0.)xb[NBKGH-1] = gbkg[NBKG-1][ich]; // special for ww 
	  }
	  else if(ich<6){ // nunubb_s, nunubb_d

	  for( Int_t j=0; j<NBKGH; ++j){ 	   
	    if(j<2&&fabs(sbkg[j][im][ich])>0.){
	      xb[j] =gbkg[j][ich]; 
	    }
	    else if(j<(NBKGH-1)&&(fabs(sbkg[j][im][ich])>0.||fabs(sbkg[j+3][im][ich])>0.)){
	      xb[j]=gbkg[j][ich]*gbkg[j+3][ich]; //w+hf
	    }
	    else if(j==(NBKGH-1)&&(fabs(sbkg[2][im][ich])>0.||fabs(sbkg[j+3][im][ich])>0.)){
	      xb[j] = gbkg[2][ich]*gbkg[j+3][ich]; //z+hf
	    }
	    if(j==3) xb[j] *=gcomm[3]; //common top xsec 
	    if(j==4) xb[j] *=gcomm[4]; // common diboson
	    float aj = gjes_nnbb[0][j]; 
	    if(ich==5)aj = gjes_nnbb[1][j];
	    if(j>1)xb[j] *=aj*gcommb_lum*geffb[0][ich]*geffb[5][ich]*geffb[1][ich]*geffb[2][ich]*geffb[4][ich]; //no addition sys for mistag/qcd
	  }

	  }
	  for(Int_t k =1; k<XDATA[im][ich]->GetNbinsX()+1; ++k){ 

	    Double_t mc = 0.; 
	    //background 
	    for( Int_t j=0; j<NBKGH; ++j){ 
	      if(xb[j]>0.)mc +=xb[j]*XBKG[j][im][ich]->GetBinContent(k); // xb[j]==0, the hist is not defined
	    }
	    Int_t nx = gRand.Poisson(mc); 
	    XMC[ich]->SetBinContent(k,nx); 
	    XMC[ich]->SetBinError(k,sqrt(nx+0.001)); 
	  } //loop over the histogram bins 
	  hevt[ich][im]->Fill(XMC[ich]->Integral());
	  } // for cdf channels
	  else if(ich==(NCHA-1)){
            //*************************************************** 
	    // for special treat of 2-d llbb fit from cdf 
	    //***************************************************
	    float xb[NBKGH]={0.}; // background normalization
	    for( Int_t j=0; j<NBKGH; ++j){ 	   
	      
	    if(j<(NBKGH-1)&&fabs(sbkg[j][im][ich])>0.){
	      xb[j] =gbkg[j][ich];
	      if(j>1)xb[j] *=geffb[1][ich]; // btag
	    }
	    else if(j==(NBKGH-1)&&fabs(sbkg[2][im][ich])>0.){ 
	      xb[j] =gbkg[2][ich]*gtagc; //ctag
	    }
	    if(j==3) xb[j] *=gcomm[3]; //common top xsec 
	    if(j==4) xb[j] *=gcomm[4]; // common diboson 
	    if(j>1)xb[j] *=gcommb_lum*geffb[0][ich]*geffb[2][ich]*geffb[3][ich]*geffb[4][ich]; // add additional systematic:lum,lep id, jes, isr/fsr+pdf   
	    }
	    // convert 1-d to 2d; 
	    TH2F* xd2d =(TH2F*)XMC[ich]; 
	    TH2F* xb2d[NBKGH]; 
	    for(Int_t i=0; i<NBKGH; ++i) xb2d[i] =(TH2F*)XBKG[i][im][ich]; 
	    if(NoSYS==10){
	      xb2d[2] =XZHshape[0]; 
	      xb2d[5] =XZHshape[1]; 
	    }

	    for(Int_t k =1; k<xd2d->GetNbinsX()+1; ++k){ 
	      for(Int_t k1=1; k1<xd2d->GetNbinsY()+1; ++k1){  
		
		Double_t mc = 0.; 
		//background 
		for( Int_t j=0; j<NBKGH; ++j){ 
		  if(xb[j]>0.)mc +=xb[j]*xb2d[j]->GetBinContent(k,k1); // xb[j]==0, the hist is not defined
		}
		Int_t nx =gRand.Poisson(mc); 
		xd2d->SetBinContent(k,k1,nx); 
		xd2d->SetBinError(k,k1,sqrt(nx+0.001)); 
	      } //loop over the histogram bins y 
	    } // loop over the histogram bins x 
	    XMC[ich] = (TH1F*)xd2d; 
	    hevt[ich][im]->Fill(xd2d->Integral());
	  }
	  //*****************************************************************************
          // liklihood for D0 channels 
          //*****************************************************************************
	  else{ // for d0 channels 
	  // eff 

	  Double_t acc_bkg(1.0); 
	  for( Int_t i=0; i<NACCD0; ++i) {
	    acc_bkg *=gbkg_d0[i][ich-NCHA];
	      }
	  acc_bkg *=gcommb_lum; // common luminosity 
	  
	  float xb[NBKGHD0]={0.}; // background normalization 
	  // signal id
	  float eids=geff_d0[8][ich-NCHA]*geff_d0[9][ich-NCHA]*geff_d0[10][ich-NCHA]; 
	  float mids=geff_d0[12][ich-NCHA]*geff_d0[13][ich-NCHA];
	  float bids = geff_d0[5][ich-NCHA]*geff_d0[6][ich-NCHA]; 
	  // background id
	  float eid=gbkg_d0[8][ich-NCHA]*gbkg_d0[9][ich-NCHA]*gbkg_d0[10][ich-NCHA]; 
	  float mid=gbkg_d0[12][ich-NCHA]*gbkg_d0[13][ich-NCHA];
          float bid =gbkg_d0[5][ich-NCHA]*gbkg_d0[6][ich-NCHA]; 	   
	  // dileptons extra id (ee,mm), double btags; ignore for now by reset to 1.  
	  eids =1.; 
	  mids =1.; 
	  bids =1.; 
	  eid = 1.; 
	  mid = 1.; 
	  bid = 1.; 
	  if(txt[ich].Contains("d0wwee")){
	    for( Int_t j=0; j<NBKGHD0; ++j){ 	   
	      if(j <3){
		xb[j] = eid*acc_bkg*gcomm[4]; //*Xsec_EM
	      }
	      else if( j==3){
		xb[j] = eid*acc_bkg*gcomm[3]; // *Xsec_top
	      }
	      else if(j<NBKGHD0-1){
		xb[j] = eid*acc_bkg*gcomm[4]; //*Xsec_EM 		
	      }
	      else{
		xb[j] = eid*acc_bkg*gbkg_d0[22][ich-NCHA]; //*Xsec_qcd 	
	      }
	    }	    
	  }
	  else if(txt[ich].Contains("d0wwem")){
	    for( Int_t j=0; j<NBKGHD0; ++j){ 	   
	      if(j <3){
		xb[j] = acc_bkg*gcomm[4]; //*Xsec_EM
	      }
	      else if( j==3){
		xb[j] = acc_bkg*gcomm[3]; // *Xsec_top
	      }
	      else if(j<NBKGHD0-1){ 
		xb[j] = acc_bkg*gcomm[4]; //*Xsec_EM 	
	      }
	      else{
		xb[j] = acc_bkg*gbkg_d0[22][ich-NCHA]; //*Xsec_qcd 
	      }
	    }	    
	  }
	  else if(txt[ich].Contains("d0wwmm")){
	    for( Int_t j=0; j<NBKGHD0-2; ++j){ 	   
	      if(j==0||j==3){
		xb[j] = mid*acc_bkg*gcomm[4]; //*Xsec_EM
	      }
	      else if(j==1){ 
		xb[j] = mid*acc_bkg*gbkg_d0[20][ich-NCHA]; //*Xsec_HF
	      }
	      else if(j==2){ 
		xb[j] = mid*acc_bkg*gbkg_d0[21][ich-NCHA]; //*Xsec_LF
	      }
              else if(j==4){ 
		xb[j] = mid*acc_bkg*gcomm[3]; // *Xsec_top
	      }
	      else{
		xb[j] = mid*acc_bkg*gbkg_d0[22][ich-NCHA]; //*Xsec_qcd 	
	      }
	    }	    
	  }
	  else if(txt[ich].Contains("d0wh")){
	    for( Int_t j=0; j<4; ++j){ 	   
	      if(j<2){
		xb[j] = acc_bkg*gcomm[4]; //*Xsec_EM
	      }
	      else{
		xb[j] = acc_bkg*gbkg_d0[22][ich-NCHA]; //*Xsec_qcd 	
	      }
	    }	    
	  }
 	  else if((txt[ich]=="d0zh")){ // d0zheebb or d0zhmmbb 2tags
 	    for( Int_t j=0; j<NBKGHD0-2; ++j){ 	   
 	      if(j <2){
 		xb[j] = acc_bkg*gcomm[4]; //eid*Lum*trigger*Xsec_EM
 	      }
 	      else if(j<3){
 		xb[j] = acc_bkg*gbkg_d0[20][ich-NCHA]; //*Xsec_HF
 	      }
 	      else if( j<4){
 		xb[j] = acc_bkg*gbkg_d0[19][ich-NCHA]; //*Xsec_bTag_LF
 	      }
 	      else if( j<5){
 		xb[j] = acc_bkg*gcomm[3]; // *Xsec_top
 	      }
 	      else{ 
		xb[j] = acc_bkg*gbkg_d0[22][ich-NCHA]; //*Xsec_qcd 	
 	      }	      	    
 	    }
	  }  
	  else if((txt[ich]=="d0evbb_s")||(txt[ich]=="d0mvbb_s")||(txt[ich]=="d0vvbb_s")){ // single tag
	    for( Int_t j=0; j<NBKGHD0-1; ++j){ 	   
	      if(j <1){
		xb[j] = acc_bkg*gcomm[4]; //eid*Lum*trigger*Xsec_EM
	      }
	      else if(j<2){
		xb[j] = acc_bkg*gbkg_d0[20][ich-NCHA]; //*Xsec_HF
	      }
	      else if( j<4){
		xb[j] = acc_bkg*gbkg_d0[19][ich-NCHA]; //*Xsec_bTag_LF
	      }
	      else if( j<6){
		xb[j] = acc_bkg*gcomm[3]; // *Xsec_top
	      }
	      else{ 
		xb[j] = acc_bkg*gbkg_d0[22][ich-NCHA]; //*Xsec_qcd 	
	      }	      	    
	    }
	  }  
	  else{ // double tags 
	    for( Int_t j=0; j<NBKGHD0-1; ++j){ 	   
	      if(j <1){
		xb[j] = bid*acc_bkg*gcomm[4]; //*Xsec_EM
	      }
	      else if( j<2){ 
		xb[j] = bid*acc_bkg*gbkg_d0[20][ich-NCHA]; //*Xsec_HF
	      }
	      else if( j<4){
		xb[j] = bid*acc_bkg*gbkg_d0[19][ich-NCHA]; //*Xsec_bTag_LF
	      }
	      else if( j<6){
		xb[j] = bid*acc_bkg*gcomm[3]; // *Xsec_top
	      }
	      else{ 
		xb[j] = bid*acc_bkg*gbkg_d0[22][ich-NCHA]; //*Xsec_qcd 	
	      }	      	    
	    }
	  }  
	  //		

	  for(Int_t k =1; k<XDATA[im][ich]->GetNbinsX()+1; ++k){ 

	    Double_t mc = 0.; 
	    //background 
	    Int_t nx = NBKGHD0-1; 
	    if(txt[ich].Contains("d0wwmm")){
	      nx = NBKGHD0-2; 
	    }
	    else if(txt[ich].Contains("d0ww")){
	      nx = NBKGHD0;
	    }
	    if(txt[ich].Contains("d0wh")) nx = NBKGHD0-4; 
 	    if(txt[ich].Contains("d0zh"))nx = NBKGHD0-2;
	    for( Int_t j=0; j<nx; ++j){
	      if(xb[j]>0.)mc +=xb[j]*XBKG_d0[j][im][ich-NCHA]->GetBinContent(k); // xb[j]==0, the hist is not defined
	    }
	    nx =gRand.Poisson(mc); 

	    XMC[ich]->SetBinContent(k, nx); 
	    XMC[ich]->SetBinError(k,sqrt(nx+0.001)); 	    
	  } //loop over the histogram bins
	  hevt[ich][im]->Fill(XMC[ich]->Integral());
	  }
	} // channels in consideration 	
      } // loop over all the channels
      //
    // start sampling loop
    Double_t lh[nXsecR] = {0.};
    Double_t nc[nXsecR] =  {0.};
    for ( Int_t Iter = 0 ; Iter < nIter ; Iter++ ) {
      // get fluctuation: for eff 
      //*****************************************************
       float gcomm[NCOMM]; // common systematic between cdf and d0 
       float gcommb_lum; // lum for backgrounds
      float geff_d0[NACCD0][NCHAD0]; // d0 acceptances of systematic 
      float gbkg_d0[NBKGD0][NCHAD0]; // d0 backgrounds of systematic 
      float geff[NACC][NCHA]; //cdf acceptance of systematic 
      float geffb[NACC][NCHA]; // possible to break down the signal and backgrounds 
      float gbkg[NBKG][NCHA]; // cdf backgrounds of systematic 
      float gjes_nnbb[2][NBKGH]; // for nunubb jes on the backgrounds  
      float gtagc; // tagging c probability at cdf 
      // modeling efficiencies 
      for(int i = 0; i<NCOMM; ++i){ 
	gcomm[i] = getGauss(gRand, 1., fabs(scomm[i]) );
      }
      if(NoSYS==4){
	gcommb_lum = getGauss(gRand, 1., fabs(scomm[0]) );
      }
      else {
	gcommb_lum = gcomm[0]; 
      }

      for(int i =0; i<NACCD0; ++i){  // for D0 treat the signal and background correlated 
	float val; 
	do{
	  float a = gRand.Gaus(0.,1.0);
	  float ab = a; 
	  if(NoSYS==4)ab = gRand.Gaus(0.,1.0);
	  val =1.0; 
	  for(int j =0; j<NCHAD0; ++j){ 
	      geff_d0[i][j] =(1.0+a*seff_d0[i][j]);
	      gbkg_d0[i][j] =(1.0+ab*sbkg_d0[i][im][j]); 
	      if(geff_d0[i][j]<=0.)val = geff_d0[i][j];
	      if(gbkg_d0[i][j]<=0.)val = gbkg_d0[i][j]; 
	  }
	}while(val<=0.); 
      }
	
      for(int i =NACCD0; i<NBKGD0; ++i){ 
	float val; 
	do{
	  float a = gRand.Gaus(0.,1.0);
	  val =1.0; 
	  for(int j =0; j<NCHAD0; ++j){ 
	      gbkg_d0[i][j] =(1.0+a*sbkg_d0[i][im][j]);
	      if(gbkg_d0[i][j]<=0.)val = gbkg_d0[i][j]; 
	  }
	}while(val<=0.); 
      }

      // modeling efficiency >0 for CDF 
      for(int i =0; i<NACC; ++i){  // for CDF treat the signal correlated 
	float val; 
	do{
	  float a = gRand.Gaus(0.,1.0);
	  float ab = a; 
	  if(NoSYS==4)ab = gRand.Gaus(0.,1.0);
	  val =1.0; 
	  for(int j =0; j<NCHA; ++j){ 
	      if(seff[i][j]>0){
		geff[i][j] =(1.0+a*seff[i][j]);
		geffb[i][j] =(1.0+ab*seff[i][j]);
	      }
	      else{
		geff[i][j] = getGauss(gRand, 1., fabs(seff[i][j]));
		geffb[i][j] = geff[i][j]; 
		if(NoSYS==4)geffb[i][j] = getGauss(gRand, 1., fabs(seff[i][j]));
	      }
	      if(geff[i][j]<=0.)val = geff[i][j];
	      if(geffb[i][j]<=0.)val = geffb[i][j];
	  }
	  if(i==1){ 
	    gtagc = 1.0+0.16*ab; 
	    if(gtagc<=0.)val = gtagc; 
	  }
	  if(i==3){ 
	    for(int j=0; j<5; ++j){ // jes uncertainties for nunubb
	      if(sjes_nnbb[0][j]>0){
		gjes_nnbb[0][j] = (1.0 +a*sjes_nnbb[0][j]); 
	      }
	      else{ 
		gjes_nnbb[0][j] = getGauss(gRand, 1.0, fabs(sjes_nnbb[0][j])); 
	      }
              if(sjes_nnbb[1][j]>0){
		gjes_nnbb[1][j] = (1.0 +a*sjes_nnbb[1][j]);
	      }
	      else{
		gjes_nnbb[1][j] = getGauss(gRand, 1.0, fabs(sjes_nnbb[1][j])); 
	      } 
	      //gjes_nnbb[0][j] = (1.0 +a*sjes_nnbb[0][j]); 
	      //gjes_nnbb[1][j] = (1.0 +a*sjes_nnbb[1][j]); 
	      if(gjes_nnbb[0][j]<=0.)val = gjes_nnbb[0][j];
	      if(gjes_nnbb[1][j]<=0.)val = gjes_nnbb[1][j]; 
	    }
	  }
	}while(val<=0.); 
      }
	
      // modeling background >0 
      // for CDF treat the background systematic and including some of systematic from acceptance too
      // For example some of uncertainties are partial only in nunubb and zh->llbb 
      for(int i =0; i<NBKG; ++i){  
	float val; 
	do{
	  float a = gRand.Gaus(0.,1.0);
	  val =1.0; 
	  for(int j =0; j<NCHA; ++j){ 
	      if(sbkg[i][im][j]>0){
		gbkg[i][j] =(1.0+a*sbkg[i][im][j]);
	      }
	      else{
		gbkg[i][j] = getGauss(gRand, 1., fabs(sbkg[i][im][j]));
	      }
	      if(gbkg[i][j]<=0.)val = gbkg[i][j];
	  }
	}while(val<=0.); 
      }

      // additional smear of ww xsec 
      float xww = getGauss(gRand,1.0, 0.1); //assign additional 10% 
      if(NoSYS==3||NoSYS==1) xww =1.; 

      // shape sampling  between two with prob = 0,1
      float xsampling = gRand.Uniform(); 
      // get cross section Ratio from 0 to nXsecR*step  
      Double_t tXsec = gRand.Uniform()*nXsecR*step; 
      Int_t tBin = (Int_t)(tXsec/step);
      if(NoSYS==3) {
	tXsec =(Iter+0.5)*step;
	tBin = Iter; 
      }
      nc[tBin]++;

      // likelihood calculation
      Double_t loglhS = 0.;
      Double_t efftot = 0.; 
      for( Int_t ich=0; ich<NCHAT; ++ich){ 
	if(icha[ich]>0&&im<Nmax[ich]&&(ich<(NCHAT-8)||(ich>=(NCHAT-8)&&ich<(NCHAT-5)&&im>1)||ich>=(NCHAT-5))){ // there is no 110 for D0 ww analysis
	  if(ich<NCHA-1){ // cdf channels 
	  // eff 
	  Double_t acc(1.); 
	  for( Int_t i=0; i<NACC; ++i) {
	    acc *=geff[i][ich]; 
	  }
	  acc *=gcomm[0]; // commom tev luminosity 
	  acc *=sfnnlo[ich]; // correct the nnlo factor for Higgs normalization 
	  if(txt[ich]=="ww")acc *=xww; // assign additional 10% uncertainties 
	  float xb[NBKGH]={0.}; // background normalization
	  if(ich<4){ //lvbb_s, lvbb_d, nnbb, ww 
	  for( Int_t j=0; j<NBKGH-1; ++j){ 	   
	    if(fabs(sbkg[j][im][ich])>0.)xb[j] =gbkg[j][ich]; 
	    if(j==3) xb[j] *=gcomm[3]*gcommb_lum*geffb[0][ich]; //common top xsec 
	    if(j==4) xb[j] *=gcomm[4]*gcommb_lum*geffb[0][ich]; // common diboson 
	  }
	  if(fabs(sbkg[NBKG-1][im][ich])>0.)xb[NBKGH-1] = gbkg[NBKG-1][ich]; // special for ww 
	  }
	  else if(ich<6){ // nunubb_s, nunubb_d

	  for( Int_t j=0; j<NBKGH; ++j){ 	   
	    if(j<2&&fabs(sbkg[j][im][ich])>0.){
	      xb[j] =gbkg[j][ich]; 
	    }
	    else if(j<(NBKGH-1)&&(fabs(sbkg[j][im][ich])>0.||fabs(sbkg[j+3][im][ich])>0.)){
	      xb[j]=gbkg[j][ich]*gbkg[j+3][ich]; //w+hf
	    }
	    else if(j==(NBKGH-1)&&(fabs(sbkg[2][im][ich])>0.||fabs(sbkg[j+3][im][ich])>0.)){
	      xb[j] = gbkg[2][ich]*gbkg[j+3][ich]; //z+hf
	    }
	    if(j==3) xb[j] *=gcomm[3]; //common top xsec 
	    if(j==4) xb[j] *=gcomm[4]; // common diboson 
	    float aj = gjes_nnbb[0][j];  // single tag
	    if(ich==5)aj = gjes_nnbb[1][j]; //double tag 
	    if(j>1)xb[j] *=aj*gcommb_lum*geffb[0][ich]*geffb[5][ich]*geffb[1][ich]*geffb[2][ich]*geffb[4][ich]; //no addition sys for mistag/qcd
	  }

	  }
	  // No SYS applied here 
	    if(NoSYS==3||NoSYS==1)acc = sfnnlo[ich]; 
	    if(NoSYS==3||NoSYS==2){
	    for(Int_t j =0; j<NBKGH; ++j){
	      if(xb[j]>0)xb[j] = 1.;
	    }
	    }
	  for(Int_t k =1; k<XDATA[im][ich]->GetNbinsX()+1; ++k){ 

	    Double_t mc = 0.; 
	    mc += tXsec*acc*XHIGGS[im][ich]->GetBinContent(k); 
	    efftot +=acc*XHIGGS[im][ich]->GetBinContent(k);

	    //background 
	    for( Int_t j=0; j<NBKGH; ++j){ 
	      if(xb[j]>0.)mc +=xb[j]*XBKG[j][im][ich]->GetBinContent(k); // xb[j]==0, the hist is not defined
	    }
	    //Int_t nx = (Int_t) XDATA[im][ich]->GetBinContent(k); 
	    Int_t nx = (Int_t) XMC[ich]->GetBinContent(k); 
	    if(mc>0) loglhS +=nx*TMath::Log(mc)-mc-lgm[nx]; 
	  } //loop over the histogram bins 
	  } // for cdf channels
	  else if(ich==(NCHA-1)){
            //*************************************************** 
	    // for special treat of 2-d llbb fit from cdf 
	    //***************************************************
	    // eff 
	    Double_t acc(1.); 
	    for( Int_t i=0; i<NACC; ++i) {
	      acc *=geff[i][ich]; 
	    }
	    acc *=gcomm[0]; // commom tev luminosity 
	    acc *=sfnnlo[ich]; // correct the nnlo factor for Higgs normalization 
	    if(txt[ich]=="ww")acc *=xww; // assign additional 10% uncertainties 
	    float xb[NBKGH]={0.}; // background normalization
	    for( Int_t j=0; j<NBKGH; ++j){ 	   
	      
	    if(j<(NBKGH-1)&&fabs(sbkg[j][im][ich])>0.){
	      xb[j] =gbkg[j][ich];
	      if(j>1)xb[j] *=geffb[1][ich]; // btag
	    }
	    else if(j==(NBKGH-1)&&fabs(sbkg[2][im][ich])>0.){ 
	      xb[j] =gbkg[2][ich]*gtagc; //ctag
	    }
	    if(j==3) xb[j] *=gcomm[3]; //common top xsec 
	    if(j==4) xb[j] *=gcomm[4]; // common diboson 
	    if(j>1)xb[j] *=gcommb_lum*geffb[0][ich]*geffb[2][ich]*geffb[3][ich]*geffb[4][ich]; // add additional systematic:lum,lep id, jes, isr/fsr+pdf   
	    }
	    // No SYS applied here 
	    if(NoSYS==3||NoSYS==1)acc = sfnnlo[ich]; 
	    if(NoSYS==3||NoSYS==2){
	      for(Int_t j =0; j<NBKGH; ++j){
		if(xb[j]>0)xb[j] = 1.;
	      }
	    }
	    // convert 1-d to 2d; 
	    TH2F* xd2d =(TH2F*)XMC[ich]; 
	    TH2F* xh2d =(TH2F*)XHIGGS[im][ich]; 
	    TH2F* xb2d[NBKGH]; 
	    for(Int_t i=0; i<NBKGH; ++i) xb2d[i] =(TH2F*)XBKG[i][im][ich]; 

	    for(Int_t k =1; k<xd2d->GetNbinsX()+1; ++k){ 
	      for(Int_t k1=1; k1<xd2d->GetNbinsY()+1; ++k1){  
		
		Double_t mc = 0.; 
		mc += tXsec*acc*xh2d->GetBinContent(k,k1); 
		efftot +=acc*xh2d->GetBinContent(k,k1);

		//background 
		for( Int_t j=0; j<NBKGH; ++j){ 
		  float a1 = xb2d[j]->GetBinContent(k,k1); 
		  if(NoSYS==11){ 
		    if(j==2)a1 =a1*xsampling +(1-xsampling)*XZHshape[0]->GetBinContent(k,k1); 
		    if(j==5)a1 =a1*xsampling+(1-xsampling)*XZHshape[1]->GetBinContent(k,k1); 
		  }
		  mc +=xb[j]*a1; // xb[j]==0, the hist is not defined
		}
		Int_t nx = (Int_t) xd2d->GetBinContent(k,k1); 
		if(mc>0) loglhS +=nx*TMath::Log(mc)-mc-lgm[nx]; 
	      } //loop over the histogram bins x 
	    } //
	  }
	  //*****************************************************************************
          // liklihood for D0 channels 
          //*****************************************************************************
	  else{ // for d0 channels 
	  // eff 

	  Double_t acc(1.); 
	  
	  for( Int_t i=0; i<NACCD0; ++i) {
	    acc *=geff_d0[i][ich-NCHA];
	      }
	  acc *=gcomm[0]; // common luminosity 
	  if(txt[ich].Contains("d0ww"))acc *=xww; // assign additional 10% uncertainties 

	  Double_t acc_bkg(1.0); 
	  for( Int_t i=0; i<NACCD0; ++i) {
	    acc_bkg *=gbkg_d0[i][ich-NCHA];
	      }
	  acc_bkg *=gcommb_lum; // common luminosity 
	  
	  float xb[NBKGHD0]={0.}; // background normalization 
	  // signal id
	  float eids=geff_d0[8][ich-NCHA]*geff_d0[9][ich-NCHA]*geff_d0[10][ich-NCHA]; 
	  float mids=geff_d0[12][ich-NCHA]*geff_d0[13][ich-NCHA];
	  float bids = geff_d0[5][ich-NCHA]*geff_d0[6][ich-NCHA]; 
	  // background id
	  float eid=gbkg_d0[8][ich-NCHA]*gbkg_d0[9][ich-NCHA]*gbkg_d0[10][ich-NCHA]; 
	  float mid=gbkg_d0[12][ich-NCHA]*gbkg_d0[13][ich-NCHA];
          float bid =gbkg_d0[5][ich-NCHA]*gbkg_d0[6][ich-NCHA]; 	   
	  // dileptons extra id (ee,mm), double btags; ignore for now by reset to 1.  
	  eids =1.; 
	  mids =1.; 
	  bids =1.; 
	  eid = 1.; 
	  mid = 1.; 
	  bid = 1.; 
	  if(txt[ich].Contains("d0wwee")){
	    acc *=eids;
	    for( Int_t j=0; j<NBKGHD0; ++j){ 	   
	      if(j <3){
		xb[j] = eid*acc_bkg*gcomm[4]; //*Xsec_EM
	      }
	      else if( j==3){
		xb[j] = eid*acc_bkg*gcomm[3]; // *Xsec_top
	      }
	      else if(j<NBKGHD0-1){
		xb[j] = eid*acc_bkg*gcomm[4]; //*Xsec_EM 		
	      }
	      else{
		xb[j] = eid*acc_bkg*gbkg_d0[22][ich-NCHA]; //*Xsec_qcd 	
	      }
	    }	    
	  }
	  else if(txt[ich].Contains("d0wwem")){
	    for( Int_t j=0; j<NBKGHD0; ++j){ 	   
	      if(j <3){
		xb[j] = acc_bkg*gcomm[4]; //*Xsec_EM
	      }
	      else if( j==3){
		xb[j] = acc_bkg*gcomm[3]; // *Xsec_top
	      }
	      else if(j<NBKGHD0-1){ 
		xb[j] = acc_bkg*gcomm[4]; //*Xsec_EM 	
	      }
	      else{
		xb[j] = acc_bkg*gbkg_d0[22][ich-NCHA]; //*Xsec_qcd 
	      }
	    }	    
	  }
	  else if(txt[ich].Contains("d0wwmm")){
	    acc *=mids;
	    for( Int_t j=0; j<NBKGHD0-2; ++j){ 	   
	      if(j==0||j==3){
		xb[j] = mid*acc_bkg*gcomm[4]; //*Xsec_EM
	      }
	      else if(j==1){ 
		xb[j] = mid*acc_bkg*gbkg_d0[20][ich-NCHA]; //*Xsec_HF
	      }
	      else if(j==2){ 
		xb[j] = mid*acc_bkg*gbkg_d0[21][ich-NCHA]; //*Xsec_LF
	      }
              else if(j==4){ 
		xb[j] = mid*acc_bkg*gcomm[3]; // *Xsec_top
	      }
	      else{
		xb[j] = mid*acc_bkg*gbkg_d0[22][ich-NCHA]; //*Xsec_qcd 	
	      }
	    }	    
	  }
	  else if(txt[ich].Contains("d0wh")){
	    for( Int_t j=0; j<4; ++j){ 	   
	      if(j<2){
		xb[j] = acc_bkg*gcomm[4]; //*Xsec_EM
	      }
	      else{
		xb[j] = acc_bkg*gbkg_d0[22][ich-NCHA]; //*Xsec_qcd 	
	      }
	    }	    
	  }
 	  else if((txt[ich]=="d0zh")){ // d0zheebb or d0zhmmbb 2tags
 	    for( Int_t j=0; j<NBKGHD0-2; ++j){ 	   
 	      if(j <2){
 		xb[j] = acc_bkg*gcomm[4]; //eid*Lum*trigger*Xsec_EM
 	      }
 	      else if(j<3){
		xb[j] = acc_bkg*gbkg_d0[20][ich-NCHA]; //*Xsec_HF
 	      }
 	      else if( j<4){
 		xb[j] = acc_bkg*gbkg_d0[19][ich-NCHA]; //*Xsec_bTag_LF
 	      }
 	      else if( j<5){
 		xb[j] = acc_bkg*gcomm[3]; // *Xsec_top
 	      }
 	      else{ 
 		xb[j] = acc_bkg*gbkg_d0[22][ich-NCHA]; //*Xsec_qcd 	
 	      }	      	    
 	    }
 	  }  
	  else if((txt[ich]=="d0evbb_s")||(txt[ich]=="d0mvbb_s")||(txt[ich]=="d0vvbb_s")){ // single tag
	    for( Int_t j=0; j<NBKGHD0-1; ++j){ 	   
	      if(j <1){
		xb[j] = acc_bkg*gcomm[4]; //eid*Lum*trigger*Xsec_EM
	      }
	      else if(j<2){
		xb[j] = acc_bkg*gbkg_d0[20][ich-NCHA]; //*Xsec_HF
	      }
	      else if( j<4){
		xb[j] = acc_bkg*gbkg_d0[19][ich-NCHA]; //*Xsec_bTag_LF
	      }
	      else if( j<6){
		xb[j] = acc_bkg*gcomm[3]; // *Xsec_top
	      }
	      else{ 
		xb[j] = acc_bkg*gbkg_d0[22][ich-NCHA]; //*Xsec_qcd 	
	      }	      	    
	    }
	  }  
	  else{ // double tags 
	    acc *=bids; 
	    for( Int_t j=0; j<NBKGHD0-1; ++j){ 	   
	      if(j <1){
		xb[j] = bid*acc_bkg*gcomm[4]; //*Xsec_EM
	      }
	      else if( j<2){ 
		xb[j] = bid*acc_bkg*gbkg_d0[20][ich-NCHA]; //*Xsec_HF
	      }
	      else if( j<4){
		xb[j] = bid*acc_bkg*gbkg_d0[19][ich-NCHA]; //*Xsec_bTag_LF
	      }
	      else if( j<6){
		xb[j] = bid*acc_bkg*gcomm[3]; // *Xsec_top
	      }
	      else{ 
		xb[j] = bid*acc_bkg*gbkg_d0[22][ich-NCHA]; //*Xsec_qcd 	
	      }	      	    
	    }
	  }  
	  //		
	  // No SYS applied here 

	    if(NoSYS==3||NoSYS==1)acc = sfnnlo[ich]; 
	    if(NoSYS==3||NoSYS==2){ 
	    for(Int_t j =0; j<NBKGHD0; ++j){
	      if(xb[j]>0)xb[j] = 1.;
	    }
	    }
	  for(Int_t k =1; k<XDATA[im][ich]->GetNbinsX()+1; ++k){ 

	    //if(txt[ich].Contains("d0wwmm")&&XDATA[im][ich]->GetBinCenter(k)>2.5)continue; // requiring dphi<2.5 for ww->mm channel to remove Z->mumu events 

	    Double_t mc = 0.; 
	    mc += tXsec*acc*XHIGGS[im][ich]->GetBinContent(k); 
	    efftot +=acc*XHIGGS[im][ich]->GetBinContent(k);

	    //background 
	    Int_t nx = NBKGHD0-1; 
	    if(txt[ich].Contains("d0wwmm")){
	      nx = NBKGHD0-2; 
	    }
	    else if(txt[ich].Contains("d0ww")){
	      nx = NBKGHD0;
	    }
	    if(txt[ich].Contains("d0wh"))nx = NBKGHD0-4; 
	    if(txt[ich].Contains("d0zh"))nx = NBKGHD0-2;
	    for( Int_t j=0; j<nx; ++j){
	      if(xb[j]>0.)mc +=xb[j]*XBKG_d0[j][im][ich-NCHA]->GetBinContent(k); // xb[j]==0, the hist is not defined
	    }
	    
	    nx = (Int_t) XMC[ich]->GetBinContent(k); 
	    if(mc>0)loglhS +=nx*TMath::Log(mc)-mc-lgm[nx]; 
	    
	  } //loop over the histogram bins
	  }
	} // channels in consideration 	
      } // loop all the channels
      lh[tBin] +=efftot*TMath::Exp(loglhS);
    }   // end of sampling loop
    // average of likelihood
    Double_t avelh [nXsecR] = {0.};
    Double_t sumlh(0.); 
    for ( Int_t i = 0 ; i < nXsecR ; i++ ) {
      if ( nc[i] > 0. ) {
	avelh[i] = lh[i]/nc[i];
	sumlh +=avelh[i]; 
      }
    }
    // normalize
    Double_t minchi(999.);
    Double_t avechi[nXsecR] = {0.};
    Int_t XsecMax(0); 
    for ( Int_t Xsec = 0 ; Xsec < nXsecR ; Xsec++ ) {
      avelh[Xsec] /= sumlh;
      if ( avelh[Xsec] > 0. ){
	avechi[Xsec] = -2.*TMath::Log(avelh[Xsec]);
      }
      else { avechi[Xsec] = 999.;}
      if ( minchi  > avechi[Xsec] ) {
	fit[im][0] = (Xsec+0.5)*step;
	minchi = avechi[Xsec];
	XsecMax = Xsec; 
      }

      hchi[im]->Fill(Xsec*step+step/2,avechi[Xsec]);
      hlh [im]->Fill(Xsec*step+step/2,avelh [Xsec]);
    }
    // compute the +- sigma
    Int_t di(0); 
    for(int i =XsecMax; i<nXsecR; ++i) {
      if((avechi[i]-avechi[XsecMax])>1.0){ 	
	di = i;
	break; 
      }
    }
    if(di==0) di = nXsecR; 
    fit[im][1] = (di-XsecMax)*step; 
    di=0;
    for(int i =XsecMax; i>-1; --i) { 
      if((avechi[i]-avechi[XsecMax])>1.0){ 
	di = i;
	break; 
      }
    }
    fit[im][2] = (XsecMax-di)*step; 

    // 95% C.L. upper limit
    Double_t A = 0.;
    for ( Int_t Xsec = 0 ; Xsec < nXsecR ; Xsec++ ) {
      A += avelh[Xsec];
      if ( A > 0.95 ) {
	limit[im] = Xsec*step+step/2;
	break;
      }
    }

    if(ie%100==0)std::cout <<" Pseudo experiments: "<<ie
			   <<" Mass="<< HMass[im]
			   <<" XsecR="<< fit  [im][0]
			   <<"+"<< fit  [im][1]
			   <<"-"<< fit  [im][2]	
			   << " 95% CL limit ="<<limit[im]
			   << std::endl;
  }

  // set 95% C.L. upper limit line for likelihood

  for(Int_t i=0; i<NMASSX; ++i){ 
    hl[i]->Fill(limit[i]); 
  }

  gROOT->Time();
  } // loop over all the experiments 

  // draw histograms
  TString mytxt("_"); 
  if(nchannel==NCHAT){ 
    mytxt += "all"; }
  else{
    for(Int_t i=0; i<NCHAT; ++i){ 
      if(icha[i]>0){
	mytxt +="1";
      }
      else{
	mytxt +="0";
      }
    }
  }
  TString psFile("R_combined_limit_pseudoexperiments_cdfd0_new");
  if(NoSYS==3)psFile +="_NoSys";
  if(NoSYS==1)psFile +="_NoSysEff";
  if(NoSYS==2)psFile +="_NoSysBkg";
  if(NoSYS==4)psFile +="_NoCor"; // no correlation between signal and backgrounds
  if(NoSYS==10)psFile +="_shape"; // ZH->llbb with different zh shape 
  if(NoSYS==11)psFile +="_combshape"; // ZH->llbb with different zh shape with alternative combined shape 
  if(NoSYS==20)psFile +="_uncorr"; // treat most of systematic uncorrelated, except lum, btag, top ew xsec   
  if(NoSYS>100){
    psFile +="_mass";
    psFile +=(Int_t)NoSYS;
  }

  psFile +=mytxt; 
  psFile += ".ps";

  TCanvas cx("cx","cx",640,640);
  TPostScript ps(psFile);

    cx.Clear();
    Int_t nx0 =(NMASSX+1)/2; 
    cx.Divide(nx0,2);
    for ( Int_t i = 0 ; i < NMASSX ; i++ ) {
      cx.cd(i+1);
      hl[i]->SetLineColor(1);
      hl[i]->Draw();
    }
    cx.Update();
 
  ps.Close();
  cx.Close();

  // output data file
  TString dataFile("R_combined_limit_pseudoexperiments_cdfd0_new");
  if(NoSYS==3)dataFile +="_NoSys";
  if(NoSYS==1)dataFile +="_NoSysEff";
  if(NoSYS==2)dataFile +="_NoSysBkg";
  if(NoSYS==4)dataFile +="_NoCor"; // no correlation between signal and backgrounds
  if(NoSYS==10)dataFile +="_shape"; // ZH->llbb with different zh shape 
  if(NoSYS==11)dataFile +="_combshape"; // ZH->llbb with different zh shape with alternative combined shape 
  if(NoSYS==20)dataFile +="_uncorr"; // treat most of systematic uncorrelated, except lum, btag, top ew xsec   
  if(NoSYS>100){
    dataFile +="_mass";
    dataFile +=(Int_t)NoSYS;
  }

  dataFile += mytxt; 
  dataFile += ".dat";
  ofstream fout(dataFile);
  fout  <<" Channels Combined="<<nchannel<<std::endl; 
  if(NoSYS==1||NoSYS==2||NoSYS==3) fout<<" With no systematic included "<<" "<<NoSYS<<std::endl; 

  for ( Int_t i = 0 ; i < NMASSX ; i++ ) {

    fout <<" Mass="<< HMass[i]
	 <<" Mean of 95% CL="<< hl[i]->GetMean(1) // 1=x axis
	 <<" RMS "<< hl[i]->GetRMS(1)
	 << std::endl;
    for( Int_t ix=0; ix<NCHAT; ix++){ 
      if(hevt[ix][i]->Integral()>0){
	fout<<" mass "<< HMass[i]<<" "<<txt[ix]<<" events="<<hevt[ix][i]->GetMean(1)
	    <<" RMS="<<hevt[ix][i]->GetRMS(1)<<std::endl;
      }
    }
  }
  fout.close();

  // output root file
  TString rootFile("R_combined_limit_pseudoexperiments_cdfd0_new");
  if(NoSYS==3)rootFile +="_NoSys";
  if(NoSYS==1)rootFile +="_NoSysEff";
  if(NoSYS==2)rootFile +="_NoSysBkg";
  if(NoSYS==4)rootFile +="_NoCor"; // no correlation between signal and backgrounds
  if(NoSYS==10)rootFile +="_shape"; // ZH->llbb with different zh shape 
  if(NoSYS==11)rootFile +="_combshape"; // ZH->llbb with different zh shape with alternative combined shape 
  if(NoSYS==20)rootFile +="_uncorr"; // treat most of systematic uncorrelated, except lum, btag, top ew xsec   
  if(NoSYS>100){
    rootFile +="_mass";
    rootFile +=(Int_t)NoSYS;
  }

  rootFile += mytxt; 
  rootFile += ".root";
  TFile froot(rootFile,"recreate");
  for ( Int_t i = 0 ; i < NMASSX ; i++ ) {
      hl [i]->Write();
      for(Int_t j=0; j<NCHAT; j++){
	if(hevt[j][i]->Integral()>0)hevt[j][i]->Write(); 
      }
  }

  froot.Close();

  // delete histograms
  for ( Int_t i = 0 ; i < NMASSX ; i++ ) {
      hchi[i]->Delete();
      hlh [i] ->Delete();
      hl[i]->Delete();
      for(Int_t j=0; j<NCHAT; j++)hevt[j][i]->Delete(); 
  }
}