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