Logo ROOT   6.10/09
Reference Guide
RooGrid.cxx
Go to the documentation of this file.
1 /*****************************************************************************
2  * Project: RooFit *
3  * Package: RooFitCore *
4  * @(#)root/roofitcore:$Id$
5  * Authors: *
6  * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
7  * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
8  * *
9  * Copyright (c) 2000-2005, Regents of the University of California *
10  * and Stanford University. All rights reserved. *
11  * *
12  * Redistribution and use in source and binary forms, *
13  * with or without modification, are permitted according to the terms *
14  * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
15  *****************************************************************************/
16 
17 /**
18 \file RooGrid.cxx
19 \class RooGrid
20 \ingroup Roofitcore
21 
22 RooGrid is a utility class for RooMCIntegrator which
23 implements an adaptive multi-dimensional Monte Carlo numerical
24 integration, following the VEGAS algorithm.
25 **/
26 
27 #include "RooFit.h"
28 
29 #include "RooGrid.h"
30 #include "RooGrid.h"
31 #include "RooAbsFunc.h"
32 #include "RooNumber.h"
33 #include "RooRandom.h"
34 #include "TMath.h"
35 #include "RooMsgService.h"
36 #include "TClass.h"
37 
38 #include <math.h>
39 #include "Riostream.h"
40 #include <iomanip>
41 
42 using namespace std;
43 
45 ;
46 
47 
48 ////////////////////////////////////////////////////////////////////////////////
49 /// Default constructor
50 
52  _valid(kFALSE), _dim(0), _bins(0), _boxes(0), _vol(0), _xl(0), _xu(0), _delx(0), _d(0), _xi(0), _xin(0), _weight(0)
53 {
54 }
55 
56 
57 ////////////////////////////////////////////////////////////////////////////////
58 /// Constructor with given function binding
59 
60 RooGrid::RooGrid(const RooAbsFunc &function)
61  : _valid(kTRUE), _xl(0),_xu(0),_delx(0),_xi(0)
62 {
63  // check that the input function is valid
64  if(!(_valid= function.isValid())) {
65  oocoutE((TObject*)0,InputArguments) << ClassName() << ": cannot initialize using an invalid function" << endl;
66  return;
67  }
68 
69  // allocate workspace memory
70  _dim= function.getDimension();
71  _xl= new Double_t[_dim];
72  _xu= new Double_t[_dim];
73  _delx= new Double_t[_dim];
74  _d= new Double_t[_dim*maxBins];
75  _xi= new Double_t[_dim*(maxBins+1)];
76  _xin= new Double_t[maxBins+1];
77  _weight= new Double_t[maxBins];
78  if(!_xl || !_xu || !_delx || !_d || !_xi || !_xin || !_weight) {
79  oocoutE((TObject*)0,Integration) << ClassName() << ": memory allocation failed" << endl;
80  _valid= kFALSE;
81  return;
82  }
83 
84  // initialize the grid
85  _valid= initialize(function);
86 }
87 
88 
89 ////////////////////////////////////////////////////////////////////////////////
90 /// Destructor
91 
93 {
94  if(_xl) delete[] _xl;
95  if(_xu) delete[] _xu;
96  if(_delx) delete[] _delx;
97  if(_d) delete[] _d;
98  if(_xi) delete[] _xi;
99  if(_xin) delete[] _xin;
100  if(_weight) delete[] _weight;
101 }
102 
103 
104 ////////////////////////////////////////////////////////////////////////////////
105 /// Calculate and store the grid dimensions and volume using the
106 /// specified function, and initialize the grid using a single bin.
107 /// Return kTRUE, or else kFALSE if the range is not valid.
108 
110 {
111  _vol= 1;
112  _bins= 1;
113  for(UInt_t index= 0; index < _dim; index++) {
114  _xl[index]= function.getMinLimit(index);
115  if(RooNumber::isInfinite(_xl[index])) {
116  oocoutE((TObject*)0,Integration) << ClassName() << ": lower limit of dimension " << index << " is infinite" << endl;
117  return kFALSE;
118  }
119  _xu[index]= function.getMaxLimit(index);
120  if(RooNumber::isInfinite(_xl[index])) {
121  oocoutE((TObject*)0,Integration) << ClassName() << ": upper limit of dimension " << index << " is infinite" << endl;
122  return kFALSE;
123  }
124  Double_t dx= _xu[index] - _xl[index];
125  if(dx <= 0) {
126  oocoutE((TObject*)0,Integration) << ClassName() << ": bad range for dimension " << index << ": [" << _xl[index]
127  << "," << _xu[index] << "]" << endl;
128  return kFALSE;
129  }
130  _delx[index]= dx;
131  _vol*= dx;
132  coord(0,index) = 0;
133  coord(1,index) = 1;
134  }
135  return kTRUE;
136 }
137 
138 
139 ////////////////////////////////////////////////////////////////////////////////
140 /// Adjust the subdivision of each axis to give the specified
141 /// number of bins, using an algorithm that preserves relative
142 /// bin density. The new binning can be finer or coarser than
143 /// the original binning.
144 
146 {
147  // is there anything to do?
148  if(bins == _bins) return;
149 
150  // weight is ratio of bin sizes
151  Double_t pts_per_bin = (Double_t) _bins / (Double_t) bins;
152 
153  // loop over grid dimensions
154  for (UInt_t j = 0; j < _dim; j++) {
155  Double_t xold,xnew(0),dw(0);
156  Int_t i = 1;
157  // loop over bins in this dimension and load _xin[] with new bin edges
158 
159  UInt_t k;
160  for(k = 1; k <= _bins; k++) {
161  dw += 1.0;
162  xold = xnew;
163  xnew = coord(k,j);
164  while(dw > pts_per_bin) {
165  dw -= pts_per_bin;
166  newCoord(i++)= xnew - (xnew - xold) * dw;
167  }
168  }
169  // copy the new edges into _xi[j]
170  for(k = 1 ; k < bins; k++) {
171  coord(k, j) = newCoord(k);
172  }
173  coord(bins, j) = 1;
174  }
175  _bins = bins;
176 }
177 
178 
179 ////////////////////////////////////////////////////////////////////////////////
180 /// Reset the values associated with each grid cell.
181 
183 {
184  for(UInt_t i = 0; i < _bins; i++) {
185  for (UInt_t j = 0; j < _dim; j++) {
186  value(i,j)= 0.0;
187  }
188  }
189 }
190 
191 
192 ////////////////////////////////////////////////////////////////////////////////
193 /// Generate a random vector in the specified box and and store its
194 /// coordinates in the x[] array provided, the corresponding bin
195 /// indices in the bin[] array, and the volume of this bin in vol.
196 /// The box is specified by the array box[] of box integer indices
197 /// that each range from 0 to getNBoxes()-1.
198 
200  Bool_t useQuasiRandom) const
201 {
202  vol= 1;
203 
204  // generate a vector of quasi-random numbers to use
205  if(useQuasiRandom) {
207  }
208  else {
210  }
211 
212  // loop over coordinate axes
213  for(UInt_t j= 0; j < _dim; ++j) {
214 
215  // generate a random point uniformly distributed (in box space)
216  // within the box[j]-th box of coordinate axis j.
217  Double_t z= ((box[j] + x[j])/_boxes)*_bins;
218 
219  // store the bin in which this point lies along the j-th
220  // coordinate axis and calculate its width and position y
221  // in normalized bin coordinates.
222  Int_t k= (Int_t)z;
223  bin[j] = k;
224  Double_t y, bin_width;
225  if(k == 0) {
226  bin_width= coord(1,j);
227  y= z * bin_width;
228  }
229  else {
230  bin_width= coord(k+1,j) - coord(k,j);
231  y= coord(k,j) + (z-k)*bin_width;
232  }
233  // transform from normalized bin coordinates to x space.
234  x[j] = _xl[j] + y*_delx[j];
235 
236  // update this bin's calculated volume
237  vol *= bin_width;
238  }
239 }
240 
241 
242 
243 ////////////////////////////////////////////////////////////////////////////////
244 /// Reset the specified array of box indices to refer to the first box
245 /// in the standard traversal order.
246 
248 {
249  for(UInt_t i= 0; i < _dim; i++) box[i]= 0;
250 }
251 
252 
253 
254 ////////////////////////////////////////////////////////////////////////////////
255 /// Update the specified array of box indices to refer to the next box
256 /// in the standard traversal order and return kTRUE, or else return
257 /// kFALSE if we the indices already refer to the last box.
258 
260 {
261  // try incrementing each index until we find one that does not roll
262  // over, starting from the last index.
263  Int_t j(_dim-1);
264  while (j >= 0) {
265  box[j]= (box[j] + 1) % _boxes;
266  if (0 != box[j]) return kTRUE;
267  j--;
268  }
269  // if we get here, then there are no more boxes
270  return kFALSE;
271 }
272 
273 
274 
275 ////////////////////////////////////////////////////////////////////////////////
276 /// Print info about this object to the specified stream.
277 
278 void RooGrid::printMultiline(ostream& os, Int_t /*contents*/, Bool_t verbose, TString indent) const
279 {
280  os << ClassName() << ": volume = " << getVolume() << endl;
281  os << indent << " Has " << getDimension() << " dimension(s) each subdivided into "
282  << getNBins() << " bin(s) and sampled with " << _boxes << " box(es)" << endl;
283  for(UInt_t index= 0; index < getDimension(); index++) {
284  os << indent << " (" << index << ") ["
285  << setw(10) << _xl[index] << "," << setw(10) << _xu[index] << "]" << endl;
286  if(!verbose) continue;
287  for(UInt_t bin= 0; bin < _bins; bin++) {
288  os << indent << " bin-" << bin << " : x = " << coord(bin,index) << " , y = "
289  << value(bin,index) << endl;
290  }
291  }
292 }
293 
294 
295 ////////////////////////////////////////////////////////////////////////////////
296 /// Print name of grid object
297 
298 void RooGrid::printName(ostream& os) const
299 {
300  os << GetName() ;
301 }
302 
303 
304 ////////////////////////////////////////////////////////////////////////////////
305 /// Print title of grid object
306 
307 void RooGrid::printTitle(ostream& os) const
308 {
309  os << GetTitle() ;
310 }
311 
312 
313 ////////////////////////////////////////////////////////////////////////////////
314 /// Print class name of grid object
315 
316 void RooGrid::printClassName(ostream& os) const
317 {
318  os << IsA()->GetName() ;
319 }
320 
321 
322 
323 ////////////////////////////////////////////////////////////////////////////////
324 /// Add the specified amount to bin[j] of the 1D histograms associated
325 /// with each axis j.
326 
327 void RooGrid::accumulate(const UInt_t bin[], Double_t amount)
328 {
329  for(UInt_t j = 0; j < _dim; j++) value(bin[j],j) += amount;
330 }
331 
332 
333 ////////////////////////////////////////////////////////////////////////////////
334 /// Refine the grid using the values that have been accumulated so far.
335 /// The parameter alpha controls the stiffness of the rebinning and should
336 /// usually be between 1 (stiffer) and 2 (more flexible). A value of zero
337 /// prevents any rebinning.
338 
340 {
341  for (UInt_t j = 0; j < _dim; j++) {
342 
343  // smooth this dimension's histogram of grid values and calculate the
344  // new sum of the histogram contents as grid_tot_j
345  Double_t oldg = value(0,j);
346  Double_t newg = value(1,j);
347  value(0,j)= (oldg + newg)/2;
348  Double_t grid_tot_j = value(0,j);
349  // this loop implements value(i,j) = ( value(i-1,j)+value(i,j)+value(i+1,j) ) / 3
350 
351  UInt_t i;
352  for (i = 1; i < _bins - 1; i++) {
353  Double_t rc = oldg + newg;
354  oldg = newg;
355  newg = value(i+1,j);
356  value(i,j)= (rc + newg)/3;
357  grid_tot_j+= value(i,j);
358  }
359  value(_bins-1,j)= (newg + oldg)/2;
360  grid_tot_j+= value(_bins-1,j);
361 
362  // calculate the weights for each bin of this dimension's histogram of values
363  // and their sum
364  Double_t tot_weight(0);
365  for (i = 0; i < _bins; i++) {
366  _weight[i] = 0;
367  if (value(i,j) > 0) {
368  oldg = grid_tot_j/value(i,j);
369  /* damped change */
370  _weight[i] = TMath::Power(((oldg-1.0)/oldg/log(oldg)), alpha);
371  }
372  tot_weight += _weight[i];
373  }
374 
375  Double_t pts_per_bin = tot_weight / _bins;
376 
377  Double_t xold;
378  Double_t xnew = 0;
379  Double_t dw = 0;
380 
381  UInt_t k;
382  i = 1;
383  for (k = 0; k < _bins; k++) {
384  dw += _weight[k];
385  xold = xnew;
386  xnew = coord(k+1,j);
387 
388  while(dw > pts_per_bin) {
389  dw -= pts_per_bin;
390  newCoord(i++) = xnew - (xnew - xold) * dw / _weight[k];
391  }
392  }
393 
394  for (k = 1 ; k < _bins ; k++) {
395  coord( k, j) = newCoord(k);
396  }
397 
398  coord(_bins, j) = 1;
399  }
400 }
virtual void printClassName(std::ostream &os) const
Print class name of grid object.
Definition: RooGrid.cxx:316
UInt_t getDimension() const
Definition: RooGrid.h:41
virtual void printTitle(std::ostream &os) const
Print title of grid object.
Definition: RooGrid.cxx:307
UInt_t _bins
Definition: RooGrid.h:75
RooGrid()
Default constructor.
Definition: RooGrid.cxx:51
Bool_t nextBox(UInt_t box[]) const
Update the specified array of box indices to refer to the next box in the standard traversal order an...
Definition: RooGrid.cxx:259
Double_t & newCoord(Int_t i)
Definition: RooGrid.h:70
Double_t * _weight
Internal workspace.
Definition: RooGrid.h:84
Double_t * _delx
Internal workspace.
Definition: RooGrid.h:80
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
static Int_t isInfinite(Double_t x)
Return true if x is infinite by RooNumBer internal specification.
Definition: RooNumber.cxx:58
static Bool_t quasi(UInt_t dimension, Double_t vector[], RooQuasiRandomGenerator *generator=quasiGenerator())
Return a quasi-random number in the range (0,1) using the Niederreiter base 2 generator described in ...
Definition: RooRandom.cxx:122
STL namespace.
Double_t * _xu
Internal workspace.
Definition: RooGrid.h:79
void box(Int_t pat, Double_t x1, Double_t y1, Double_t x2, Double_t y2)
Definition: fillpatterns.C:1
Bool_t _valid
Definition: RooGrid.h:74
Double_t * _d
Internal workspace.
Definition: RooGrid.h:81
Double_t coord(Int_t i, Int_t j) const
Definition: RooGrid.h:65
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:628
virtual void printMultiline(std::ostream &os, Int_t contents, Bool_t verbose=kFALSE, TString indent="") const
Print info about this object to the specified stream.
Definition: RooGrid.cxx:278
UInt_t getNBins() const
Definition: RooGrid.h:43
void generatePoint(const UInt_t box[], Double_t x[], UInt_t bin[], Double_t &vol, Bool_t useQuasiRandom=kTRUE) const
Generate a random vector in the specified box and and store its coordinates in the x[] array provided...
Definition: RooGrid.cxx:199
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition: TObject.cxx:135
Double_t x[n]
Definition: legend1.C:17
#define oocoutE(o, a)
Definition: RooMsgService.h:47
Bool_t isValid() const
Definition: RooGrid.h:40
Double_t _vol
Definition: RooGrid.h:76
virtual ~RooGrid()
Destructor.
Definition: RooGrid.cxx:92
UInt_t _dim
Definition: RooGrid.h:75
void firstBox(UInt_t box[]) const
Reset the specified array of box indices to refer to the first box in the standard traversal order...
Definition: RooGrid.cxx:247
Double_t * _xin
Internal workspace.
Definition: RooGrid.h:83
unsigned int UInt_t
Definition: RtypesCore.h:42
bool verbose
const Bool_t kFALSE
Definition: RtypesCore.h:92
void refine(Double_t alpha=1.5)
Refine the grid using the values that have been accumulated so far.
Definition: RooGrid.cxx:339
#define ClassImp(name)
Definition: Rtypes.h:336
static Double_t uniform(TRandom *generator=randomGenerator())
Return a number uniformly distributed from (0,1)
Definition: RooRandom.cxx:84
Double_t value(Int_t i, Int_t j) const
Definition: RooGrid.h:66
double Double_t
Definition: RtypesCore.h:55
Bool_t initialize(const RooAbsFunc &function)
Calculate and store the grid dimensions and volume using the specified function, and initialize the g...
Definition: RooGrid.cxx:109
UInt_t _boxes
Definition: RooGrid.h:75
Double_t * _xi
Internal workspace.
Definition: RooGrid.h:82
Double_t y[n]
Definition: legend1.C:17
Double_t * _xl
Definition: RooGrid.h:78
void accumulate(const UInt_t bin[], Double_t amount)
Add the specified amount to bin[j] of the 1D histograms associated with each axis j...
Definition: RooGrid.cxx:327
Mother of all ROOT objects.
Definition: TObject.h:37
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
virtual const char * GetTitle() const
Returns title of object.
Definition: TObject.cxx:408
void resetValues()
Reset the values associated with each grid cell.
Definition: RooGrid.cxx:182
Double_t getVolume() const
Definition: RooGrid.h:42
virtual void printName(std::ostream &os) const
Print name of grid object.
Definition: RooGrid.cxx:298
virtual const char * GetName() const
Returns name of object.
Definition: TObject.cxx:364
void resize(UInt_t bins)
Adjust the subdivision of each axis to give the specified number of bins, using an algorithm that pre...
Definition: RooGrid.cxx:145
const Bool_t kTRUE
Definition: RtypesCore.h:91
RooGrid is a utility class for RooMCIntegrator which implements an adaptive multi-dimensional Monte C...
Definition: RooGrid.h:24
Abstract interface for evaluating a real-valued function of one real variable and performing numerica...
Definition: RooAbsFunc.h:23
double log(double)