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