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 example2dGausContour() {

   // draw 1-sigma and 2 s-gma contour of a 2d correlated gaussian
   // Definition of 1 sigma contour is the contour which has a prob content of 0.708
   // which is normal_cdf(1,1)^2
   // 


   TCanvas *c  = new TCanvas("c","Contours",0,0,800,700);


#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, 0,1,0,1, 6);
#endif


   f2->SetFillStyle(1000);
   f2->SetLineWidth(1);
   f2->SetRange(-5,-5,5,5);
   f2->SetParameters(1,0,1,0,1,-0.7);
   f2->SetNpx(1000);
   f2->SetNpy(1000);
   double contours[2];
   // determine level of gaussian at 1 and 2 sigma using a non-correlated 2d gaussian
   // corr factor needed because the constant term is different in the case of a correlated 2d-gaussian

#if ROOT_VERSION_CODE >= ROOT_VERSION(6,7,6)
    double corr = f2->Eval(0,0)/ROOT::Math::bigaussian_pdf(0,0,1,1,0);
    contours[0] = corr*ROOT::Math::bigaussian_pdf(2,2,1,1,0); //0.05;
    contours[1] = corr*ROOT::Math::bigaussian_pdf(1,1,1,1,0);
#else
   double corr = f2->Eval(0,0)/(ROOT::Math::normal_pdf(0,1)*ROOT::Math::normal_pdf(0,1));
   contours[0] = corr*ROOT::Math::normal_pdf(2,1)*ROOT::Math::normal_pdf(2,1);
   contours[1] = corr*ROOT::Math::normal_pdf(1,1)*ROOT::Math::normal_pdf(1,1);
#endif


   f2->SetContour(2,contours);


   Int_t colors[2] = {kYellow, kGreen};
   gStyle->SetPalette(2,colors);

   f2->Draw("cont0 z");
   f2->Draw("cont3 same");

}