const Int_t MAXBINS = 200;
const Double_t sigmas=0.88; // nlo; 0.6 at lo?
const Double_t sigmatd=21.0;
const Double_t sigmats=5.4;
const Double_t sigmatb=1.96;

void vti(
	 Int_t nbins=100,
	 bool useTotalSingleTop=true,
	 Double_t beta=1.
	 )
{
  histos(nbins,"d0",useTotalSingleTop,0.79,4.7,1.3,1.3,1.0,0.9,0.9,4.2,1.8,1.4,beta); // D0
  histos(nbins,"cdf",useTotalSingleTop,0.61,2.2,0.7,0.7,1.6,0.9,0.8,0.8,0.7,0.6,beta); // CDF
}

void histos(
	    Int_t nbins=10,
	    TString name="d0",
	    bool useTotalSingleTop=true,
	    Double_t r_min=0.79,
	    Double_t expst=4.7,
	    Double_t expst_err_up=1.3,
	    Double_t expst_err_down=1.3,
	    Double_t exps=1.0,
	    Double_t exps_err_up=0.9,
	    Double_t exps_err_down=0.9,
	    Double_t expt=4.2,
	    Double_t expt_err_up=1.8,
	    Double_t expt_err_down=1.4,
	    Double_t beta=1., // probably 0.5 is more realistic
	    Double_t ab=1.,
	    Double_t ad=1.,
	    Double_t as=1.
	    )
{
  //  gROOT->Reset();
  gStyle->SetOptStat(0000);
  gStyle->SetCanvasBorderMode(-1);
  gStyle->SetCanvasBorderSize(1);
  gStyle->SetCanvasColor(10);

  TH3F* ckm = new TH3F("histo","allowed by R, s-channel, t-channel", nbins,0.,1.,nbins,0.,1.,nbins,0.,1.);
  TH3F* ckm_r = new TH3F("histo_r","allowed by R", nbins,0.,1.,nbins,0.,1.,nbins,0.,1.);
  TH3F* ckm_st = new TH3F("histo_st","allowed by single top", nbins,0.,1.,nbins,0.,1.,nbins,0.,1.);
  TH3F* ckm_s = new TH3F("histo_s","allowed by s-channel", nbins,0.,1.,nbins,0.,1.,nbins,0.,1.);
  TH3F* ckm_t = new TH3F("histo_t","allowed by t-channel", nbins,0.,1.,nbins,0.,1.,nbins,0.,1.);

  TH2F* xy = new TH2F("xy","allowed by R, s-channel, t-channel", nbins,0.,1.,nbins,0.,1.);
  TH2F* xy_r = new TH2F("xy_r","allowed by R", nbins,0.,1.,nbins,0.,1.);
  TH2F* xy_st = new TH2F("xy_st","allowed by single top", nbins,0.,1.,nbins,0.,1.);
  TH2F* xy_s = new TH2F("xy_s","allowed by s-channel", nbins,0.,1.,nbins,0.,1.);
  TH2F* xy_t = new TH2F("xy_t","allowed by t-channel", nbins,0.,1.,nbins,0.,1.);

  TH2F* zx = new TH2F("zx","allowed by R, s-channel, t-channel", nbins,0.,1.,nbins,0.,1.);
  TH2F* zx_r = new TH2F("zx_r","allowed by R", nbins,0.,1.,nbins,0.,1.);
  TH2F* zx_st = new TH2F("zx_st","allowed by single top", nbins,0.,1.,nbins,0.,1.);
  TH2F* zx_s = new TH2F("zx_s","allowed by s-channel", nbins,0.,1.,nbins,0.,1.);
  TH2F* zx_t = new TH2F("zx_t","allowed by t-channel", nbins,0.,1.,nbins,0.,1.);

  TH2F* zy = new TH2F("zy","allowed by R, s-channel, t-channel", nbins,0.,1.,nbins,0.,1.);
  TH2F* zy_r = new TH2F("zy_r","allowed by R", nbins,0.,1.,nbins,0.,1.);
  TH2F* zy_st = new TH2F("zy_st","allowed by single top", nbins,0.,1.,nbins,0.,1.);
  TH2F* zy_s = new TH2F("zy_s","allowed by s-channel", nbins,0.,1.,nbins,0.,1.);
  TH2F* zy_t = new TH2F("zy_t","allowed by t-channel", nbins,0.,1.,nbins,0.,1.);

  Float_t prob[MAXBINS][MAXBINS][MAXBINS];
  Float_t prob_r[MAXBINS][MAXBINS][MAXBINS];
  Float_t prob_st[MAXBINS][MAXBINS][MAXBINS];
  Float_t prob_s[MAXBINS][MAXBINS][MAXBINS];
  Float_t prob_t[MAXBINS][MAXBINS][MAXBINS];
  Float_t prob_trivial[MAXBINS][MAXBINS][MAXBINS];

  for (Int_t i=0; i<nbins; i++) {
    Float_t vtd = (Float_t)i/(Float_t)nbins;
    for (Int_t j=0; j<nbins; j++) {
      Float_t vts = (Float_t)j/(Float_t)nbins;
      for (Int_t k=0; k<nbins; k++) {
	Float_t vtb = (Float_t)k/(Float_t)nbins;

	// initialize
	prob[i][j][k]=0.;
	prob_r[i][j][k]=0.;
	prob_st[i][j][k]=0.;
	prob_s[i][j][k]=0.;
	prob_t[i][j][k]=0.;
	prob_trivial[i][j][k]=0.;

	// "trivial" constraint (radius <= 1)
	Float_t radius=sqrt(vtd**2+vts**2+vtb**2);
	if (radius < 1.) prob_trivial[i][j][k]=1.;
	// constraint from R
	Float_t r=0.;
	if (radius>0.) { // just to avoid a division by 0 in the first bin
	  r = vtb**2/radius**2;
	  if (r > r_min) prob_r[i][j][k]=1.;
	}
	// constraint from s-channel
	Float_t s = r*sigmas*vtb**2;
	if (s < exps+2*exps_err_up && s > exps-2*exps_err_down) prob_s[i][j][k]=1.;
	// constraint from t-channel
	Float_t t = r*(
			 ad*sigmatd*vtd**2 + as*sigmats*vts**2 + ab*sigmatb*vtb**2
			 + 2*beta*(vtd**2 + vts**2)*sigmas
			 );
	if (t < expt+2*expt_err_up && t > expt-2*expt_err_down) prob_t[i][j][k]=1.;
	// constraint from single top
	Float_t st = s+t;
	if (st < expst+2*expst_err_up && st > expst-2*expst_err_down) prob_st[i][j][k]=1.;

	// combined constraint
	if (useTotalSingleTop) {
	  prob[i][j][k]=prob_r[i][j][k]*prob_st[i][j][k]*prob_trivial[i][j][k];
	} else {
	  prob[i][j][k]=prob_r[i][j][k]*prob_s[i][j][k]*prob_t[i][j][k]*prob_trivial[i][j][k];
	}

	ckm->Fill(vtd,vts,vtb,prob[i][j][k]);
	ckm_r->Fill(vtd,vts,vtb,prob_r[i][j][k]);
	ckm_st->Fill(vtd,vts,vtb,prob_st[i][j][k]);
	ckm_s->Fill(vtd,vts,vtb,prob_s[i][j][k]);
	ckm_t->Fill(vtd,vts,vtb,prob_t[i][j][k]);

	xy->Fill(vtd,vts,prob[i][j][k]);
	xy_r->Fill(vtd,vts,prob_r[i][j][k]);
	xy_st->Fill(vtd,vts,prob_st[i][j][k]);
	xy_s->Fill(vtd,vts,prob_s[i][j][k]);
	xy_t->Fill(vtd,vts,prob_t[i][j][k]);

	zx->Fill(vtd,vtb,prob[i][j][k]);
	zx_r->Fill(vtd,vtb,prob_r[i][j][k]);
	zx_st->Fill(vtd,vtb,prob_st[i][j][k]);
	zx_s->Fill(vtd,vtb,prob_s[i][j][k]);
	zx_t->Fill(vtd,vtb,prob_t[i][j][k]);

	zy->Fill(vts,vtb,prob[i][j][k]);
	zy_r->Fill(vts,vtb,prob_r[i][j][k]);
	zy_st->Fill(vts,vtb,prob_st[i][j][k]);
	zy_s->Fill(vts,vtb,prob_s[i][j][k]);
	zy_t->Fill(vts,vtb,prob_t[i][j][k]);
      }
    }
  }


  TString canvasname1(name);
  TString canvasname2(name);
  canvasname2 += "_r";
  TString canvasname3(name);
  canvasname3 += "_st";
  TString canvasname4(name);
  canvasname4 += "_s";
  TString canvasname5(name);
  canvasname5 += "_t";

  plot(ckm,xy,zx,zy,canvasname1);
  plot(ckm_r,xy_r,zx_r,zy_r,canvasname2);
  plot(ckm_st,xy_st,zx_st,zy_st,canvasname3);
  plot(ckm_s,xy_s,zx_s,zy_s,canvasname4);
  plot(ckm_t,xy_t,zx_t,zy_t,canvasname5);

  delete ckm;
  delete ckm_r;
  delete ckm_st;
  delete ckm_s;
  delete ckm_t;
  delete xy;
  delete xy_r;
  delete xy_st;
  delete xy_s;
  delete xy_t;
  delete zx;
  delete zx_r;
  delete zx_st;
  delete zx_s;
  delete zx_t;
  delete zy;
  delete zy_r;
  delete zy_st;
  delete zy_s;
  delete zy_t;
}

