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 
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 "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 
41 using 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 
59 RooGrid::RooGrid(const RooAbsFunc &function)
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];
76  _weight= new Double_t[maxBins];
77  if(!_xl || !_xu || !_delx || !_d || !_xi || !_xin || !_weight) {
78  oocoutE((TObject*)0,Integration) << ClassName() << ": memory allocation failed" << endl;
79  _valid= kFALSE;
80  return;
81  }
82 
83  // initialize the grid
84  _valid= initialize(function);
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 
277 void 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 
297 void RooGrid::printName(ostream& os) const
298 {
299  os << GetName() ;
300 }
301 
302 
303 ////////////////////////////////////////////////////////////////////////////////
304 /// Print title of grid object
305 
306 void RooGrid::printTitle(ostream& os) const
307 {
308  os << GetTitle() ;
309 }
310 
311 
312 ////////////////////////////////////////////////////////////////////////////////
313 /// Print class name of grid object
314 
315 void 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 
326 void 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 }
RooGrid::_valid
Bool_t _valid
Definition: RooGrid.h:74
RooGrid::value
Double_t value(Int_t i, Int_t j) const
Definition: RooGrid.h:66
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:91
RooGrid::newCoord
Double_t & newCoord(Int_t i)
Definition: RooGrid.h:70
RooGrid.h
RooMsgService.h
RooGrid::accumulate
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
RooGrid::getDimension
UInt_t getDimension() const
Definition: RooGrid.h:41
RooFit.h
RooFit::InputArguments
@ InputArguments
Definition: RooGlobalFunc.h:68
RooGrid::_delx
Double_t * _delx
Internal workspace.
Definition: RooGrid.h:80
RooGrid::resetValues
void resetValues()
Reset the values associated with each grid cell.
Definition: RooGrid.cxx:181
ClassImp
#define ClassImp(name)
Definition: Rtypes.h:364
RooGrid::_d
Double_t * _d
Internal workspace.
Definition: RooGrid.h:81
RooGrid::_xu
Double_t * _xu
Internal workspace.
Definition: RooGrid.h:79
RooGrid::nextBox
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
log
double log(double)
Int_t
int Int_t
Definition: RtypesCore.h:45
box
void box(Int_t pat, Double_t x1, Double_t y1, Double_t x2, Double_t y2)
Definition: fillpatterns.C:1
x
Double_t x[n]
Definition: legend1.C:17
RooAbsFunc.h
TClass.h
indent
static void indent(ostringstream &buf, int indent_level)
Definition: TClingCallFunc.cxx:87
TObject::GetTitle
virtual const char * GetTitle() const
Returns title of object.
Definition: TObject.cxx:403
oocoutE
#define oocoutE(o, a)
Definition: RooMsgService.h:48
RooGrid::_bins
UInt_t _bins
Definition: RooGrid.h:75
TString
Definition: TString.h:136
RooGrid::isValid
Bool_t isValid() const
Definition: RooGrid.h:40
bool
RooAbsFunc
Definition: RooAbsFunc.h:23
TObject::GetName
virtual const char * GetName() const
Returns name of object.
Definition: TObject.cxx:359
RooGrid::initialize
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
RooGrid::_xi
Double_t * _xi
Internal workspace.
Definition: RooGrid.h:82
RooGrid::_xin
Double_t * _xin
Internal workspace.
Definition: RooGrid.h:83
RooGrid::_xl
Double_t * _xl
Definition: RooGrid.h:78
RooNumber::isInfinite
static Int_t isInfinite(Double_t x)
Return true if x is infinite by RooNumBer internal specification.
Definition: RooNumber.cxx:58
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:92
RooRandom.h
RooGrid::printMultiline
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
RooGrid::coord
Double_t coord(Int_t i, Int_t j) const
Definition: RooGrid.h:65
TMath::Power
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:724
RooGrid::refine
void refine(Double_t alpha=1.5)
Refine the grid using the values that have been accumulated so far.
Definition: RooGrid.cxx:338
RooGrid::firstBox
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
RooGrid::getNBins
UInt_t getNBins() const
Definition: RooGrid.h:43
RooGrid::printName
virtual void printName(std::ostream &os) const
Print name of grid object.
Definition: RooGrid.cxx:297
y
Double_t y[n]
Definition: legend1.C:17
RooGrid::generatePoint
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
unsigned int
RooGrid::resize
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
RooNumber.h
RooGrid::~RooGrid
virtual ~RooGrid()
Destructor.
Definition: RooGrid.cxx:91
RooGrid::_vol
Double_t _vol
Definition: RooGrid.h:76
Double_t
double Double_t
Definition: RtypesCore.h:59
RooRandom::quasi
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
RooGrid::printTitle
virtual void printTitle(std::ostream &os) const
Print title of grid object.
Definition: RooGrid.cxx:306
RooGrid::maxBins
@ maxBins
Definition: RooGrid.h:61
TObject
Definition: TObject.h:37
RooGrid::RooGrid
RooGrid()
Default constructor.
Definition: RooGrid.cxx:50
RooGrid::_boxes
UInt_t _boxes
Definition: RooGrid.h:75
genreflex::verbose
bool verbose
Definition: rootcling_impl.cxx:133
RooGrid
Definition: RooGrid.h:24
RooGrid::getVolume
Double_t getVolume() const
Definition: RooGrid.h:42
RooFit::Integration
@ Integration
Definition: RooGlobalFunc.h:67
Riostream.h
RooGrid::_weight
Double_t * _weight
Internal workspace.
Definition: RooGrid.h:84
RooRandom::uniform
static Double_t uniform(TRandom *generator=randomGenerator())
Return a number uniformly distributed from (0,1)
Definition: RooRandom.cxx:83
TObject::ClassName
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition: TObject.cxx:130
RooGrid::printClassName
virtual void printClassName(std::ostream &os) const
Print class name of grid object.
Definition: RooGrid.cxx:315
TMath.h
int
RooGrid::_dim
UInt_t _dim
Definition: RooGrid.h:75