#ifndef __CINT__
#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"
#include "TMath.h"
#include <sys/time.h>

#define MAX(x,y)	((x) > (y) ? (x) : (y))
#define MIN(x,y)	((x) < (y) ? (x) : (y))

Double_t GetMedian(TH1F* th1);
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 technicolor_limit_cdf(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=30+Luminosity rescale; NoSYS>100 for individual masses
void technicolor_limit_pseudoexperiment_cdf(Int_t NoSYS=0); // doing pseudo-experiments 

//Int_t Nbins = 0; 
static const bool debug = true; // for debug purpose 
static const Int_t NCHA = 6; // 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 = 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 NMASS = 9; // Higgs mass 110-180 with 10 GeV steps, we can extrapolate in 5 GeV  
static const Int_t RhoMass[NMASS]={180,190,200,200,210,210,220,220,220}; // higgs mass for this 10 points
static const Int_t PiMass[NMASS]={95,105,105,115,115,125,115,125,135}; // higgs mass for this 10 points
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] ={NMASS, NMASS, NMASS};
static const Int_t Npseudo=1000; // 1000 doing psudo-experiment 
TH1F* XDATA[NMASS][NCHA] = {{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* XHIGGS[NMASS][NCHA] ={{0}}; // higgs signal
TH1F* XMC[NCHA] ={0}; // histograms for pseudo-experiments 
static Double_t lgm[10000]={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 technicolor_limit_cdf(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);
  // dt1=st-st; dt2=st-jpb
  TString txt[NCHA] ={"lnubb_dt1c","lnubb_dt2c","lnubb_stnnc","lnubb_dt1n","lnubb_dt2n","lnubb_stnnn"};

  Double_t nfact(0);
  Double_t klum(1.0); // luminosity rescaling from the current 1 fb^-1 
  if(NoSYS>30&&NoSYS<40)klum = NoSYS-30;
  lgm[0] = 0.; 
  for(Int_t i =1; i<10000; ++i){ 
    nfact +=TMath::Log(i); 
    lgm[i] = nfact; 
  }
  char junk[120] =" ";
  FILE * infile = fopen("technicolor_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,120, 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,120, 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; 
    if(NoSYS==20)sbkg[j][0][i] = -1.*fabs(sbkg[j][0][i]); 
  }
  fgets(junk,120, infile);
  }
  
  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]; 
      }
    }
  }


  //get cdf data wh->lnubb with nnoutput
   TString Chanx[3]={"STST","STJP","1tagNN"};
   TString Chans[NMASS]={"1a","2b","3b","3c","4c","4d","5c","5d","5e"};
   for(int i =0; i<3; ++i){ 
      for(int j = 0; j<Nmax[i]; ++j){ 
	TFile* f = new TFile(Chanx[i]+"_ALL.root","read");
	TString Higgsx = "h_Mjj_vs_Qval_"+Chanx[i]; 
	XHIGGS[j][i] =(TH1F*)f->Get(Higgsx+"_cexo"+Chans[j]+"_2jet");
	XDATA[j][i] = (TH1F*)f->Get(Higgsx+"_Data_2jet"); 
	XBKG[0][j][i] = (TH1F*)f->Get(Higgsx+"_W+npAll_2jet"); 
	XBKG[1][j][i] = (TH1F*)f->Get(Higgsx+"_AE_2jet"); 
	XBKG[2][j][i] = (TH1F*)f->Get(Higgsx+"_WbbAll_2jet");
	XBKG[2][j][i]->Add((TH1F*)f->Get(Higgsx+"_WccAll_2jet"));
	XBKG[3][j][i] = (TH1F*)f->Get(Higgsx+"_ttop75_2jet");
	XBKG[3][j][i]->Add((TH1F*)f->Get(Higgsx+"_stop-s_2jet"));
	XBKG[3][j][i]->Add((TH1F*)f->Get(Higgsx+"_stop-t_2jet")); 
	XBKG[4][j][i] = (TH1F*)f->Get(Higgsx+"_WW_2jet"); 
	XBKG[4][j][i]->Add((TH1F*)f->Get(Higgsx+"_WZ_2jet"));
	XBKG[4][j][i]->Add((TH1F*)f->Get(Higgsx+"_ZZ_2jet"));
	XBKG[4][j][i]->Add((TH1F*)f->Get(Higgsx+"_ZtauAll_2jet"));
	if(j==0){ 
	  std::cout<<Chanx[i]<<" Data "<<XDATA[j][i]->Integral()<<" Charged Signal "<<XHIGGS[j][i]->Integral()<<" Bkg0 "<<XBKG[0][j][i]->Integral()<<" bkg1 "<<XBKG[1][j][i]->Integral()<<" Bkg2 "<<XBKG[2][j][i]->Integral()<<" Bkg3 "<<XBKG[3][j][i]->Integral()<<" Bkg4 "<<XBKG[4][j][i]->Integral()<<std::endl;
	}
      }
   }

   for(int i =0; i<3; ++i){ 
      for(int j = 0; j<Nmax[i]; ++j){ 
	TFile* f = new TFile(Chanx[i]+"_ALL.root","read");
	TString Higgsx = "h_Mjj_vs_Qval_"+Chanx[i]; 
	XHIGGS[j][i+3] =(TH1F*)f->Get(Higgsx+"_nexo"+Chans[j]+"_2jet");
	XDATA[j][i+3] = (TH1F*)f->Get(Higgsx+"_Data_2jet"); 
	XBKG[0][j][i+3] = (TH1F*)f->Get(Higgsx+"_W+npAll_2jet"); 
	XBKG[1][j][i+3] = (TH1F*)f->Get(Higgsx+"_AE_2jet"); 
	XBKG[2][j][i+3] = (TH1F*)f->Get(Higgsx+"_WbbAll_2jet");
	XBKG[2][j][i+3]->Add((TH1F*)f->Get(Higgsx+"_WccAll_2jet"));
	XBKG[3][j][i+3] = (TH1F*)f->Get(Higgsx+"_ttop75_2jet");
	XBKG[3][j][i+3]->Add((TH1F*)f->Get(Higgsx+"_stop-s_2jet")); 
	XBKG[3][j][i+3]->Add((TH1F*)f->Get(Higgsx+"_stop-t_2jet")); 
	XBKG[4][j][i+3] = (TH1F*)f->Get(Higgsx+"_WW_2jet"); 
	XBKG[4][j][i+3]->Add((TH1F*)f->Get(Higgsx+"_WZ_2jet"));
	XBKG[4][j][i+3]->Add((TH1F*)f->Get(Higgsx+"_ZZ_2jet"));
	XBKG[4][j][i+3]->Add((TH1F*)f->Get(Higgsx+"_ZtauAll_2jet"));
	if(j==1){ 
	  std::cout<<Chanx[i]<<" Data "<<XDATA[j][i+3]->Integral()<<" Neutral Signal "<<XHIGGS[j][i+3]->Integral()<<" Bkg0 "<<XBKG[0][j][i+3]->Integral()<<" bkg1 "<<XBKG[1][j][i+3]->Integral()<<" Bkg2 "<<XBKG[2][j][i+3]->Integral()<<" Bkg3 "<<XBKG[3][j][i+3]->Integral()<<" Bkg4 "<<XBKG[4][j][i+3]->Integral()<<std::endl;
	}
      }
   }

   std::cout<<" finishing reading WH "<<std::endl; 

  // Rescaling the luminosity for future sensitivity studies  
  for(int i=0; i<NCHA; ++i){ 
    for(int j=0; j<Nmax[i]; ++j){ 
      if(XDATA[j][i])XDATA[j][i]->Scale(klum);
      if(XHIGGS[j][i])XHIGGS[j][i]->Scale(klum);
	for(int k=0; k<NBKGH; ++k){ 
	  if(XBKG[k][j][i]){
	    XBKG[k][j][i]->Scale(klum);
	  }
      }
    }
  }


  /////////////////////

  const Int_t nXsecR = 100;
  static Int_t nXsecRx = nXsecR;
  static Int_t nIter =1000000; //3000000; // 100k or 500k 
  const float step =0.05; 
  // 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]; 
    }
  }
   
  if(NoSYS==3) nIter = 2*nXsecRx; 
  // save the histograms for efficiency and background uncertainties 
  TH1F* heff[NCHA][3]; 
  for(Int_t i=0; i<NCHA; ++i){ 
    if(icha[i]>0){ 
      TString txte("Eff_"); 
      txte +=txt[i];
      TString txtb("BKG_"); 
      txtb +=txt[i];
      TString txtc("ToTBKG_"); 
      txtc +=txt[i];
      heff[i][0] = new TH1F(txte, txte, 50, 0., 3.); 
      heff[i][1] = new TH1F(txtb, txtb, 50, 0., 3.); 
      heff[i][2] = new TH1F(txtc, txtc, 500, 0., 2000.); 
    }
  }      
  // combined chi^{2} and likelihood histograms
  
  TH1F *hchi[NMASS];
  TH1F *hlh [NMASS];
  Int_t nmbin = NMASS*2; 
  TH1F *hl = new TH1F("Limit","Limit", nmbin, 0., NMASSX);
  for ( Int_t i = 0 ; i < NMASSX ; i++ ) {
      TString nchi("hchi_");
      Int_t j = RhoMass[i]; 
      Int_t j1 = PiMass[i];
      nchi += (Int_t)j;
      nchi +="_";
      nchi += (Int_t)j1;
      TString nlh("hlh_");
      nlh += (Int_t)j;
      nlh += "_"; 
      nlh += (Int_t)j1; 

      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++ ) {
    Double_t xoff(-99999.); // to make the likelihood calculatable exp(x-xoff)
    if(NoSYS>=100&&im!=(NoSYS-100))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 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 
      
      // 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];
	  }
	}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.); 
      }

      // get cross section Ratio from 0 to nXsecR*step  
      Double_t tXsec = getExponential(gRand, nXsecRx, step, nXsecRx*step); //.Uniform()*nXsecR*step; 
      if(Iter<nXsecRx)tXsec = step*(Iter+0.5); // for computer xoff for scan the xsec   
      Int_t tBin = (Int_t)(tXsec/step);
      if(NoSYS==3&&Iter>=nXsecRx){
	tXsec =(Iter+0.5-nXsecRx)*step; 
	tBin = Iter-nXsecRx; 
      }
      if(Iter>=nXsecRx)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]){
            //*************************************************** 
	    // for special treat of 2-d Qvalue vs Mjj 
	    //***************************************************
	    // eff 
	    Double_t acc(1.);
	    Double_t acn(1.0);
	    for( Int_t i=0; i<NACC; ++i) {
	      acc *=geff[i][ich]; 
	      acn *=geff[i][ich+3];
	    }
	    if(im==3)heff[ich][0]->Fill(acc); 
	    float xb[NBKGH]={0.}; // background normalization
	    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] *=geffb[0][ich]; //common top lum 
	      if(j==4) xb[j] *=geffb[0][ich]; // common diboson lum
	    }
	    // save the background uncertainties 
	    if(im==3){
	      float sumb(0.); 
	      float sumbv(0.);
	      for( Int_t j=0; j<NBKGH-1; ++j){ 
		if(xb[j]>0.){
		  sumb +=XBKG[j][im][ich]->Integral(); 
		  sumbv +=xb[j]*XBKG[j][im][ich]->Integral(); 
		}
	      }
	      if(sumb>0.){
		heff[ich][1]->Fill(sumbv/sumb); 
		heff[ich][2]->Fill(sumbv);
	      }
	    }
	    
	    // No SYS applied here 
	    if(NoSYS==3||NoSYS==1){
	      acc =1.; 
	      acn =1.;
	    }
	    if(NoSYS==3||NoSYS==2){
	      for(Int_t j =0; j<NBKGH-1; ++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* xh2dn =(TH2F*)XHIGGS[im][ich+3];
	    TH2F* xb2d[NBKGH]; 
	    Double_t totb(0.); 
	    for(Int_t i=0; i<NBKGH-1; ++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)+acn*xh2dn->GetBinContent(k,k1) ); 
		efftot +=(acc*xh2d->GetBinContent(k,k1)+acn*xh2dn->GetBinContent(k,k1));

		//background 
		for( Int_t j=0; j<NBKGH-1; ++j){ 
		  if(xb[j]>0.){
		    float a1 = xb2d[j]->GetBinContent(k,k1); 
		    mc +=xb[j]*a1; // xb[j]==0, the hist is not defined
		    totb +=a1;
		  }
		}
		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 
	    } //
	    if(Iter==0||Iter==nXsecRx)std::cout<<" I am here "<< Iter<<" im "<< im <<" ich "<< ich <<" loglhS "<< loglhS<< " Data "<<xd2d->Integral()<<" Signal "<<xh2d->Integral()<<" Tot Bkg "<< totb<<std::endl;
	} // channels in consideration 
      } // loop all the channels
      if(Iter<nXsecRx){ // scan xsec to pick the best likelihood
	if(xoff<loglhS) xoff= loglhS; // pick the first loglhS as a reference point for likelihood calculation 
	if(im==3&&tBin<15)std::cout<<" I am here tBin "<< tBin<<" efftot "<< efftot <<" loglhS "<<loglhS<<" lh "<<lh[tBin]<<" dlh "<<efftot*TMath::Exp(loglhS-xoff)<<" n try "<< nc[tBin]<<" xoff "<<xoff<<std::endl; 
      }
      else{
	lh[tBin] +=efftot*TMath::Exp(loglhS-xoff);
      }
    }   // 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+ (0.95-A+avelh[Xsec])/avelh[Xsec]*step;
	break;
      }
    }

    std::cout<<" Rho Mass="<< RhoMass[im]
             <<" Pi Mass="<< PiMass[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(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 Run II Preliminary (1.9 fb^{-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 +="1";
      }
      else{
	mytxt +="0";
      }
    }
  }


  TString psFile("Technicolor_limit_cdf_2fb");
  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==20)psFile +="_uncorr"; // treat most of systematic uncorrelated, except lum, btag, top ew xsec   

  if(NoSYS>=100){
    psFile +="_mass";
    psFile +=(Int_t)RhoMass[NoSYS-100]; 
    psFile +="_";
    psFile +=(Int_t)PiMass[NoSYS-100];
  }
  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 Mass1 = RhoMass[i];
      Int_t Mass2 = PiMass[i];
      sprintf(cLimit,"Rho:%3i Pi:%3i GeV/c^{2} Limit = %3.1f pb",Mass1,Mass2,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();
    cx.Clear(); 
    gStyle->SetOptStat();
    gStyle->SetOptTitle();
    for(Int_t i=0; i<NCHA; ++i){ 
      if(icha[i]>0){ 
	if(heff[i][0]->Integral()>0.&&heff[i][1]->Integral()>0.){ 
	  cx.Divide(1,3); 
	  cx.cd(1); 
	  heff[i][0]->GetXaxis()->SetTitle(heff[i][0]->GetTitle());
	  heff[i][0]->Draw(); 
	  cx.cd(2);
	  heff[i][1]->GetXaxis()->SetTitle(heff[i][1]->GetTitle());
	  heff[i][1]->Draw(); 
	  cx.cd(3);
	  heff[i][2]->GetXaxis()->SetTitle(heff[i][2]->GetTitle());
	  heff[i][2]->Draw(); 
	  cx.Update(); 
	  cx.Clear();
	}
      }
    }
    //
  ps.Close();

  // output data file
  TString dataFile("Technicolor_limit_cdf_2fb");
  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==20)dataFile +="_uncorr"; // treat most of systematic uncorrelated, except lum, btag, top ew xsec   

  if(NoSYS>=100){
    dataFile +="_mass";
    dataFile +=(Int_t)RhoMass[NoSYS-100]; 
    dataFile +="_";
    dataFile +=(Int_t)PiMass[NoSYS-100];
  }

  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 <<" RhoMass="<< RhoMass[i]
         <<" PiMass="<<PiMass[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("Technicolor_limit_cdf_2fb");
  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==20)rootFile +="_uncorr"; // treat most of systematic uncorrelated, except lum, btag, top ew xsec   
 
  if(NoSYS>=100){
    rootFile +="_mass";
    rootFile +=(Int_t)RhoMass[NoSYS-100]; 
    rootFile +="_";
    rootFile +=(Int_t)PiMass[NoSYS-100];
  }
  rootFile += mytxt; 
  rootFile += ".root";
  TFile froot(rootFile,"recreate");
  for ( Int_t i = 0 ; i < NMASSX ; i++ ) {
      hchi[i]->Write();
      hlh [i]->Write();
  }
  hl->Write(); 

  for(Int_t i=0; i<NCHA; ++i){ 
    if(icha[i]>0){ 
      if(heff[i][0]->Integral()>0.&&heff[i][1]->Integral()>0.){ 
	heff[i][0]->Write();
	heff[i][1]->Write();
	heff[i][2]->Write();
      }
    }
  }

  froot.Close();

  // delete histograms
  for(Int_t i=0; i<NCHA; ++i){ 
    if(icha[i]>0){ 
      if(heff[i][0]->Integral()>0.&&heff[i][1]->Integral()>0.){ 
	heff[i][0]->Delete();
	heff[i][1]->Delete();
	heff[i][2]->Delete();
      }
    }
  }
  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 technicolor_limit_pseudoexperiment_cdf(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);
  //
  // dt1=st-st; dt2=st-jpb
  TString txt[NCHA] ={"lnubb_dt1c","lnubb_dt2c","lnubb_stnnc","lnubb_dt1n","lnubb_dt2n","lnubb_stnnn"};

  Double_t nfact(0);
  Double_t klum(1.0); // luminosity rescaling from the current 1 fb^-1 
  if(NoSYS>30&&NoSYS<40)klum = NoSYS-30;
  lgm[0] = 0.; 
  for(Int_t i =1; i<10000; ++i){ 
    nfact +=TMath::Log(i); 
    lgm[i] = nfact; 
  }
  char junk[120] =" ";
  FILE * infile = fopen("technicolor_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,120, 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,120, 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; 
    if(NoSYS==20)sbkg[j][0][i] = -1.*fabs(sbkg[j][0][i]); 
  }
  fgets(junk,120, infile);
  }
  
 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]; 
      }
    }
  }

  //get cdf data wh->lnubb with nnoutput
   TString Chanx[3]={"STST","STJP","1tagNN"};
   TString Chans[NMASS]={"1a","2b","3b","3c","4c","4d","5c","5d","5e"};
   for(int i =0; i<3; ++i){ 
      for(int j = 0; j<Nmax[i]; ++j){ 
	TFile* f = new TFile(Chanx[i]+"_ALL.root","read");
	TString Higgsx = "h_Mjj_vs_Qval_"+Chanx[i]; 
	XHIGGS[j][i] =(TH1F*)f->Get(Higgsx+"_cexo"+Chans[j]+"_2jet");
	XDATA[j][i] = (TH1F*)f->Get(Higgsx+"_Data_2jet"); 
	XBKG[0][j][i] = (TH1F*)f->Get(Higgsx+"_W+npAll_2jet"); 
	XBKG[1][j][i] = (TH1F*)f->Get(Higgsx+"_AE_2jet"); 
	XBKG[2][j][i] = (TH1F*)f->Get(Higgsx+"_WbbAll_2jet");
	XBKG[2][j][i]->Add((TH1F*)f->Get(Higgsx+"_WccAll_2jet"));
	XBKG[3][j][i] = (TH1F*)f->Get(Higgsx+"_ttop75_2jet");
	XBKG[3][j][i]->Add((TH1F*)f->Get(Higgsx+"_stop-s_2jet")); 
	XBKG[3][j][i]->Add((TH1F*)f->Get(Higgsx+"_stop-t_2jet")); 
	XBKG[4][j][i] = (TH1F*)f->Get(Higgsx+"_WW_2jet"); 
	XBKG[4][j][i]->Add((TH1F*)f->Get(Higgsx+"_WZ_2jet"));
	XBKG[4][j][i]->Add((TH1F*)f->Get(Higgsx+"_ZZ_2jet"));
	XBKG[4][j][i]->Add((TH1F*)f->Get(Higgsx+"_ZtauAll_2jet"));
      }
   }

   for(int i =0; i<3; ++i){ 
      for(int j = 0; j<Nmax[i]; ++j){ 
	TFile* f = new TFile(Chanx[i]+"_ALL.root","read");
	TString Higgsx = "h_Mjj_vs_Qval_"+Chanx[i]; 
	XHIGGS[j][i+3] =(TH1F*)f->Get(Higgsx+"_nexo"+Chans[j]+"_2jet");
	XDATA[j][i+3] = (TH1F*)f->Get(Higgsx+"_Data_2jet"); 
	XBKG[0][j][i+3] = (TH1F*)f->Get(Higgsx+"_W+npAll_2jet"); 
	XBKG[1][j][i+3] = (TH1F*)f->Get(Higgsx+"_AE_2jet"); 
	XBKG[2][j][i+3] = (TH1F*)f->Get(Higgsx+"_WbbAll_2jet");
	XBKG[2][j][i+3]->Add((TH1F*)f->Get(Higgsx+"_WccAll_2jet"));
	XBKG[3][j][i+3] = (TH1F*)f->Get(Higgsx+"_ttop75_2jet");
	XBKG[3][j][i+3]->Add((TH1F*)f->Get(Higgsx+"_stop-s_2jet")); 
	XBKG[3][j][i+3]->Add((TH1F*)f->Get(Higgsx+"_stop-t_2jet")); 
	XBKG[4][j][i+3] = (TH1F*)f->Get(Higgsx+"_WW_2jet"); 
	XBKG[4][j][i+3]->Add((TH1F*)f->Get(Higgsx+"_WZ_2jet"));
	XBKG[4][j][i+3]->Add((TH1F*)f->Get(Higgsx+"_ZZ_2jet"));
	XBKG[4][j][i+3]->Add((TH1F*)f->Get(Higgsx+"_ZtauAll_2jet"));
      }
   }


   std::cout<<" finishing reading WH "<<std::endl; 


  // Rescaling the luminosity for future sensitivity studies  
  for(int i=0; i<NCHA; ++i){ 
    for(int j=0; j<Nmax[i]; ++j){ 
      if(XDATA[j][i])XDATA[j][i]->Scale(klum);
      if(XHIGGS[j][i])XHIGGS[j][i]->Scale(klum);
	for(int k=0; k<NBKGH; ++k){ 
	  if(XBKG[k][j][i]){
	    XBKG[k][j][i]->Scale(klum);
	  }
	}
    }
  }

  const Int_t nXsecR = 50; // 300
  static Int_t nXsecRx = nXsecR;
  static Int_t nIter = 5000; //10k,  20000,100k or 500k 
  const float step =0.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]; 
    }
  }

  if(NoSYS==3) nIter = 2*nXsecRx; 
  // make copy of data: (starting 120 GeV) 
  for(Int_t i =0; i<NCHA; ++i){ 
    Int_t ix(0); 
    for(int j =0; j<Nmax[i]; ++j){ 
      if(XDATA[j][i]){
	ix = j;
	break; 
      }
    }
    if(XDATA[ix][i]){ 
      if(XDATA[ix][i]->GetDimension()<2){ 
	XMC[i] = (TH1F*)XDATA[ix][i]->Clone(); 
      }
      else{ 
	TH2F* x2d = (TH2F*)XDATA[ix][i];
	XMC[i] = (TH1F*)x2d->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_t j = RhoMass[i];
      Int_t j1 = PiMass[i];
      nchi += (Int_t)j;
      nchi +="_";
      nchi +=(Int_t)j1;

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

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

      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 +=(Int_t)j1; 
	ne +="_";
	ne +=txt[ix];
	hevt[ix][i] = new TH1F(ne,ne, 1000,0., 2000.*klum);
      }
  }

  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++ ) {
    Double_t xoff(-99999.); // to make the likelihood calculatable exp(x-xoff)
    if(NoSYS>=100&&im!=(NoSYS-100))continue;
    hchi[im]->Reset(); 
    hlh[im]->Reset(); 
    // make data: 
      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 

      // modeling efficiencies 

      // 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];
	  }
	}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<NCHA; ++ich){ 
	if(icha[ich]>0&&im<Nmax[ich]){

	  XMC[ich]->Reset();

	  float xb[NBKGH]={0.}; // background normalization
	  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] *=geffb[0][ich]; //common top luminosity 
	    if(j==4) xb[j] *=geffb[0][ich]; // common diboson luminosity
	  }

	  // No SYS applied here 
	  if(NoSYS==3||NoSYS==2){
	    for(Int_t j =0; j<NBKGH-1; ++j){
	      if(xb[j]>0)xb[j] = 1.;
	    }
	  }
	  // convert 1-d to 2d; 
	  TH2F* xd2d =(TH2F*)XMC[ich]; 
	  TH2F* xb2d[NBKGH]; 
	  for(Int_t i=0; i<NBKGH-1; ++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.; 
	      //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 = 0; 
	      if(mc>0)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());
	} // 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 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 

      // modeling efficiencies 

      // 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];
	  }
	}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.); 
      }


      // get cross section Ratio from 0 to nXsecR*step  
      Double_t tXsec = getExponential(gRand, nXsecRx, step, nXsecRx*step); //.Uniform()*nXsecR*step; 
      if(Iter<nXsecRx)tXsec = step*(0.5+Iter); // computer xoff 
      Int_t tBin = (Int_t)(tXsec/step);
      if(NoSYS==3&&Iter>=nXsecRx) {
	tXsec =(Iter+0.5-nXsecRx)*step;
	tBin = Iter-nXsecRx; 
      }
      if(Iter>=nXsecRx)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.);
	  Double_t acn(1.);
	  for( Int_t i=0; i<NACC; ++i) {
	    acc *=geff[i][ich];
	    acn *=geff[i][ich+3];
	  }
	  float xb[NBKGH]={0.}; // background normalization
	  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] *=geffb[0][ich]; //common top lum 
	    if(j==4) xb[j] *=geffb[0][ich]; // common diboson lum 
	  }

	  // No SYS applied here 
	  if(NoSYS==3||NoSYS==1){
	    acc = 1.; 
	    acn = 1.;
	  }
	  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* xh2dn =(TH2F*)XHIGGS[im][ich+3];
	  TH2F* xb2d[NBKGH]; 
	  for(Int_t i=0; i<NBKGH-1; ++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)+acn*xh2dn->GetBinContent(k,k1)); 
	      efftot +=(acc*xh2d->GetBinContent(k,k1)+acn*xh2dn->GetBinContent(k,k1));
	      
	      //background 
	      for( Int_t j=0; j<NBKGH-1; ++j){ 
		  if(xb[j]>0){
		    float a1 = xb2d[j]->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 
	  } //
	} // channels in consideration 	
      } // loop all the channels
      if(Iter<nXsecRx){  // scan the xsec first to pick the best likelihood 
	if(xoff<loglhS)xoff = loglhS; // pick the first loglhS as a reference point for likelihood calculation
      }
      else{
	lh[tBin] +=efftot*TMath::Exp(loglhS-xoff);
	if(Iter<(3+nXsecRx)&&ie%10==0)std::cout<<" Iter "<< Iter<<" tBin "<<tBin<<" efftot "<< efftot <<" loglhS "<<loglhS<<" lh "<<lh[tBin]<<" dlh "<<efftot*TMath::Exp(loglhS-xoff)<<" n try "<< nc[tBin]<<" xoff "<< xoff<<std::endl;
      }
    }   // 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+(0.95-A+avelh[Xsec])/avelh[Xsec]*step;
	break;
      }
    }

    if(ie%10==0)std::cout <<" Pseudo experiments: "<<ie
			   <<" Rho Mass="<< RhoMass[im]
			  <<" Pi Mass="<<PiMass[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==NCHA){ 
    mytxt += "all"; }
  else{
    for(Int_t i=0; i<NCHA; ++i){ 
      if(icha[i]>0){
	mytxt +="1";
      }
      else{
	mytxt +="0";
      }
    }
  }
  TString psFile("Technicolor_limit_pseudoexperiments_cdf_2fb");
  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==20)psFile +="_uncorr"; // treat most of systematic uncorrelated, except lum, btag, top ew xsec   
  if(NoSYS>30&&NoSYS<40){
    psFile +="_Lum_"; 
    Int_t i0 = NoSYS-30; 
    psFile +=(Int_t)i0; 
    psFile +="invfb"; 
  }
  if(NoSYS>=100){
    psFile +="_mass";
    psFile +=(Int_t)RhoMass[NoSYS-100];
    psFile +="_";
    psFile +=(Int_t)PiMass[NoSYS-100];
  }

  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("Technicolor_limit_pseudoexperiments_cdf_2fb");
  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==20)dataFile +="_uncorr"; // treat most of systematic uncorrelated, except lum, btag, top ew xsec   

  if(NoSYS>30&&NoSYS<40){
    dataFile +="_Lum_"; 
    Int_t i0 = NoSYS-30; 
    dataFile +=(Int_t)i0; 
    dataFile +="invfb"; 
  }

  if(NoSYS>=100){
    dataFile +="_mass";
    dataFile +=(Int_t)RhoMass[NoSYS-100];
    dataFile +="_";
    dataFile +=(Int_t)PiMass[NoSYS-100];
  }

  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 <<" Rho Mass="<< RhoMass[i]
         <<" Pi Mass="<<PiMass[i]
	 <<" Mean of 95% CL="<< hl[i]->GetMean(1) // 1=x axis
	 <<" RMS "<< hl[i]->GetRMS(1)<<" Median "<<GetMedian(hl[i])
	 << std::endl;
    for( Int_t ix=0; ix<NCHA; ix++){ 
      if(hevt[ix][i]->Integral()>0){
	fout<<" Rho mass "<< RhoMass[i]<<" Pi Mass "<<PiMass[i]<<" "<<txt[ix]<<" events="<<hevt[ix][i]->GetMean(1)
	    <<" RMS="<<hevt[ix][i]->GetRMS(1)<<std::endl;
      }
    }
  }
  fout.close();

  // output root file
  TString rootFile("Technicolor_limit_pseudoexperiments_cdf_2fb");
  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==20)rootFile +="_uncorr"; // treat most of systematic uncorrelated, except lum, btag, top ew xsec  

  if(NoSYS>30&&NoSYS<40){
    rootFile +="_Lum_"; 
    Int_t i0 = NoSYS-30; 
    rootFile +=(Int_t)i0; 
    rootFile +="invfb"; 
  } 
  if(NoSYS>=100){
    rootFile +="_mass";
    rootFile +=(Int_t)RhoMass[NoSYS-100];
    rootFile +="_";
    rootFile +=(Int_t)PiMass[NoSYS-100];
  }

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