void plot(TH3F* h, TH2F* h_xy, TH2F* h_zy, TH2F* h_zx, TString canvasname)
{
  // 2D projections
  //  TH2D* h_xy = h->Project3D("xy");
  //  TH2D* h_zy = h->Project3D("zy");
  //  TH2D* h_zx = h->Project3D("zx");
  
  // axes
  h->GetXaxis()->SetTitle("Vtd");
  h->GetYaxis()->SetTitle("Vts");
  h->GetZaxis()->SetTitle("Vtb");

  h_xy->GetXaxis()->SetTitle("Vtd");
  h_xy->GetYaxis()->SetTitle("Vts");

  h_zx->GetXaxis()->SetTitle("Vtd");
  h_zx->GetYaxis()->SetTitle("Vtb");

  h_zy->GetXaxis()->SetTitle("Vts");
  h_zy->GetYaxis()->SetTitle("Vtb");


  // plot
  TCanvas* myCanvas = new TCanvas(canvasname,canvasname);
  myCanvas->Divide(2,2);
  myCanvas->cd(1);
  h->DrawCopy();
  myCanvas->cd(2);
  h_xy->DrawCopy();
  myCanvas->cd(3);
  h_zy->DrawCopy();
  myCanvas->cd(4);
  h_zx->DrawCopy();
  TString gif(canvasname);
  gif += ".gif";
  myCanvas->Print(gif);

}