/***************************************************************************** 
 * Project: RooFit                                                           * 
 *                                                                           * 
 * This code was autogenerated by RooClassFactory                            * 
 *****************************************************************************/ 

// Your description goes here... 

#include "Riostream.h" 

#include "Roo2DimGauss.h" 
#include "RooAbsReal.h" 
#include "RooAbsCategory.h" 
#include <math.h> 
#include "TMath.h" 

ClassImp(Roo2DimGauss) 

Roo2DimGauss::Roo2DimGauss(const char *name, const char *title, 
			   RooAbsReal& _x,
			   RooAbsReal& _y,
			   RooAbsReal& _rho,
			   RooAbsReal& _sigX,
			   RooAbsReal& _sigY,
			   RooAbsReal& _muX,
			   RooAbsReal& _muY) :
RooAbsPdf(name,title), 
  x("x","x",this,_x),
  y("y","y",this,_y),
  rho("rho","rho",this,_rho),
  sigX("sigX","sigX",this,_sigX),
  sigY("sigY","sigY",this,_sigY),
  muX("muX","muX",this,_muX),
  muY("muY","muY",this,_muY)
{ 
} 


Roo2DimGauss::Roo2DimGauss(const Roo2DimGauss& other, const char* name) :  
RooAbsPdf(other,name), 
  x("x",this,other.x),
  y("y",this,other.y),
  rho("rho",this,other.rho),
  sigX("sigX",this,other.sigX),
  sigY("sigY",this,other.sigY),
  muX("muX",this,other.muX),
  muY("muY",this,other.muY)
{ 
} 


//Int_t RooG2DimGauss::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
 /*
{
  if (matchArgs(allVars,analVars,x)) return 1 ;
  return 0 ;
}


Double_t Roo2DimGauss::analyticalIntegral(Int_t code, const char* rangeName) const 
{
  assert(code==1) ;

}
 */

Double_t Roo2DimGauss::evaluate() const 
{ 
  double factor=2*TMath::Pi()*sigX*sigY;
  factor*=TMath::Sqrt(1-rho*rho);
  factor=1./factor;

  double val1=(x-muX)*(x-muX);
  val1/=sigX*sigX;;
  double val2=(y-muY)*(y-muY);
  val2/=sigY*sigY;
  double val3=2*rho*(x-muX)*(y-muY);
  val3/=sigX*sigY;

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

  double retval=factor*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<<" y:"<<y<<" rho:"<<rho<<" sigx"<<sigX<<" sigy:"<<sigY<<" muX:"<<muX<<" muY"<<muY<<"\n";
    }
  */
  return factor*exp(exponent);
}