// -- CLASS DESCRIPTION [PDF] --
// Johnson's S_U distribution
// 

#include "myRooJohnsonSU.h"
#include "TMath.h"

//ClassImp(myRooJohnsonSU)
  
//_____________________________________________________________________________
myRooJohnsonSU::myRooJohnsonSU(const char *name, const char *title,
			   RooAbsReal& _x,
			   RooAbsReal& _mu,
			   RooAbsReal& _sigma,
			   RooAbsReal& _delta,
			   RooAbsReal& _gamma) :
RooAbsPdf(name,title), 
  x("x","x",this,_x),
  mu("mu","mu",this,_mu),
  sigma("sigma","sigma",this,_sigma),
  delta("delta","delta",this,_delta),
  gamma("gamma","gamma",this,_gamma)
{;}

//_____________________________________________________________________________
myRooJohnsonSU::myRooJohnsonSU(const myRooJohnsonSU& other, const char* name) :
  RooAbsPdf(other,name), 
  x("x",this,other.x),
  mu("mu",this,other.mu),
  sigma("sigma",this,other.sigma),
  delta("delta",this,other.delta),
  gamma("gamma",this,other.gamma)
{;}

//_____________________________________________________________________________
Double_t myRooJohnsonSU::evaluate() const
{ 
	if(sigma<=0.) return 0.;
	double z     = (x-mu)/sigma;
	double arg   = gamma+delta*TMath::ASinH(z);
	double value = TMath::Exp(-arg*arg/2.)/TMath::Sqrt(1.+z*z);
	return value;
} 

//_____________________________________________________________________________
Double_t myRooJohnsonSU::analyticalIntegral(Int_t /*code*/, const char* rangeName) const
{	
	Double_t x_min = x.min(rangeName);
	Double_t x_max = x.max(rangeName);
	if (x_max<x_min) return 0.;
	
	if(delta==0.) return 0.;
	
	double arg_max = gamma+delta*TMath::ASinH((x_max-mu)/sigma);
	double arg_min = gamma+delta*TMath::ASinH((x_min-mu)/sigma);
	double norm = TMath::Sqrt(TMath::Pi()/2.)*sigma/delta*(TMath::Erf(arg_max/TMath::Sqrt2())-TMath::Erf(arg_min/TMath::Sqrt2()));
	return norm;
}