Logo ROOT   6.18/05
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 "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
42using 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
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];
78 if(!_xl || !_xu || !_delx || !_d || !_xi || !_xin || !_weight) {
79 oocoutE((TObject*)0,Integration) << ClassName() << ": memory allocation failed" << endl;
81 return;
82 }
83
84 // initialize the grid
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
278void 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
298void RooGrid::printName(ostream& os) const
299{
300 os << GetName() ;
301}
302
303
304////////////////////////////////////////////////////////////////////////////////
305/// Print title of grid object
306
307void RooGrid::printTitle(ostream& os) const
308{
309 os << GetTitle() ;
310}
311
312
313////////////////////////////////////////////////////////////////////////////////
314/// Print class name of grid object
315
316void 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
327void 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}
#define oocoutE(o, a)
Definition: RooMsgService.h:47
int Int_t
Definition: RtypesCore.h:41
unsigned int UInt_t
Definition: RtypesCore.h:42
const Bool_t kFALSE
Definition: RtypesCore.h:88
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
const Bool_t kTRUE
Definition: RtypesCore.h:87
#define ClassImp(name)
Definition: Rtypes.h:365
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:109
@ maxBins
Definition: RooGrid.h:61
Double_t * _delx
Internal workspace.
Definition: RooGrid.h:80
virtual void printName(std::ostream &os) const
Print name of grid object.
Definition: RooGrid.cxx:298
virtual void printClassName(std::ostream &os) const
Print class name of grid object.
Definition: RooGrid.cxx:316
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:92
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:339
UInt_t _bins
Definition: RooGrid.h:75
virtual void printTitle(std::ostream &os) const
Print title of grid object.
Definition: RooGrid.cxx:307
Double_t * _xu
Internal workspace.
Definition: RooGrid.h:79
void resetValues()
Reset the values associated with each grid cell.
Definition: RooGrid.cxx:182
RooGrid()
Default constructor.
Definition: RooGrid.cxx:51
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:278
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:145
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
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:327
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
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
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:122
static Double_t uniform(TRandom *generator=randomGenerator())
Return a number uniformly distributed from (0,1)
Definition: RooRandom.cxx:84
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:58
@ Integration
Definition: RooGlobalFunc.h:57
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:723