#include "TMath.h"
#include <iostream>

Double_t gauss2d(Double_t* x, Double_t* par)
{ 
  //par=[sigX,sigY,rho,meanX,meanY,amp]
  double factor=2*TMath::Pi()*par[0]*par[1];
  factor*=TMath::Sqrt(1-par[2]*par[2]);
  factor=1./factor;

  double val1=(x[0]-par[3])*(x[0]-par[3]);
  val1/=par[0]*par[0];;
  double val2=(x[1]-par[4])*(x[1]-par[4]);
  val2/=par[1]*par[1];
  double val3=2*par[2]*(x[0]-par[3])*(x[1]-par[4]);
  val3/=par[0]*par[1];

  double exponent=-1./(2*(1-par[2]*par[2]));
  exponent*=(val1+val2-val3);

  //double retval=factor*TMath::Exp(exponent);
  double retval=par[5]*TMath::Exp(exponent);

  /*
  if (retval!=retval)
    std::cout<<"not a number\n";
  if (retval<=0)
    {
      std::cout<<"less than 0!!! "<<retval<<"\n";
  
    }
  std::cout<<"x:"<<x[0]<<" y:"<<x[1]<<" par[2]:"<<par[2]<<" sigx"<<par[0]<<" sigy:"<<par[1]<<" par[3]:"<<par[3]<<" par[4]"<<par[4]<<"\n";
  */
  return retval;
}