// @(#)root/unuran:$Id: TUnuranContDist.cxx 28968 2009-06-12 16:41:43Z moneta $
// Authors: L. Moneta, J. Leydold Wed Feb 28 2007

/**********************************************************************
 *                                                                    *
 * Copyright (c) 2006  LCG ROOT Math Team, CERN/PH-SFT                *
 *                                                                    *
 *                                                                    *
 **********************************************************************/

// Implementation file for class TUnuranContDist

#include "TUnuranContDist.h"
#include "Math/RichardsonDerivator.h"
#include "Math/WrappedTF1.h"

#include "Math/Integrator.h"

#include "TF1.h"
#include <cassert>
#include <cmath>

ClassImp(TUnuranContDist)

TUnuranContDist::TUnuranContDist (const ROOT::Math::IGenFunction & pdf, const ROOT::Math::IGenFunction * deriv, bool isLogPdf  ) : 
   fPdf(&pdf),
   fDPdf(deriv),
   fCdf(0), 
   fXmin(1.), 
   fXmax(-1.), 
   fMode(0), 
   fArea(0),
   fIsLogPdf(isLogPdf),
   fHasDomain(0),
   fHasMode(0),
   fHasArea(0), 
   fOwnFunc(false)
{
   // Constructor from generic function interfaces
} 


TUnuranContDist::TUnuranContDist (TF1 * pdf, TF1 * deriv, bool isLogPdf  ) : 
   fPdf(  (pdf) ? new ROOT::Math::WrappedTF1 ( *pdf) : 0 ), 
   fDPdf( (deriv) ?  new ROOT::Math::WrappedTF1 ( *pdf) : 0 ),
   fCdf(0), 
   fXmin(1.), 
   fXmax(-1.), 
   fMode(0), 
   fArea(0),
   fIsLogPdf(isLogPdf),
   fHasDomain(0),
   fHasMode(0),
   fHasArea(0), 
   fOwnFunc(true)
{
   // Constructor from a TF1 objects
   // function pointers are managed by class 
} 


TUnuranContDist::TUnuranContDist(const TUnuranContDist & rhs) : 
   TUnuranBaseDist(), 
   fPdf(0), 
   fDPdf(0),
   fCdf(0)
{
   // Implementation of copy constructor 
   operator=(rhs);
}

TUnuranContDist & TUnuranContDist::operator = (const TUnuranContDist &rhs) 
{
   // Implementation of assignment operator.
   if (this == &rhs) return *this;  // time saving self-test
   fXmin  = rhs.fXmin;  
   fXmax  = rhs.fXmax;  
   fMode  = rhs.fMode; 
   fArea  = rhs.fArea;
   fIsLogPdf  = rhs.fIsLogPdf;
   fHasDomain = rhs.fHasDomain;
   fHasMode   = rhs.fHasMode;
   fHasArea   = rhs.fHasArea;
   fOwnFunc   = rhs.fOwnFunc;
   if (!fOwnFunc) { 
      fPdf   = rhs.fPdf;
      fDPdf  = rhs.fDPdf;
      fCdf   = rhs.fCdf; 
   }
   else {
      if (fPdf) delete fPdf;
      if (fDPdf) delete fDPdf; 
      if (fCdf) delete fCdf; 
      fPdf  = (rhs.fPdf)  ? rhs.fPdf->Clone()  : 0;  
      fDPdf = (rhs.fDPdf) ? rhs.fDPdf->Clone() : 0;  
      fCdf  = (rhs.fCdf)  ? rhs.fCdf->Clone()  : 0;  
   }

   return *this;
}

TUnuranContDist::~TUnuranContDist() { 
   // destructor implementation
   if (fOwnFunc) { 
      if (fPdf) delete fPdf; 
      if (fDPdf) delete fDPdf; 
      if (fCdf) delete fCdf; 
   }
}

void TUnuranContDist::SetCdf(const ROOT::Math::IGenFunction & cdf) {  
   //  set cdf distribution using a generic function interface
   fCdf = (fOwnFunc) ? cdf.Clone() : &cdf; 
}


void TUnuranContDist::SetCdf(TF1 *  cdf) { 
   // set cumulative distribution function from a TF1
   if (!fOwnFunc) { 
      // need to manage all functions now
      assert (fPdf != 0); 
      fPdf = fPdf->Clone(); 
      if (fDPdf) fDPdf->Clone(); 
   }
   else 
      if (fCdf) delete fCdf;

   fCdf = (cdf) ? new ROOT::Math::WrappedTF1 ( *cdf) : 0;    
   fOwnFunc = true; 
}

double TUnuranContDist::Pdf ( double x) const { 
   // evaluate the pdf of the distribution    
   assert(fPdf != 0);
   //fX[0] = x; 
   return (*fPdf)(x); 
}

double TUnuranContDist::DPdf( double x) const { 
   // evaluate the derivative of the pdf
   // if derivative function is not given is evaluated numerically
   if (fDPdf != 0) { 
      //fX[0] = x; 
      return (*fDPdf)(x); 
   }
   // do numerical derivation using numerical derivation
   ROOT::Math::RichardsonDerivator rd;
   static double gEps = 0.001;
   double h = ( std::abs(x) > 0 ) ?  gEps * std::abs(x) : gEps; 
   assert (fPdf != 0);
   return rd.Derivative1( *fPdf, x, h);
}

double TUnuranContDist::Cdf(double x) const {   
   // evaluate the integral (cdf)  on the domain
   if (fCdf != 0) {  
     // fX[0] = x; 
      return (*fCdf)(x);
   }
   // do numerical integration
   ROOT::Math::Integrator ig; 
   if (fXmin > fXmax) return ig.Integral( *fPdf ); 
   else 
      return ig.Integral( *fPdf, fXmin, fXmax );
   
}