Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
RooGrid.cxx
Go to the documentation of this file.
1/// \cond ROOFIT_INTERNAL
2
3/*****************************************************************************
4 * Project: RooFit *
5 * Package: RooFitCore *
6 * @(#)root/roofitcore:$Id$
7 * Authors: *
8 * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
9 * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
10 * *
11 * Copyright (c) 2000-2005, Regents of the University of California *
12 * and Stanford University. All rights reserved. *
13 * *
14 * Redistribution and use in source and binary forms, *
15 * with or without modification, are permitted according to the terms *
16 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
17 *****************************************************************************/
18
19/**
20\file RooGrid.cxx
21\class RooGrid
22\ingroup Roofitcore
23
24Utility class for RooMCIntegrator which
25implements an adaptive multi-dimensional Monte Carlo numerical
26integration, following the VEGAS algorithm.
27**/
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
36#include <cmath>
37#include "Riostream.h"
38#include <iomanip>
39
40
41
42////////////////////////////////////////////////////////////////////////////////
43/// Constructor with given function binding
44
45RooGrid::RooGrid(const RooAbsFunc &function)
46 : _valid(true)
47{
48 // check that the input function is valid
49 if(!(_valid= function.isValid())) {
50 oocoutE(nullptr,InputArguments) << "RooGrid: cannot initialize using an invalid function" << std::endl;
51 return;
52 }
53
54 // allocate workspace memory
55 _dim= function.getDimension();
56 _xl.resize(_dim);
57 _xu.resize(_dim);
58 _delx.resize(_dim);
59 _d.resize(_dim*maxBins);
60 _xi.resize(_dim*(maxBins+1));
61 _xin.resize(maxBins+1);
62 _weight.resize(maxBins);
63
64 // initialize the grid
65 _valid= initialize(function);
66}
67
68
69////////////////////////////////////////////////////////////////////////////////
70/// Calculate and store the grid dimensions and volume using the
71/// specified function, and initialize the grid using a single bin.
72/// Return true, or else false if the range is not valid.
73
74bool RooGrid::initialize(const RooAbsFunc &function)
75{
76 _vol= 1;
77 _bins= 1;
78 for(UInt_t index= 0; index < _dim; index++) {
79 _xl[index]= function.getMinLimit(index);
81 oocoutE(nullptr,Integration) << "RooGrid: lower limit of dimension " << index << " is infinite" << std::endl;
82 return false;
83 }
84 _xu[index]= function.getMaxLimit(index);
86 oocoutE(nullptr,Integration) << "RooGrid: upper limit of dimension " << index << " is infinite" << std::endl;
87 return false;
88 }
89 double dx= _xu[index] - _xl[index];
90 if(dx <= 0) {
91 oocoutE(nullptr,Integration) << "RooGrid: bad range for dimension " << index << ": [" << _xl[index]
92 << "," << _xu[index] << "]" << std::endl;
93 return false;
94 }
95 _delx[index]= dx;
96 _vol*= dx;
97 coord(0,index) = 0;
98 coord(1,index) = 1;
99 }
100 return true;
101}
102
103
104////////////////////////////////////////////////////////////////////////////////
105/// Adjust the subdivision of each axis to give the specified
106/// number of bins, using an algorithm that preserves relative
107/// bin density. The new binning can be finer or coarser than
108/// the original binning.
109
110void RooGrid::resize(UInt_t bins)
111{
112 // is there anything to do?
113 if(bins == _bins) return;
114
115 // weight is ratio of bin sizes
116 double pts_per_bin = (double) _bins / (double) bins;
117
118 // loop over grid dimensions
119 for (UInt_t j = 0; j < _dim; j++) {
120 double xold;
121 double xnew(0);
122 double 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
149void RooGrid::resetValues()
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 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) {
173 RooRandom::quasi(_dim,x);
174 }
175 else {
176 RooRandom::uniform(_dim,x);
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= static_cast<Int_t>(z);
190 bin[j] = k;
191 double y;
192 double bin_width;
193 if(k == 0) {
194 bin_width= coord(1,j);
195 y= z * bin_width;
196 }
197 else {
198 bin_width= coord(k+1,j) - coord(k,j);
199 y= coord(k,j) + (z-k)*bin_width;
200 }
201 // transform from normalized bin coordinates to x space.
202 x[j] = _xl[j] + y*_delx[j];
203
204 // update this bin's calculated volume
205 vol *= bin_width;
206 }
207}
208
209
210
211////////////////////////////////////////////////////////////////////////////////
212/// Reset the specified array of box indices to refer to the first box
213/// in the standard traversal order.
214
215void RooGrid::firstBox(UInt_t box[]) const
216{
217 for(UInt_t i= 0; i < _dim; i++) box[i]= 0;
218}
219
220
221
222////////////////////////////////////////////////////////////////////////////////
223/// Update the specified array of box indices to refer to the next box
224/// in the standard traversal order and return true, or else return
225/// false if we the indices already refer to the last box.
226
227bool RooGrid::nextBox(UInt_t box[]) const
228{
229 // try incrementing each index until we find one that does not roll
230 // over, starting from the last index.
231 Int_t j(_dim-1);
232 while (j >= 0) {
233 box[j]= (box[j] + 1) % _boxes;
234 if (0 != box[j]) return true;
235 j--;
236 }
237 // if we get here, then there are no more boxes
238 return false;
239}
240
241
242
243////////////////////////////////////////////////////////////////////////////////
244/// Print info about this object to the specified stream.
245
246void RooGrid::print(std::ostream& os, bool verbose, std::string const& indent) const
247{
248 os << "RooGrid: volume = " << getVolume() << std::endl;
249 os << indent << " Has " << getDimension() << " dimension(s) each subdivided into "
250 << getNBins() << " bin(s) and sampled with " << _boxes << " box(es)" << std::endl;
251 for(std::size_t index= 0; index < getDimension(); index++) {
252 os << indent << " (" << index << ") ["
253 << std::setw(10) << _xl[index] << "," << std::setw(10) << _xu[index] << "]" << std::endl;
254 if(!verbose) continue;
255 for(std::size_t bin= 0; bin < _bins; bin++) {
256 os << indent << " bin-" << bin << " : x = " << coord(bin,index) << " , y = "
257 << value(bin,index) << std::endl;
258 }
259 }
260}
261
262
263////////////////////////////////////////////////////////////////////////////////
264/// Add the specified amount to bin[j] of the 1D histograms associated
265/// with each axis j.
266
267void RooGrid::accumulate(const UInt_t bin[], double amount)
268{
269 for(UInt_t j = 0; j < _dim; j++) value(bin[j],j) += amount;
270}
271
272
273////////////////////////////////////////////////////////////////////////////////
274/// Refine the grid using the values that have been accumulated so far.
275/// The parameter alpha controls the stiffness of the rebinning and should
276/// usually be between 1 (stiffer) and 2 (more flexible). A value of zero
277/// prevents any rebinning.
278
279void RooGrid::refine(double alpha)
280{
281 for (UInt_t j = 0; j < _dim; j++) {
282
283 // smooth this dimension's histogram of grid values and calculate the
284 // new sum of the histogram contents as grid_tot_j
285 double oldg = value(0,j);
286 double newg = value(1,j);
287 value(0,j)= (oldg + newg)/2;
288 double grid_tot_j = value(0,j);
289 // this loop implements value(i,j) = ( value(i-1,j)+value(i,j)+value(i+1,j) ) / 3
290
291 UInt_t i;
292 for (i = 1; i < _bins - 1; i++) {
293 double rc = oldg + newg;
294 oldg = newg;
295 newg = value(i+1,j);
296 value(i,j)= (rc + newg)/3;
297 grid_tot_j+= value(i,j);
298 }
299 value(_bins-1,j)= (newg + oldg)/2;
300 grid_tot_j+= value(_bins-1,j);
301
302 // calculate the weights for each bin of this dimension's histogram of values
303 // and their sum
304 double tot_weight(0);
305 for (i = 0; i < _bins; i++) {
306 _weight[i] = 0;
307 if (value(i,j) > 0) {
308 oldg = grid_tot_j/value(i,j);
309 /* damped change */
310 _weight[i] = TMath::Power(((oldg-1.0)/oldg/log(oldg)), alpha);
311 }
312 tot_weight += _weight[i];
313 }
314
315 double pts_per_bin = tot_weight / _bins;
316
317 double xold;
318 double xnew = 0;
319 double dw = 0;
320
321 i = 1;
322 for (UInt_t k = 0; k < _bins; k++) {
323 dw += _weight[k];
324 xold = xnew;
325 xnew = coord(k+1,j);
326
327 while(dw > pts_per_bin) {
328 dw -= pts_per_bin;
329 newCoord(i++) = xnew - (xnew - xold) * dw / _weight[k];
330 }
331 }
332
333 for (UInt_t k = 1 ; k < _bins ; k++) {
334 coord( k, j) = newCoord(k);
335 }
336
337 coord(_bins, j) = 1;
338 }
339}
340
341/// \endcond
RooAbsReal & function()
bool initialize()
#define oocoutE(o, a)
int Int_t
Definition RtypesCore.h:45
unsigned int UInt_t
Definition RtypesCore.h:46
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
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void value
Abstract interface for evaluating a real-valued function of one real variable and performing numerica...
Definition RooAbsFunc.h:27
bool isValid() const override
Check if current value is valid.
Definition RooAbsReal.h:452
static constexpr int isInfinite(double x)
Return true if x is infinite by RooNumber internal specification.
Definition RooNumber.h:27
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 ...
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:1804
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition TMath.h:721