#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooStats/FeldmanCousins.h"
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooArgSet.h"
#include "RooWorkspace.h"
#include "RooPoisson.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TLatex.h"
#include "TMarker.h"
#include "TGraph2D.h"
#include "TObjArray.h"

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

void histInvSm(){

  string hmass = "100";
  string amass = "10";

  RooRealVar* mu = new RooRealVar("mu", "mu", 1., 0, 2);
  RooRealVar* BRt = new RooRealVar("BRt", "BRt", 1, 0, 1);
  RooRealVar* BRg = new RooRealVar("BRg", "BRg", 1, 0, 1);
  RooRealVar* BRc = new RooRealVar("BRc", "BRc", 1, 0, 1);

  RooDataSet* poiToTest = new RooDataSet("poiToTest","",RooArgSet(*BRt,*BRg,*BRc,*mu));
  RooDataSet* BRsToTest = new RooDataSet("BRsToTest","",RooArgSet(*BRt,*BRg,*BRc));
  double poi1 = 0.0;
  double poi2 = 0.0;
  double poi3 = 0.0;
  int steps = 11;
  int musteps = 11;
  double stepsize = 1.0/(steps-1);
  double mustepsize = 2.0/(musteps-1);

  int brPts  = 0;
  int totPts = 0;

  for(int kk=0; kk<steps; ++kk){
    poi1 = kk*stepsize;
    for(int nn=0; nn<steps; ++nn){
      poi2 = nn*stepsize;
      poi3 = 1. - poi1 - poi2;
      if(poi3 < 1.e-10) poi3 = 0.0;
      if(poi1+poi2+poi3<=1.0 && poi1+poi2+poi3>=0.0 && poi1>=0.0 && poi2>=0.0 && poi3>=0.0){
	BRg->setVal(poi1);
	BRt->setVal(poi2);
	BRc->setVal(poi3);
	BRsToTest->add(RooArgSet(*BRt,*BRg,*BRc));
	++brPts;
	for(int mm=0; mm<musteps; ++mm){
	  mu->setVal(mm*mustepsize);
	  poiToTest->add(RooArgSet(*BRt,*BRg,*BRc,*mu));
	  ++totPts;
	}
      }
    }
  }

  TH2F* altHist = new TH2F("altHist", "title", 11, -.05, 1.05, 11, -.05, 1.05);
  altHist->SetTitle("");
  altHist->SetStats(kFALSE);

  double muULs[66] = {1.2,1.2,1.0,1.0,0.6,0.4,0.4,0.4,0.4,0.2,0.2,1.6,1.4,1.0,0.8,0.6,0.4,0.4,0.4,0.2,0.2,2.0,1.6,1.2,0.8,0.6,0.4,0.4,0.4,0.2,2.0,1.6,1.2,0.8,0.6,0.6,0.4,0.2,2.0,1.6,1.2,0.8,0.6,0.4,0.4,2.0,1.8,1.4,0.8,0.6,0.4,2.0,1.4,1.0,0.8,0.6,1.8,1.4,1.0,0.6,1.4,1.2,0.8,1.2,1.0,1.0};

  int ww = 0;
  int numOnes = 0;

  for(int kk=0; kk<steps; ++kk){
    poi1 = kk*stepsize;
    for(int nn=0; nn<steps; ++nn){
      poi2 = nn*stepsize;
      poi3 = 1. - poi1 - poi2;
      if(poi3 < 1.e-10) poi3 = 0.0;
      if(poi1+poi2+poi3<=1.0 && poi1+poi2+poi3>=0.0 && poi1>=0.0 && poi2>=0.0 && poi3>=0.0){
        BRg->setVal(poi1);
        BRt->setVal(poi2);
        BRc->setVal(poi3);
	altHist->SetBinContent(kk+1,nn+1,muULs[ww]);
	if(muULs[ww]==1.0) ++numOnes;
	++ww;
      }
    }
  }

  double muX[numOnes];
  double muY[numOnes];
  double dummies[numOnes];
  int rrr = 0;
  int xx = 0;

  for(int kk=0; kk<steps; ++kk){
    poi1 = kk*stepsize;
    for(int nn=0; nn<steps; ++nn){
      poi2 = nn*stepsize;
      poi3 = 1. - poi1 - poi2;
      if(poi3 < 1.e-10) poi3 = 0.0;
      if(poi1+poi2+poi3<=1.0 && poi1+poi2+poi3>=0.0 && poi1>=0.0 && poi2>=0.0 && poi3>=0.0){
        BRg->setVal(poi1);
        BRt->setVal(poi2);
        BRc->setVal(poi3);
	if(muULs[xx]==1.0){
	  muX[rrr] = BRg->getVal();
	  muY[rrr] = BRt->getVal();
	  dummies[rrr] = 1.0;
	  ++rrr;
	}
	++xx;
      }
    }
  }

  //////////////////////////////
  // PLOTTING ISSUE STARTS HERE
  //////////////////////////////

  TGraph2D* histGr = new TGraph2D(altHist);

  gROOT->SetStyle("Plain");
  gStyle->SetPalette(1);

  TCanvas* tempCanv = new TCanvas("","",2);
  altHist->Draw("cont,list");
  gPad->Update();

  TObjArray* contours = (TObjArray*) gROOT->GetListOfSpecials()->FindObject("contours");
  int ncontours = contours->GetSize();
  cout << "ncontours is " << ncontours << endl;
  TList* list = (TList*) contours->At(10);
  TGraph* gr1 = (TGraph*) list->First();

  TCanvas* canv = new TCanvas("canv","canv",2);
  
  TPad* z1 = new TPad("z1","z1",0.0,0.0,1.0,1.0);
  z1->Draw();
  TPad* z2 = new TPad("z2","z2",0.0,0.0,1.0,1.0);
  z2->Draw();
  z2->SetFillStyle(4000);

  histGr->SetMaximum(2.0);

  z1->cd();
  histGr->Draw("cont4z");

  z2->cd();
  gPad->DrawFrame(0.0,0.0,1.0,1.0);
  gr1->Draw();


}