Logo ROOT   6.08/07
Reference Guide
TFoamSampler.cxx
Go to the documentation of this file.
1 // @(#)root/unuran:$Id$
2 // Authors: L. Moneta, Dec 2010
3 
4 /**********************************************************************
5  * *
6  * Copyright (c) 2010 LCG ROOT Math Team, CERN/PH-SFT *
7  * *
8  * *
9  **********************************************************************/
10 
11 // Implementation file for class TFoamSampler
12 
13 #include "TFoamSampler.h"
15 
16 #include "TFoam.h"
17 #include "TFoamIntegrand.h"
19 #include "Math/IOptions.h"
20 #include "Fit/DataRange.h"
21 
22 #include "TRandom.h"
23 #include "TError.h"
24 
25 #include "TF1.h"
26 #include <cassert>
27 #include <cmath>
28 
29 class FoamDistribution : public TFoamIntegrand {
30 
31 public:
32 
33  FoamDistribution(const ROOT::Math::IMultiGenFunction & f, const ROOT::Fit::DataRange & range) :
34  fFunc(f),
35  fX(std::vector<double>(f.NDim() ) ),
36  fMinX(std::vector<double>(f.NDim() ) ),
37  fDeltaX(std::vector<double>(f.NDim() ) )
38  {
39  assert(f.NDim() == range.NDim() );
40  std::vector<double> xmax(f.NDim() );
41  for (unsigned int i = 0; i < range.NDim(); ++i) {
42  if (range.Size(i) == 0)
43  Error("FoamDistribution","Range is not set for coordinate dim %d",i);
44  else if (range.Size(i)>1)
45  Warning("FoamDistribution","Using only first range in coordinate dim %d",i);
46 
47  std::pair<double,double> r = range(i);
48  fMinX[i] = r.first;
49  fDeltaX[i] = r.second - r.first;
50  }
51  }
52  // in principle function does not need to be cloned
53 
54  virtual double Density(int ndim, double * x) {
55  assert(ndim == (int) fFunc.NDim() );
56  for (int i = 0; i < ndim; ++i)
57  fX[i] = fMinX[i] + x[i] * fDeltaX[i];
58 
59  return (fFunc)(&fX[0]);
60  }
61 
62  double MinX(unsigned int i) { return fMinX[i]; }
63  double DeltaX(unsigned int i) { return fDeltaX[i]; }
64 
65 private:
66 
67  const ROOT::Math::IMultiGenFunction & fFunc;
68  std::vector<double> fX;
69  std::vector<double> fMinX;
70  std::vector<double> fDeltaX;
71 
72 };
73 
74 
75 
77 
78 
79 ////////////////////////////////////////////////////////////////////////////////
80 
81 /**
82  TFoamSampler class
83  class implementing the ROOT::Math::DistSampler interface using FOAM
84  for sampling arbitrary distributions.
85 
86 
87 */
88 TFoamSampler::TFoamSampler() : ROOT::Math::DistSampler(),
89 // fOneDim(false)
90 // fDiscrete(false),
91 // fHasMode(false), fHasArea(false),
92 // fMode(0), fArea(0),
93  fFunc1D(0),
94  fFoam(new TFoam("FOAM") ),
95  fFoamDist(0)
96 {}
97 
99  assert(fFoam != 0);
100  delete fFoam;
101  if (fFoamDist) delete fFoamDist;
102 }
103 
104 bool TFoamSampler::Init(const char *) {
105 
106  // initialize using default options
109  if (foamOpt) opt.SetExtraOptions(*foamOpt);
110  return Init(opt);
111 }
112 
114  // initialize foam classes using the given algorithm
115  assert (fFoam != 0 );
116  if (NDim() == 0) {
117  Error("TFoamSampler::Init","Distribution function has not been set ! Need to call SetFunction first.");
118  return false;
119  }
120 
121  // initialize the foam
122  fFoam->SetkDim(NDim() );
123 
124  // initialize random number
125  if (!GetRandom()) SetRandom(gRandom);
126 
127  // create TFoamIntegrand class
128  if (fFoamDist) delete fFoamDist;
129  fFoamDist = new FoamDistribution(ParentPdf(),PdfRange());
130 
131  fFoam->SetRho(fFoamDist);
132  // set print level
133  fFoam->SetChat(opt.PrintLevel());
134 
135  // get extra options
136  ROOT::Math::IOptions * fopt = opt.ExtraOptions();
137  if (fopt) {
138  int nval = 0;
139  double fval = 0;
140  if (fopt->GetIntValue("nCells", nval) ) fFoam->SetnCells(nval);
141  if (fopt->GetIntValue("nCell1D", nval) && NDim() ==1) fFoam->SetnCells(nval);
142  if (fopt->GetIntValue("nCellND", nval) && NDim() >1) fFoam->SetnCells(nval);
143  if (fopt->GetIntValue("nCell2D", nval) && NDim() ==2) fFoam->SetnCells(nval);
144  if (fopt->GetIntValue("nCell3D", nval) && NDim() ==3) fFoam->SetnCells(nval);
145 
146  if (fopt->GetIntValue("nSample", nval) ) fFoam->SetnSampl(nval);
147  if (fopt->GetIntValue("nBin", nval) ) fFoam->SetnBin(nval);
148  if (fopt->GetIntValue("OptDrive",nval) ) fFoam->SetOptDrive(nval);
149  if (fopt->GetIntValue("OptRej",nval) ) fFoam->SetOptRej(nval);
150  if (fopt->GetRealValue("MaxWtRej",fval) ) fFoam->SetMaxWtRej(fval);
151 
152 
153  if (fopt->GetIntValue("chatLevel", nval) ) fFoam->SetChat(nval);
154  }
155  fFoam->Initialize();
156 
157  return true;
158 
159 }
160 
161 
163  // set function from a TF1 pointer
164  SetFunction<TF1>(*pdf, pdf->GetNdim());
165 }
166 
168  // set random generator (must be called before Init to have effect)
169  fFoam->SetPseRan(r);
170 }
171 
172 void TFoamSampler::SetSeed(unsigned int seed) {
173  // set random generator seed (must be called before Init to have effect)
174  TRandom * r = fFoam->GetPseRan();
175  if (r) r->SetSeed(seed);
176 }
177 
179  // get random generator used
180  return fFoam->GetPseRan();
181 }
182 
183 // double TFoamSampler::Sample1D() {
184 // // sample 1D distributions
185 // return (fDiscrete) ? (double) fFoam->SampleDiscr() : fFoam->Sample();
186 // }
187 
188 bool TFoamSampler::Sample(double * x) {
189  // sample multi-dim distributions
190 
191  fFoam->MakeEvent();
192  fFoam->GetMCvect(x);
193  // adjust for the range
194  for (unsigned int i = 0; i < NDim(); ++i)
195  x[i] = ( (FoamDistribution*)fFoamDist)->MinX(i) + ( ( (FoamDistribution*) fFoamDist)->DeltaX(i))*x[i];
196 
197  return true;
198 }
199 
200 
201 bool TFoamSampler::SampleBin(double prob, double & value, double *error) {
202  // sample a bin according to Poisson statistics
203 
204  TRandom * r = GetRandom();
205  if (!r) return false;
206  value = r->Poisson(prob);
207  if (error) *error = std::sqrt(value);
208  return true;
209 }
static ROOT::Math::IOptions * FindDefault(const char *name)
This namespace contains pre-defined functions to be used in conjuction with TExecutor::Map and TExecu...
Definition: StringConv.hxx:21
bool SampleBin(double prob, double &value, double *error=0)
sample one bin given an estimated of the pdf in the bin (this can be function value at the center or ...
const double * Sample()
sample one event and rerturning array x with coordinates
Definition: DistSampler.h:173
void SetFunction(const ROOT::Math::IGenFunction &func)
set the parent function distribution to use for random sampling (one dim case)
Definition: TFoamSampler.h:65
virtual ~TFoamSampler()
virtual destructor
STL namespace.
virtual bool GetRealValue(const char *, double &) const
Definition: IOptions.h:85
double sqrt(double)
Double_t x[n]
Definition: legend1.C:17
TRandom * GetRandom()
Get the random engine used by the sampler.
This is the base class for the ROOT Random number generators.
Definition: TRandom.h:31
virtual void SetSeed(ULong_t seed=0)
Set the random generator seed.
Definition: TRandom.cxx:568
void Init(TClassEdit::TInterpreterLookupHelper *helper)
Definition: TClassEdit.cxx:119
void SetExtraOptions(const IOptions &opt)
set extra options (in this case pointer is cloned)
virtual Int_t GetNdim() const
Definition: TF1.h:350
int PrintLevel() const
non-static methods for retrivieng options
TRandom2 r(17)
virtual Double_t Density(Int_t ndim, Double_t *)=0
DistSampler options class.
virtual unsigned int NDim() const =0
Retrieve the dimension of the function.
unsigned int Size(unsigned int icoord=0) const
return range size for coordinate icoord (starts from zero) Size == 0 indicates no range is present [-...
Definition: DataRange.h:70
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:925
bool Init(const char *="")
initialize the generators with the default options
virtual bool GetIntValue(const char *, int &) const
Definition: IOptions.h:87
float xmax
Definition: THbookFile.cxx:93
R__EXTERN TRandom * gRandom
Definition: TRandom.h:66
#define ClassImp(name)
Definition: Rtypes.h:279
class describing the range in the coordinates it supports multiple range in a coordinate.
Definition: DataRange.h:34
double f(double x)
Namespace for new Math classes and functions.
Generic interface for defining configuration options of a numerical algorithm.
Definition: IOptions.h:32
unsigned int NDim() const
get range dimension
Definition: DataRange.h:64
1-Dim function class
Definition: TF1.h:149
IOptions * ExtraOptions() const
return extra options (NULL pointer if they are not present)
virtual Int_t Poisson(Double_t mean)
Generates a random integer N according to a Poisson law.
Definition: TRandom.cxx:362
void SetRandom(TRandom *r)
Set the random engine to be used Needs to be called before Init to have effect.
Documentation for the abstract class IBaseFunctionMultiDim.
Definition: IFunction.h:63
void SetSeed(unsigned int seed)
Set the random seed for the TRandom instances used by the sampler classes Needs to be called before I...
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:911
Definition: TFoam.h:29
TFoamSampler class class implementing the ROOT::Math::DistSampler interface using FOAM for sampling a...
Definition: TFoamSampler.h:50