#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(); // calculate the combined limits 
void higgs_combined_limit_pseudoexperiment( ); // doing pseudo-experiments 

//Int_t Nbins = 0; 
static const bool debug = false; // for debug purpose 
static const Int_t NCHA = 4; // number of channels to combine (lnubb_s,lnubb_d,nunubb,ww)  
static const Int_t NACC = 6; // source of systematic for acceptance(lum/btag/lepton/JES/I(F)SR+PDF/Other)
static const Int_t NBKG = 6; // source of backgrounds(HF/mistag/top/nonW/diboson/other) 
static const Int_t NMASS = 10; // Higgs mass 110-200 with 10 GeV steps, we can extrapolate in 5 GeV  
static Int_t icha[NCHA]={0}; // the channel to include in combination 
static float seff[NACC][NCHA]={0.}; // values of acceptance systematic  
static float sbkg[NBKG][NMASS][NCHA]={0.}; // values of background systematic  
static const Int_t Nmax[NCHA] ={5, 5, 3, 10}; 
static const float sfnnlo[NCHA]={0.96,0.96,0.99*1.3,1.45}; // 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][NCHA] = {0}; // the point to data histograms 
TH1F* XBKG[NBKG][NMASS][NCHA] = {0}; // the point to background histograms 
TH1F* XHIGGS[NMASS][NCHA] ={0}; // higgs signal 
TH1F* XMC[NCHA] = {0}; // histogram of 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( )
{
  // 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[NCHA] ={"lnubb_s","lnubb_d","nunubb","ww"};
  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][NCHA-1]);
    if(debug)std::cout<<" sbkg["<<NBKG-1<<"]["<<i<<"]["<<NCHA-1<<"]="<<sbkg[NBKG-1][i][NCHA-1]<<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]; 
      }
    }
  }
  // 
  TFile* f[2];
  f[0] = new TFile("shape_snntag.root","read"); 
  f[1] = new TFile("shape_secvtx.root","read");
  TFile* f3 = new TFile("shape_nunubb.root","read");   
  TFile* f4 = new TFile("shape_hww.root","read"); 
  
  for(int i=0; i<2; ++i){ 
    XDATA[0][i]  = (TH1F*)f[i]->Get("M2j_data"); 
    XHIGGS[0][i] = (TH1F*)f[i]->Get("M2j_WH110_CTEQ5L_ISRon_FSRon"); 
    XHIGGS[1][i] = (TH1F*)f[i]->Get("M2j_WH120_CTEQ5L_ISRon_FSRon"); 
    XHIGGS[2][i] = (TH1F*)f[i]->Get("M2j_WH130_CTEQ5L_ISRon_FSRon"); 
    XHIGGS[3][i] = (TH1F*)f[i]->Get("M2j_WH140_CTEQ5L_ISRon_FSRon"); 
    XHIGGS[4][i] = (TH1F*)f[i]->Get("M2j_WH150_CTEQ5L_ISRon_FSRon"); 
    XBKG[0][0][i] = (TH1F*)f[i]->Get("M2j_HF");     
    XBKG[1][0][i] = (TH1F*)f[i]->Get("M2j_mistag"); 
    XBKG[2][0][i] = (TH1F*)f[i]->Get("M2j_Top"); 
    XBKG[3][0][i] = (TH1F*)f[i]->Get("M2j_nonW"); 
    XBKG[4][0][i] = (TH1F*)f[i]->Get("M2j_diboson");

    for(int j =1; j<NMASS; ++j){
      XDATA[j][i] = XDATA[0][i]; 
      for(int k=0; k<NBKG; ++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[1][2] = (TH1F*)f3->Get("Dijet_mh120"); 
    XHIGGS[2][2] = (TH1F*)f3->Get("Dijet_mh130"); 
    XBKG[0][0][2] = (TH1F*)f3->Get("Dijet_hf");     
    XBKG[1][0][2] = (TH1F*)f3->Get("Dijet_mistag"); 
    XBKG[2][0][2] = (TH1F*)f3->Get("Dijet_top"); 
    XBKG[3][0][2] = (TH1F*)f3->Get("Dijet_qcd"); 
    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<NBKG; ++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 = 110+i*10; 
    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; 
  }

  const Int_t nXsecR = 600;
  const Int_t nIter = 2000000; // 100k or 500k 

  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<NCHA; ++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];
  TH1F *hl = new TH1F("Limit","Limit", NMASSX, 105., 105+NMASSX*10);
  for ( Int_t i = 0 ; i < NMASSX ; i++ ) {
      TString nchi("hchi_");
      int j = 110 + i*10; 
      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<< " " <<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++ ) {

    // 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 geff[NACC]; 
      float gbkg[NBKG]; 
      // modeling efficiency >0

      for(int i=0; i<NACC; ++i){ 
	float val; 
	float v1; 
	do{
	  geff[i] = gRand.Gaus(0.,1.0);
	  val =1.0; 
	  for(int j =0; j<NCHA; ++j){ 
	    if(icha[j]>0&&seff[i][j]>0.){ 
	      v1 =(1.0+geff[i]*seff[i][j]);
	      if(v1<=0.)val = v1; 
	    }
	  }
	}while(val<=0.); 
      }
      // modeling background >0 
      for(int i=0; i<NBKG; ++i){ 
	float val;
	float v1;
	do{
	  gbkg[i] = gRand.Gaus(0.,1.0);
	  val =1.0;
	  for(int j =0; j<NCHA; ++j){ 
	    if(icha[j]>0&&sbkg[i][im][j]>0.){ 
	      v1 =(1.0+gbkg[i]*sbkg[i][im][j]); 
	      if(v1<=0.)val = v1;
	    }
	  }
	}while(val<=0.); 
      }

      // get cross section Ratio from 0 to nXsecR*step  
      Double_t tXsec = gRand.Uniform()*nXsecR*step; 
      Int_t tBin = (Int_t)(tXsec/step);
      nc[tBin]++;


      // likelihood calculation
      Double_t loglhS = 0.;
      Double_t efftot = 0.; 
      for( Int_t ich=0; ich<NCHA; ++ich){ 
	if(icha[ich]>0&&im<Nmax[ich]){
	  // eff 

	  Double_t acc(1.); 
	  for( Int_t i=0; i<NACC; ++i) {
	    if(seff[i][ich]<0.){ // acc systematic treated uncorrected  
	      acc *=getGauss(gRand, 1., fabs(seff[i][ich])); 
	    }
	    else if(seff[i][ich]>0.){ 
	      acc *=(1.+geff[i]*seff[i][ich] ); 
	    } 
	  }
	  acc *=sfnnlo[ich]; // correct the nnlo factor for Higgs normalization 
	  if(txt[ich]=="ww")acc *=getGauss(gRand, 1., 0.1); // assign additional 10% uncertainties 
	  float xb[NBKG]; // background normalization 
	  for( Int_t j=0; j<NBKG; ++j){ 	   
	    xb[j] = 0.; 
	    if(sbkg[j][im][ich]>0.){
	      xb[j]=(1+gbkg[j]*sbkg[j][im][ich]); }
	    else if( sbkg[j][im][ich]<0.){ // background systematic treated uncorrected  
	      xb[j]=getGauss(gRand, 1., fabs(sbkg[j][im][ich])); }
	  }
	  
	  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<NBKG; ++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 
	} // 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="<< 110+im*10
	     <<" 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->SetBinContent(i+1, limit[i]); 
    hl->SetBinError(i+1,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 Run II Preliminary (695 pb^{-1})");
  pPreliminary.SetNDC(true);
  pPreliminary.SetTextSize(0.05);

  // draw histograms
  TString mytxt("_"); 
  if(nchannel==NCHA){ 
    mytxt += "all"; }
  else{
    for(Int_t i=0; i<NCHA; ++i){ 
      if(icha[i]>0){
	mytxt +=txt[i];
	mytxt +="_"; 
      }
    }
  }
  TString psFile("R_combined_limit");
  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 = 110+i*10;
      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);
      hlh[i]->Draw();
      line[i].Draw();
    }
    cx.Update();
    
    cx.Clear(); 
    hl->Draw();
    cx.Update(); 

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

  // output data file
  TString dataFile("R_combined_limit");
  dataFile += mytxt; 
  dataFile += ".dat";
  ofstream fout(dataFile);
  fout <<" Channels Combined="<<nchannel<<std::endl; 

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

    fout <<" Mass="<< 110+i*10
	 <<" 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");
  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( )
{
  // 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[NCHA] ={"lnubb_s","lnubb_d","nunubb","ww"};
  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][NCHA-1]);
    if(debug)std::cout<<" sbkg["<<NBKG-1<<"]["<<i<<"]["<<NCHA-1<<"]="<<sbkg[NBKG-1][i][NCHA-1]<<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]; 
      }
    }
  }
  // 
  TFile* f[2];
  f[0] = new TFile("shape_snntag.root","read"); 
  f[1] = new TFile("shape_secvtx.root","read");
  TFile* f3 = new TFile("shape_nunubb.root","read");   
  TFile* f4 = new TFile("shape_hww.root","read"); 
  
  for(int i=0; i<2; ++i){ 
    XDATA[0][i]  = (TH1F*)f[i]->Get("M2j_data"); 
    XDATA[0][i]->Rebin();
    XHIGGS[0][i] = (TH1F*)f[i]->Get("M2j_WH110_CTEQ5L_ISRon_FSRon");
    XHIGGS[1][i] = (TH1F*)f[i]->Get("M2j_WH120_CTEQ5L_ISRon_FSRon");
    XHIGGS[2][i] = (TH1F*)f[i]->Get("M2j_WH130_CTEQ5L_ISRon_FSRon"); 
    XHIGGS[3][i] = (TH1F*)f[i]->Get("M2j_WH140_CTEQ5L_ISRon_FSRon"); 
    XHIGGS[4][i] = (TH1F*)f[i]->Get("M2j_WH150_CTEQ5L_ISRon_FSRon"); 
    XBKG[0][0][i] = (TH1F*)f[i]->Get("M2j_HF");     
    XBKG[1][0][i] = (TH1F*)f[i]->Get("M2j_mistag"); 
    XBKG[2][0][i] = (TH1F*)f[i]->Get("M2j_Top"); 
    XBKG[3][0][i] = (TH1F*)f[i]->Get("M2j_nonW"); 
    XBKG[4][0][i] = (TH1F*)f[i]->Get("M2j_diboson");
    for(int j=0; j<5; ++j){  // rebin to be consistent 
      XHIGGS[j][i]->Rebin(); 
      XBKG[j][0][i]->Rebin(); 
    }
    for(int j =1; j<NMASS; ++j){ // copy to the other mass points 
      XDATA[j][i] = XDATA[0][i]; 
      for(int k=0; k<NBKG; ++k){
	XBKG[k][j][i] = XBKG[k][0][i];
      }
    }
  }
  // what about ZH->nunubb ?   
    XDATA[0][2]  = (TH1F*)f3->Get("Dijet_data"); 
    XDATA[0][2]->Rebin(); 
    XHIGGS[0][2] = (TH1F*)f3->Get("Dijet_mh110");
    XHIGGS[0][2]->Rebin();
    XHIGGS[1][2] = (TH1F*)f3->Get("Dijet_mh120");
    XHIGGS[1][2]->Rebin();
    XHIGGS[2][2] = (TH1F*)f3->Get("Dijet_mh130");
    XHIGGS[2][2]->Rebin();
    XBKG[0][0][2] = (TH1F*)f3->Get("Dijet_hf");     
    XBKG[1][0][2] = (TH1F*)f3->Get("Dijet_mistag"); 
    XBKG[2][0][2] = (TH1F*)f3->Get("Dijet_top"); 
    XBKG[3][0][2] = (TH1F*)f3->Get("Dijet_qcd"); 
    XBKG[4][0][2] = (TH1F*)f3->Get("Dijet_diboson");

    for(int j=0; j<5; ++j){
      XBKG[j][0][2]->Rebin(); 
    }

    for(int j =1; j<NMASS; ++j){
      XDATA[j][2] = XDATA[0][2]; 
      for(int k=0; k<NBKG; ++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 = 110+i*10; 
    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);
  }

  const Int_t nXsecR = 300;
  const Int_t nIter = 200000; // 100k or 500k 

  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<NCHA; ++i){ 
    if(icha[i]>0){
      nchannel++; 
      if(NMASSX<Nmax[i])NMASSX = Nmax[i]; 
    }
  }
   
  // make copy of data 
  for(Int_t i =0; i<NCHA; ++i){ 
    XMC[i] = (TH1F*)XDATA[0][i]->Clone(); 
  }
  // combined chi^{2} and likelihood histograms

  TH1F *hchi[NMASS];
  TH1F *hlh [NMASS];
  TH1F *hl[NMASS]; 
  TH1F *hevt[NCHA][NMASS]; 
  for ( Int_t i = 0 ; i < NMASSX ; i++ ) {
      TString nchi("hchi_");
      int j = 110 + i*10; 
      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<NCHA; ++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; 
  // 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
  float geff[NACC]; 
  float gbkg[NBKG]; 
  for ( Int_t im = 0 ; im < NMASSX ; im++ ) {
    hchi[im]->Reset(); 
    hlh[im]->Reset(); 
    // make data: 
    for(int i=0; i<NBKG; ++i){ 
      float val;
      float v1;
      do{
	gbkg[i] = gRand.Gaus(0.,1.0);
	val = 1.0;
	for(int j =0; j<NCHA; ++j){ 
	  if(icha[j]>0&&sbkg[i][im][j]>0.){ 
	    v1 =(1.0+gbkg[i]*sbkg[i][im][j]); 
	    if(v1<=0.)val = v1; 
	  }
	}
      }while(val<=0.); 
    }    
    for( Int_t ich=0; ich<NCHA; ++ich){      
      XMC[ich]->Reset(); 
      if(icha[ich]>0&&im<Nmax[ich]){
	float xb[NBKG]; // background normalization 
	for( Int_t j=0; j<NBKG; ++j){
	  xb[j]=0.;
	  if(sbkg[j][im][ich]>0.){
	    xb[j]=(1+gbkg[j]*sbkg[j][im][ich]); }
	  else if( sbkg[j][im][ich]<0.){ // background systematic treated uncorrected  
	    xb[j]=getGauss(gRand, 1., fabs(sbkg[j][im][ich])); }
	}
	for(Int_t k =1; k<XDATA[im][ich]->GetNbinsX()+1; ++k){ 	  
	    Double_t mc = 0.; 
	    //background 
	    for( Int_t j=0; j<NBKG; ++j){
	      if(xb[j]>0.)mc +=xb[j]*XBKG[j][im][ich]->GetBinContent(k);
	    }
	    Int_t nx = gRand.Poisson(mc); 
	    XMC[ich]->SetBinContent(k, nx); 
	    XMC[ich]->SetBinError(k, sqrt(nx)); 
	} // loop over 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 
      // modeling efficiency >0
      for(int i=0; i<NACC; ++i){ 
	float val; 
	float v1;
	do{
	  geff[i] = gRand.Gaus(0.,1.0);
	  val = 1.0;
	  for(int j =0; j<NCHA; ++j){ 
	    if(icha[j]>0&&seff[i][j]>0.){ 
	      v1=(1.0+geff[i]*seff[i][j]); 
	      if(v1<=0.)val = v1;
	    }
	  }
	}while(val<=0.); 
      }
      // modeling background >0 
      for(int i=0; i<NBKG; ++i){ 
	float val;
	float v1;
	do{
	  gbkg[i] = gRand.Gaus(0.,1.0);
	  val = 1.0;
	  for(int j =0; j<NCHA; ++j){ 
	    if(icha[j]>0&&sbkg[i][im][j]>0.){ 
	      v1 =(1.0+gbkg[i]*sbkg[i][im][j]);
	      if(v1<=0.)val = v1;
	    }
	  }
	}while(val<=0.); 
      }

      // get cross section Ratio from 0 to nXsecR*step  
      Double_t tXsec = gRand.Uniform()*nXsecR*step; 
      Int_t tBin = (Int_t)(tXsec/step);
      nc[tBin]++;

      // likelihood calculation
      Double_t loglhS = 0.;
      Double_t efftot = 0.; 
      for( Int_t ich=0; ich<NCHA; ++ich){ 
	if(icha[ich]>0&&im<Nmax[ich]){
	  // eff 

	  Double_t acc(1.); 
	  for( Int_t i=0; i<NACC; ++i) {
	    if(seff[i][ich]<0.){ // acc systematic treated uncorrected  
	      acc *=getGauss(gRand, 1., fabs(seff[i][ich])); 
	    }
	    else if(seff[i][ich]>0.){ 
	      acc *=(1.+geff[i]*seff[i][ich] ); 
	    } 
	  }
	  acc *=sfnnlo[ich]; // correct the nnlo factor for Higgs normalization 
	  if(txt[ich]=="ww")acc *=getGauss(gRand, 1., 0.1); // assign additional 10% uncertainties 

	  float xb[NBKG]; // background normalization 
	  for( Int_t j=0; j<NBKG; ++j){ 
	    xb[j] =0.;
	    if(sbkg[j][im][ich]>0.){
	      xb[j]=(1+gbkg[j]*sbkg[j][im][ich]); }
	    else if( sbkg[j][im][ich]<0.){ // background systematic treated uncorrected  
	      xb[j]=getGauss(gRand, 1., fabs(sbkg[j][im][ich])); }
	  }
	  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<NBKG; ++j){ 
	      if(xb[j]>0.)mc +=xb[j]*XBKG[j][im][ich]->GetBinContent(k);
	    }
	    //Int_t nx = 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 
	} // 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="<< 110+im*10
			   <<" 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==NCHA){ 
    mytxt += "all"; }
  else{
    for(Int_t i=0; i<NCHA; ++i){ 
      if(icha[i]>0){
	mytxt +=txt[i];
	mytxt +="_"; 
      }
    }
  }
  TString psFile("R_combined_limit_pseudoexperiments");
  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");
  dataFile += mytxt; 
  dataFile += ".dat";
  ofstream fout(dataFile);
  fout  <<" Channels Combined="<<nchannel<<std::endl; 

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

    fout <<" Mass="<< 110+i*10
	 <<" Mean of 95% CL="<< hl[i]->GetMean(1) // 1=x axis
	 <<" RMS "<< hl[i]->GetRMS(1)
	 << std::endl;
    for( Int_t ix=0; ix<NCHA; ix++){ 
      if(hevt[ix][i]->Integral()>0){
	fout<<" mass "<< 110+i*10<<" "<<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");
  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<NCHA; 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<NCHA; j++)hevt[j][i]->Delete(); 
  }

}