testDistSampler.cxx

Go to the documentation of this file.

// program to test distribution sampling

#include "Math/Factory.h"
#include "Math/DistSampler.h"
#include "TF1.h"
#include "TH1.h"

#include "TStopwatch.h"
#include "TRandom.h"
#include "TError.h"
#include "TCanvas.h"
// double Pdf(double x) {
// }

using namespace ROOT::Math;


int testCont1D(int n)  {

   // test gaussian

   DistSampler * sampler = Factory::CreateDistSampler("Unuran");
   if (!sampler) return -1;
   TF1 * f = new TF1("pdf","gaus");
   f->SetParameters(1,0,1);

   TH1D * h1 = new TH1D("h1","h1",100,-3,3);
   TH1D * hr = new TH1D("hr","h2",100,-3,3);

   sampler->SetFunction(*f,1);
   bool ret = sampler->Init("AUTO");
   if (!ret)      return -1;


   TStopwatch w;
   w.Start();
   for (int i = 0; i < n; ++i) {
      h1->Fill(sampler->Sample1D() );
   }
   w.Stop();
   double c = 1.E9/double(n);
   std::cout << "Unuran sampling - (ns)/call = " << c*w.RealTime() << "   " << c*w.CpuTime() << std::endl;
   new TCanvas("Continous test");
   h1->SetLineColor(kBlue);
   h1->Draw();

   // generate ref histogram
   w.Start();
   for (int i = 0; i < n; ++i) {
      hr->Fill(gRandom->Gaus(0,1) );
   }
   w.Stop();
   std::cout << "TRandom::Gauss sampling - (ns)/call = " << c*w.RealTime() << "   " << c*w.CpuTime() << std::endl;

   hr->Draw("SAME");

   // do a Chi2 test
   // switch off printing of  info messages from chi2 test
   gErrorIgnoreLevel = 1001;
   double prob = h1->Chi2Test(hr,"UU");
   if (prob < 1.E-6) {
      std::cerr << "Chi2 test of generated histogram failed" << std::endl;
      return -2;
   }
   gErrorIgnoreLevel = 0;

   return 0;
}


int testDisc1D(int n)  {

   // test a discrete distribution


   DistSampler * sampler = Factory::CreateDistSampler("Unuran");
   if (!sampler) return -1;
   TF1 * f = new TF1("pdf","TMath::Poisson(x,[0])");
   double mu = 10;
   f->SetParameter(0,mu);

   TH1D * h1 = new TH1D("h2","h1",50,0,50);
   TH1D * hr = new TH1D("hd","h2",50,0,50);

   sampler->SetFunction(*f,1);
   sampler->SetMode(mu);
   sampler->SetArea(1);
   bool ret = sampler->Init("DARI");
   if (!ret)      return -1;


   TStopwatch w;
   w.Start();
   for (int i = 0; i < n; ++i) {
      h1->Fill(sampler->Sample1D() );
   }
   w.Stop();
   double c = 1.E9/double(n);
   std::cout << "Unuran sampling - (ns)/call = " << c*w.RealTime() << "   " << c*w.CpuTime() << std::endl;
   new TCanvas("Discrete test");
   h1->SetLineColor(kBlue);
   h1->Draw();

   // generate ref histogram
   w.Start();
   for (int i = 0; i < n; ++i) {
      hr->Fill(gRandom->Poisson(mu) );
   }
   w.Stop();
   std::cout << "TRandom::Poisson sampling - (ns)/call = " << c*w.RealTime() << "   " << c*w.CpuTime() << std::endl;

   hr->Draw("SAME");

   // do a Chi2 test
   // switch off printing of  info messages from chi2 test
   gErrorIgnoreLevel = 1001;
   double prob = h1->Chi2Test(hr,"UU");
   if (prob < 1.E-6) {
      std::cerr << "Chi2 test of generated histogram failed" << std::endl;
      return -2;
   }
   gErrorIgnoreLevel = 0;

   return 0;
}


int testDistSampler(int n = 10000) {
   int iret = 0;
   iret |= testCont1D(n);
   iret |= testDisc1D(n);
   return iret;
}
int main() {
   int iret = testDistSampler();
   if (iret)  std::cerr << "\ntestDistSampler: ....  FAILED!" << std::endl;
   else std::cerr << "\ntestDistSampler: ....  OK" << std::endl;
   return iret;
}