#include "TMath.h"
#include "TH2.h"
#include "TF2.h"
#include "TGraph.h"
#include "TCanvas.h"
#include "TLatex.h"
#include "TMarker.h"
#include "TBackCompFitter.h"

double bigaus_function(double *x, double *p) {
  double u = (x[0]-p[1])/p[2];
  double v = (x[1]-p[3])/p[4];
  double rho = p[5];
  double c = 1. - rho*rho;
  double z = u*u - 2.*rho*u*v + v*v;
  return  p[0]/(2. * TMath::Pi() * p[2] * p[4] * std::sqrt(c) )* std::exp(- z / (2. * c) );
}

void exampleBigaus() {
  int n = 100000;
  int nbins = 50;
  double xmin = -10;
  double xmax = 10;
  TH2 * h2 = new TH2D("h2", "h2", nbins, xmin, xmax, nbins, xmin, xmax);
#if ROOT_VERSION_CODE >= ROOT_VERSION(6,7,6)
  // for ROOT 6.07.06 and newer
  TF2 *f2 = new TF2("bigaus", "bigaus"); // built-in function
#else
  // for older ROOT version
  TF2 *f2 = new TF2("bigaus", bigaus_function, xmin, xmax, xmin, xmax, 6);
#endif
  
  f2->SetParameters(1,2,2,-1,3,0.3);
  
  h2->FillRandom("bigaus");
  
  f2->SetParameters(100,2,2,-1,3,0.3);
  
  h2->Fit("bigaus");
  
  h2->Draw("COLZ");
  
  // draw contours of fit result
  
  TBackCompFitter *fitter =
    ((TBackCompFitter *)(TVirtualFitter::GetFitter()));
  
  int par1 = 1; // index of meanx
  int par2 = 3; // index of meany
  
  TGraph *gr1 = new TGraph( 80 );
  gr1->SetTitle(";top mass [GeV];#alpha_{S}^{}");
  gr1->SetFillColor(kGreen); gr1->SetLineColor(kBlack);
  double cl1 = 0.6827; // 1 sigma
  fitter->Contour(par1, par2, gr1, cl1);
  gr1->SetPoint(gr1->GetN(), gr1->GetX()[0], gr1->GetY()[0]); // "close" it
  
  TGraph *gr2 = new TGraph( 80 );
  gr2->SetTitle(";top mass [GeV];#alpha_{S}^{}");
  gr2->SetFillColor(kYellow); gr2->SetLineColor(kBlack);
  double cl2 = 0.9545; // 2 sigma
  fitter->Contour(par1, par2, gr2, cl2);
  gr2->SetPoint(gr2->GetN(), gr2->GetX()[0], gr2->GetY()[0]); // "close" it
  
  TMarker *m0 = new TMarker(fitter->GetParameter(par1),
                            fitter->GetParameter(par2),
                            kFullTriangleDown);
  m0->SetMarkerColor(kBlue);
  m0->SetMarkerSize(1.5);
  
  TLatex *t0 = new TLatex(fitter->GetParameter(par1),
                          fitter->GetParameter(par2),
                          TString::Format("   [%.5g GeV; %.5g]",
                                          fitter->GetParameter(par1),
                                          fitter->GetParameter(par2)));
  t0->SetTextSize(0.03333);
  
  TLatex *t1 = new TLatex(2.005, -0.980, "1 #sigma"); // "manual" position
  t1->SetTextSize(0.03333);
  
  TLatex *t2 = new TLatex(2.030, -0.940, "2 #sigma"); // "manual" position
  t2->SetTextSize(0.03333);
  
  TCanvas *c1 = new TCanvas("c1", "c1");
  gr2->Draw("AF"); gr2->Draw("L");
  gr1->Draw("F"); gr1->Draw("L");
  t0->Draw();
  t1->Draw();
  t2->Draw();
  m0->Draw();
  c1->Update();
}