#include "TUnuranSampler.h"
#include "TUnuranContDist.h"
#include "TUnuranDiscrDist.h"
#include "TUnuranMultiContDist.h"
#include "TUnuran.h"
#include "Math/OneDimFunctionAdapter.h"
#include "Math/DistSamplerOptions.h"
#include "Fit/DataRange.h"
#include "TRandom.h"
#include "TError.h"
#include "TF1.h"
#include <cassert>
#include <cmath>
ClassImp(TUnuranSampler)
TUnuranSampler::TUnuranSampler() : ROOT::Math::DistSampler(),
fOneDim(false),
fDiscrete(false),
fHasMode(false), fHasArea(false),
fMode(0), fArea(0),
fFunc1D(0),
fUnuran(new TUnuran() )
{
fLevel = ROOT::Math::DistSamplerOptions::DefaultPrintLevel();
}
TUnuranSampler::~TUnuranSampler() {
assert(fUnuran != 0);
delete fUnuran;
}
bool TUnuranSampler::Init(const char * algo) {
assert (fUnuran != 0 );
if (NDim() == 0) {
Error("TUnuranSampler::Init","Distribution function has not been set ! Need to call SetFunction first.");
return false;
}
if (fLevel < 0) fLevel = ROOT::Math::DistSamplerOptions::DefaultPrintLevel();
TString method(algo);
if (method.IsNull() ) {
if (NDim() == 1) method = ROOT::Math::DistSamplerOptions::DefaultAlgorithm1D();
else method = ROOT::Math::DistSamplerOptions::DefaultAlgorithmND();
}
method.ToUpper();
bool ret = false;
if (NDim() == 1) {
if (method.First("D") == 0) {
if (fLevel>1) Info("TUnuranSampler::Init","Initialize one-dim discrete distribution with method %s",method.Data());
ret = DoInitDiscrete1D(method);
}
else {
if (fLevel>1) Info("TUnuranSampler::Init","Initialize one-dim continuous distribution with method %s",method.Data());
ret = DoInit1D(method);
}
}
else {
if (fLevel>1) Info("TUnuranSampler::Init","Initialize multi-dim continuous distribution with method %s",method.Data());
ret = DoInitND(method);
}
if (fLevel>0) {
if (ret) Info("TUnuranSampler::Init","Successfully initailized Unuran with method %s",method.Data() );
else Error("TUnuranSampler::Init","Failed to initailize Unuran with method %s",method.Data() );
}
return ret;
}
bool TUnuranSampler::Init(const ROOT::Math::DistSamplerOptions & opt ) {
SetPrintLevel(opt.PrintLevel() );
return Init(opt.Algorithm().c_str() );
}
bool TUnuranSampler::DoInit1D(const char * method) {
fOneDim = true;
TUnuranContDist * dist = 0;
if (fFunc1D == 0) {
ROOT::Math::OneDimMultiFunctionAdapter<> function(ParentPdf() );
dist = new TUnuranContDist(function,0,false,true);
}
else {
dist = new TUnuranContDist(*fFunc1D);
}
const ROOT::Fit::DataRange & range = PdfRange();
if (range.Size(0) > 0) {
double xmin, xmax;
range.GetRange(0,xmin,xmax);
dist->SetDomain(xmin,xmax);
}
if (fHasMode) dist->SetMode(fMode);
if (fHasArea) dist->SetPdfArea(fArea);
bool ret = false;
if (method) ret = fUnuran->Init(*dist, method);
else ret = fUnuran->Init(*dist);
delete dist;
return ret;
}
bool TUnuranSampler::DoInitDiscrete1D(const char * method) {
fOneDim = true;
fDiscrete = true;
TUnuranDiscrDist * dist = 0;
if (fFunc1D == 0) {
ROOT::Math::OneDimMultiFunctionAdapter<> function(ParentPdf() );
dist = new TUnuranDiscrDist(function,true);
}
else {
dist = new TUnuranDiscrDist(*fFunc1D, false);
}
const ROOT::Fit::DataRange & range = PdfRange();
if (range.Size(0) > 0) {
double xmin, xmax;
range.GetRange(0,xmin,xmax);
if (xmin < 0) {
Warning("DoInitDiscrete1D","range starts from negative values - set minimum to zero");
xmin = 0;
}
dist->SetDomain(int(xmin+0.1),int(xmax+0.1));
}
if (fHasMode) dist->SetMode(int(fMode+0.1));
if (fHasArea) dist->SetProbSum(fArea);
bool ret = fUnuran->Init(*dist, method);
delete dist;
return ret;
}
bool TUnuranSampler::DoInitND(const char * method) {
TUnuranMultiContDist dist(ParentPdf());
const ROOT::Fit::DataRange & range = PdfRange();
if (range.IsSet()) {
std::vector<double> xmin(range.NDim() );
std::vector<double> xmax(range.NDim() );
range.GetRange(&xmin[0],&xmax[0]);
dist.SetDomain(&xmin.front(),&xmax.front());
}
fOneDim = false;
if (method) return fUnuran->Init(dist, method);
return fUnuran->Init(dist);
}
void TUnuranSampler::SetFunction(TF1 * pdf) {
SetFunction<TF1>(*pdf, pdf->GetNdim());
}
void TUnuranSampler::SetRandom(TRandom * r) {
fUnuran->SetRandom(r);
}
void TUnuranSampler::SetSeed(unsigned int seed) {
fUnuran->SetSeed(seed);
}
TRandom * TUnuranSampler::GetRandom() {
return fUnuran->GetRandom();
}
double TUnuranSampler::Sample1D() {
return (fDiscrete) ? (double) fUnuran->SampleDiscr() : fUnuran->Sample();
}
bool TUnuranSampler::Sample(double * x) {
if (!fOneDim) return fUnuran->SampleMulti(x);
x[0] = Sample1D();
return true;
}
bool TUnuranSampler::SampleBin(double prob, double & value, double *error) {
TRandom * r = fUnuran->GetRandom();
if (!r) return false;
value = r->Poisson(prob);
if (error) *error = std::sqrt(value);
return true;
}