TUnuranDiscrDist.cxx

Go to the documentation of this file.

// @(#)root/unuran:$Id: TUnuranDiscrDist.cxx 34077 2010-06-23 10:09:18Z moneta $
// Authors: L. Moneta, J. Leydold Wed Feb 28 2007 

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

// Implementation file for class TUnuranDiscrDist

#include "TUnuranDiscrDist.h"

#include "Math/IFunction.h"
#include "TF1.h"
#include "Math/WrappedTF1.h"


#include <cassert>


TUnuranDiscrDist::TUnuranDiscrDist (const ROOT::Math::IGenFunction & func, bool copyFunc) : 
   fPmf(&func), 
   fCdf(0), 
   fXmin(1), 
   fXmax(-1), 
   fMode(0), 
   fSum(0),
   fHasDomain(0),
   fHasMode(0),
   fHasSum(0), 
   fOwnFunc(copyFunc)
{
   //Constructor from a generic function object
   if (fOwnFunc) { 
      fPmf = fPmf->Clone(); 
      //if (fCdf) fCdf->Clone();
   }
} 


TUnuranDiscrDist::TUnuranDiscrDist (TF1 * func) : 
   fPmf( (func) ? new ROOT::Math::WrappedTF1 ( *func) : 0 ),
   fCdf(0), 
   fXmin(1), 
   fXmax(-1), 
   fMode(0), 
   fSum(0),
   fHasDomain(0),
   fHasMode(0),
   fHasSum(0), 
   fOwnFunc(true)
{
   //Constructor from a TF1 objects
} 


TUnuranDiscrDist::TUnuranDiscrDist(const TUnuranDiscrDist & rhs) :
   TUnuranBaseDist(), 
   fPmf(0), 
   fCdf(0)
{
   // Implementation of copy ctor using aassignment operator
   operator=(rhs);
}

TUnuranDiscrDist & TUnuranDiscrDist::operator = (const TUnuranDiscrDist &rhs) 
{
   // Implementation of assignment operator (copy only the function pointer not the function itself)
   if (this == &rhs) return *this;  // time saving self-test
   fPVec = rhs.fPVec;
   fPVecSum = rhs.fPVecSum;
   fXmin = rhs.fXmin;
   fXmax = rhs.fXmax;
   fMode = rhs.fMode;
   fSum  = rhs.fSum;
   fHasDomain = rhs.fHasDomain;
   fHasMode   = rhs.fHasMode;
   fHasSum    = rhs.fHasSum;
   fOwnFunc   = rhs.fOwnFunc;
   if (!fOwnFunc) { 
      fPmf   = rhs.fPmf;
      fCdf   = rhs.fCdf; 
   }
   else {
      if (fPmf) delete fPmf;
      if (fCdf) delete fCdf; 
      fPmf  = (rhs.fPmf)  ? rhs.fPmf->Clone()  : 0;  
      fCdf  = (rhs.fCdf)  ? rhs.fCdf->Clone()  : 0;  
   }

   return *this;
}

TUnuranDiscrDist::~TUnuranDiscrDist() { 
   // destructor implementation
   if (fOwnFunc) { 
      if (fPmf) delete fPmf; 
      if (fCdf) delete fCdf; 
   }
}

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

void TUnuranDiscrDist::SetCdf(TF1 *  cdf) { 
   // set cumulative distribution function from a TF1
   if (!fOwnFunc && fPmf) { 
      // need to manage also the pmf
      fPmf = fPmf->Clone(); 
   }
   else 
      if (fCdf) delete fCdf;

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

double TUnuranDiscrDist::Pmf ( int x) const {  
   // evaluate the distribution 
   if (!fPmf) { 
      if (x < static_cast<int>(fPVec.size()) || x >= static_cast<int>(fPVec.size()) ) return 0; 
      return fPVec[x]; 
   }
   return (*fPmf)(double(x)); 
}

double TUnuranDiscrDist::Cdf ( int x) const {  
   // evaluate the cumulative distribution 
   // otherwise evaluate from the sum of the probabilities 
   if (fHasDomain && x < fXmin) return 0; 

   if (fCdf) { 
      return (*fCdf)(double(x));
   }
 
   //estimation from sum of probability 
   int vsize = fPVecSum.size();
   if ( x < vsize ) 
      return fPVecSum[x]; 
   
   // calculate the sum 
   int x0 = ( fHasDomain) ? fXmin : 0; 
   int i0 = vsize;     // starting index
   int iN = x - x0 + 1;  // maximum index 
   fPVecSum.resize(iN); 
   double sum = ( i0 > 0 ) ? fPVecSum.back() : 0;
   for (int i = i0; i < iN; ++i) {
      sum += Pmf(i + x0); 
      fPVecSum[i] =  sum; 
   }

   return fPVecSum.back();
      
}

---