double GetMedian(TH1F* th1){ 
  double Median;
  int Nentry = th1->Integral();
  int Nhalf = Nentry/2;
 
  int Nbin = th1->GetNbinsX();
  int Xmin = th1->GetXaxis()->GetXmin();
  int Xmax = th1->GetYaxis()->GetXmax();
  int Nacc = 0;
  int CenterBin;
  for(int i=0;i<Nbin;i++){
    Nacc += th1->GetBinContent(i+1);
    if(Nacc >= Nhalf){
      CenterBin = i+1;
      break;
    }
  }
  Double_t nx = th1->GetBinContent(CenterBin);
  Double_t x = (Nhalf - Nacc + nx);
  Double_t bw = th1->GetBinWidth(CenterBin);
  Median = th1->GetBinLowEdge(CenterBin) + x/nx*bw;
  return Median;
}

int main(int argc, char* argv[]) { // two arguments limit/pseudo,SYS

  timeval a,b;
  gettimeofday(&a,NULL);
  std::cout<< argv[1]<<" "<<argv[2]<<std::endl;
  Int_t do_pseudoexperiment=atoi(argv[1]);
  Int_t NoSys=atoi(argv[2]);
  std::cout<< " do pseudoexperiment "<<do_pseudoexperiment<<" NoSys "<<NoSys<<std::endl;
  if(do_pseudoexperiment){
    technicolor_limit_pseudoexperiment_cdf(NoSys);
  }
  else{
    technicolor_limit_cdf(NoSys);
  }
  gettimeofday(&b,NULL);
  double deltat=a.tv_sec*1000.0+a.tv_usec/1000.0;
  deltat=b.tv_sec*1000.0+b.tv_usec/1000.0-deltat;
  std::cout<< argv[1]<<" "<<argv[2]<<std::endl;
  printf(" %f sec run time\n",deltat/1000);
  printf("\n");
  return 0;

}

#endif