#include "RooFit.h"

#include "Riostream.h"
#include "Riostream.h"
#include <math.h>
#include "TMath.h"

#include "RooGaussVar.h"
#include "RooAbsReal.h"
#include "RooRealVar.h"
#include "RooArgList.h"
#include "RooMsgService.h"
#include "RooMath.h"
#include "TMath.h"

using namespace std;

ClassImp(RooGaussVar)
;


//_____________________________________________________________________________
RooGaussVar::RooGaussVar()
{
  // Default constructor
}


//_____________________________________________________________________________
RooGaussVar::RooGaussVar(const char* name, const char* title, 
			     RooAbsReal& _x, RooAbsReal& _height, RooAbsReal& _mean, RooAbsReal& _sigma) :
  RooAbsReal(name, title),
  x("x", "Dependent", this, _x),
  height("height", "Height", this,_height),
  mean("mean", "Mean", this, _mean),
  sigma("sigma", "Sigma", this, _sigma)

  {
;
  }



//_____________________________________________________________________________
/*
RooGaussVar::RooGaussVar(const char* name, const char* title,
                           RooAbsReal& x) :
  RooAbsReal(name, title),
  _x("x", "Dependent", this, x),
{
  // Constructor of gaussian function
  _parIter = _parList.createIterator() ;
}                                                                                                                                 
*/


//_____________________________________________________________________________

RooGaussVar::RooGaussVar(const RooGaussVar& other, const char* name) :
  RooAbsReal(other, name), 
  x("x", this, other.x),
  height("height", this, other.height),
  mean("mean", this, other.mean),
  sigma("sigma", this, other.sigma)
{

}




//_____________________________________________________________________________
RooGaussVar::~RooGaussVar() 
{
  // Destructor
}




//_____________________________________________________________________________
Double_t RooGaussVar::evaluate() const 
{
  // Calculate and return value of gaussian

  Double_t arg(0) ;


  arg=(x-mean.arg().getVal()) / sigma.arg().getVal();
  if (arg < -39.0 || arg > 39.0) return 0.0; 

  Double_t res = height.arg().getVal() * TMath::Exp(-0.5*arg*arg);

  //std::cout<<"gausvar evaluate:"<<res<<"\n";
  return res;
}



//_____________________________________________________________________________
Int_t RooGaussVar::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const 
{
  // Advertise that we can internally integrate over x

  if (matchArgs(allVars, analVars, x)) return 1;
  if (matchArgs(allVars, analVars, mean)) return 2 ;
  return 0;
}



//_____________________________________________________________________________
Double_t RooGaussVar::analyticalIntegral(Int_t code, const char* rangeName) const 
{

	assert(code==1 || code==2) ;

	static const Double_t root2 = sqrt(2.) ;
	static const Double_t rootPiBy2 = sqrt(atan2(0.0,-1.0)/2.0);
	Double_t xscale = root2/sigma;
	Double_t ret = 0;
	if(code==1){
		ret = rootPiBy2*height*(RooMath::erf((x.max(rangeName)-mean)/xscale) - RooMath::erf((x.min(rangeName)-mean)/xscale));
	} else if(code==2) {
		ret = rootPiBy2*height*(RooMath::erf((mean.max(rangeName)-x)/xscale) - RooMath::erf((mean.min(rangeName)-x)/xscale));
	} else{
		cout << "Error in RooGaussian::analyticalIntegral" << endl;
	}
	cout<<"roogaussvar integral \n";
	return ret;
}