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