ROOT logo
/*****************************************************************************
 * Project: RooFit                                                           *
 * Package: RooFitCore                                                       *
 * @(#)root/roofitcore:$Id: RooGrid.cxx 24278 2008-06-15 15:21:16Z wouter $
 * Authors:                                                                  *
 *   WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu       *
 *   DK, David Kirkby,    UC Irvine,         dkirkby@uci.edu                 *
 *                                                                           *
 * Copyright (c) 2000-2005, Regents of the University of California          *
 *                          and Stanford University. All rights reserved.    *
 *                                                                           *
 * Redistribution and use in source and binary forms,                        *
 * with or without modification, are permitted according to the terms        *
 * listed in LICENSE (http://roofit.sourceforge.net/license.txt)             *
 *****************************************************************************/

//////////////////////////////////////////////////////////////////////////////
//
// BEGIN_HTML 
// RooGrid is a utility class for RooMCIntegrator which
// implements an adaptive multi-dimensional Monte Carlo numerical
// integration, following the VEGAS algorithm.  
// END_HTML
//

#include "RooFit.h"

#include "RooGrid.h"
#include "RooGrid.h"
#include "RooAbsFunc.h"
#include "RooNumber.h"
#include "RooRandom.h"
#include "TMath.h"
#include "RooMsgService.h"
#include "TClass.h"

#include <math.h>
#include "Riostream.h"
#include <iomanip>

ClassImp(RooGrid)
;


//_____________________________________________________________________________
RooGrid::RooGrid() :
  _xl(0),  _xu(0),  _delx(0),  _d(0),  _xi(0),  _xin(0),  _weight(0)
{
  // Default constructor
}


//_____________________________________________________________________________
RooGrid::RooGrid(const RooAbsFunc &function)
  : _valid(kTRUE), _xl(0),_xu(0),_delx(0),_xi(0)
{
  // Constructor with given function binding

  // check that the input function is valid
  if(!(_valid= function.isValid())) {
    oocoutE((TObject*)0,InputArguments) << ClassName() << ": cannot initialize using an invalid function" << endl;
    return;
  }

  // allocate workspace memory
  _dim= function.getDimension();
  _xl= new Double_t[_dim];
  _xu= new Double_t[_dim];
  _delx= new Double_t[_dim];
  _d= new Double_t[_dim*maxBins];
  _xi= new Double_t[_dim*(maxBins+1)];
  _xin= new Double_t[maxBins+1];
  _weight= new Double_t[maxBins];
  if(!_xl || !_xu || !_delx || !_d || !_xi || !_xin || !_weight) {
    oocoutE((TObject*)0,Integration) << ClassName() << ": memory allocation failed" << endl;
    _valid= kFALSE;
    return;
  }

  // initialize the grid
  _valid= initialize(function);
}


//_____________________________________________________________________________
RooGrid::~RooGrid() 
{
  // Destructor
  if(_xl)     delete[] _xl;
  if(_xu)     delete[] _xu;
  if(_delx)   delete[] _delx;
  if(_d)      delete[] _d;
  if(_xi)     delete[] _xi;
  if(_xin)    delete[] _xin;
  if(_weight) delete[] _weight;
}


//_____________________________________________________________________________
Bool_t RooGrid::initialize(const RooAbsFunc &function) 
{
  // Calculate and store the grid dimensions and volume using the
  // specified function, and initialize the grid using a single bin.
  // Return kTRUE, or else kFALSE if the range is not valid.

  _vol= 1;
  _bins= 1;
  for(UInt_t index= 0; index < _dim; index++) {
    _xl[index]= function.getMinLimit(index);
    if(RooNumber::isInfinite(_xl[index])) {
      oocoutE((TObject*)0,Integration) << ClassName() << ": lower limit of dimension " << index << " is infinite" << endl;
      return kFALSE;
    }
    _xu[index]= function.getMaxLimit(index);
    if(RooNumber::isInfinite(_xl[index])) {
      oocoutE((TObject*)0,Integration) << ClassName() << ": upper limit of dimension " << index << " is infinite" << endl;
      return kFALSE;
    }
    Double_t dx= _xu[index] - _xl[index];
    if(dx <= 0) {
      oocoutE((TObject*)0,Integration) << ClassName() << ": bad range for dimension " << index << ": [" << _xl[index]
				       << "," << _xu[index] << "]" << endl;
      return kFALSE;
    }
    _delx[index]= dx;
    _vol*= dx;
    coord(0,index) = 0;
    coord(1,index) = 1;
  }
  return kTRUE;
}


//_____________________________________________________________________________
void RooGrid::resize(UInt_t bins) 
{
  // Adjust the subdivision of each axis to give the specified
  // number of bins, using an algorithm that preserves relative
  // bin density. The new binning can be finer or coarser than
  // the original binning.

  // is there anything to do?
  if(bins == _bins) return;
  
  // weight is ratio of bin sizes
  Double_t pts_per_bin = (Double_t) _bins / (Double_t) bins;

  // loop over grid dimensions
  for (UInt_t j = 0; j < _dim; j++) {
    Double_t xold,xnew(0),dw(0);
    Int_t i = 1;
    // loop over bins in this dimension and load _xin[] with new bin edges

    UInt_t k;
    for(k = 1; k <= _bins; k++) {
      dw += 1.0;
      xold = xnew;
      xnew = coord(k,j);      
      while(dw > pts_per_bin) {
	dw -= pts_per_bin;
	newCoord(i++)= xnew - (xnew - xold) * dw;
      }
    }
    // copy the new edges into _xi[j]
    for(k = 1 ; k < bins; k++) {
      coord(k, j) = newCoord(k);
    }
    coord(bins, j) = 1;
  }  
  _bins = bins;
}


//_____________________________________________________________________________
void RooGrid::resetValues() 
{
  // Reset the values associated with each grid cell.
  
  for(UInt_t i = 0; i < _bins; i++) {
    for (UInt_t j = 0; j < _dim; j++) {
      value(i,j)= 0.0;
    }
  }
}


//_____________________________________________________________________________
void RooGrid::generatePoint(const UInt_t box[], Double_t x[], UInt_t bin[], Double_t &vol,
			    Bool_t useQuasiRandom) const 
{
  // Generate a random vector in the specified box and and store its
  // coordinates in the x[] array provided, the corresponding bin
  // indices in the bin[] array, and the volume of this bin in vol.
  // The box is specified by the array box[] of box integer indices
  // that each range from 0 to getNBoxes()-1.

  vol= 1;

  // generate a vector of quasi-random numbers to use
  if(useQuasiRandom) {
    RooRandom::quasi(_dim,x);
  }
  else {
    RooRandom::uniform(_dim,x);
  }

  // loop over coordinate axes
  for(UInt_t j= 0; j < _dim; ++j) {

    // generate a random point uniformly distributed (in box space)
    // within the box[j]-th box of coordinate axis j.
    Double_t z= ((box[j] + x[j])/_boxes)*_bins;

    // store the bin in which this point lies along the j-th
    // coordinate axis and calculate its width and position y
    // in normalized bin coordinates.
    Int_t k= (Int_t)z;
    bin[j] = k;
    Double_t y, bin_width;
    if(k == 0) {
      bin_width= coord(1,j);
      y= z * bin_width;
    }
    else {
      bin_width= coord(k+1,j) - coord(k,j);
      y= coord(k,j) + (z-k)*bin_width;
    }
    // transform from normalized bin coordinates to x space.
    x[j] = _xl[j] + y*_delx[j];

    // update this bin's calculated volume
    vol *= bin_width;
  }
}



//_____________________________________________________________________________
void RooGrid::firstBox(UInt_t box[]) const 
{
  // Reset the specified array of box indices to refer to the first box
  // in the standard traversal order.

  for(UInt_t i= 0; i < _dim; i++) box[i]= 0;
}



//_____________________________________________________________________________
Bool_t RooGrid::nextBox(UInt_t box[]) const 
{
  // Update the specified array of box indices to refer to the next box
  // in the standard traversal order and return kTRUE, or else return
  // kFALSE if we the indices already refer to the last box.

  // try incrementing each index until we find one that does not roll
  // over, starting from the last index.
  Int_t j(_dim-1);
  while (j >= 0) {
    box[j]= (box[j] + 1) % _boxes;
    if (0 != box[j]) return kTRUE;
    j--;
  }
  // if we get here, then there are no more boxes
  return kFALSE;
}



//_____________________________________________________________________________
void RooGrid::printMultiline(ostream& os, Int_t /*contents*/, Bool_t verbose, TString indent) const 
{
  // Print info about this object to the specified stream.

  os << ClassName() << ": volume = " << getVolume() << endl;
  os << indent << "  Has " << getDimension() << " dimension(s) each subdivided into "
     << getNBins() << " bin(s) and sampled with " << _boxes << " box(es)" << endl;
  for(UInt_t index= 0; index < getDimension(); index++) {
    os << indent << "  (" << index << ") ["
       << setw(10) << _xl[index] << "," << setw(10) << _xu[index] << "]" << endl;
    if(!verbose) continue;
    for(UInt_t bin= 0; bin < _bins; bin++) {
      os << indent << "    bin-" << bin << " : x = " << coord(bin,index) << " , y = "
	 << value(bin,index) << endl;
    }
  }
}


//_____________________________________________________________________________
void RooGrid::printName(ostream& os) const 
{
  // Print name of grid object
  os << GetName() ;
}


//_____________________________________________________________________________
void RooGrid::printTitle(ostream& os) const 
{
  // Print title of grid object
  os << GetTitle() ;
}


//_____________________________________________________________________________
void RooGrid::printClassName(ostream& os) const 
{
  // Print class name of grid object
  os << IsA()->GetName() ;
}



//_____________________________________________________________________________
void RooGrid::accumulate(const UInt_t bin[], Double_t amount) 
{
  // Add the specified amount to bin[j] of the 1D histograms associated
  // with each axis j.

  for(UInt_t j = 0; j < _dim; j++) value(bin[j],j) += amount;
}


//_____________________________________________________________________________
void RooGrid::refine(Double_t alpha) 
{
  // Refine the grid using the values that have been accumulated so far.
  // The parameter alpha controls the stiffness of the rebinning and should
  // usually be between 1 (stiffer) and 2 (more flexible). A value of zero
  // prevents any rebinning.
  
  for (UInt_t j = 0; j < _dim; j++) {

    // smooth this dimension's histogram of grid values and calculate the
    // new sum of the histogram contents as grid_tot_j
    Double_t oldg = value(0,j);
    Double_t newg = value(1,j);
    value(0,j)= (oldg + newg)/2;
    Double_t grid_tot_j = value(0,j);    
    // this loop implements value(i,j) = ( value(i-1,j)+value(i,j)+value(i+1,j) ) / 3

    UInt_t i;
    for (i = 1; i < _bins - 1; i++) {
      Double_t rc = oldg + newg;
      oldg = newg;
      newg = value(i+1,j);
      value(i,j)= (rc + newg)/3;
      grid_tot_j+= value(i,j);
    }
    value(_bins-1,j)= (newg + oldg)/2;
    grid_tot_j+= value(_bins-1,j);

    // calculate the weights for each bin of this dimension's histogram of values
    // and their sum
    Double_t tot_weight(0);
    for (i = 0; i < _bins; i++) {
      _weight[i] = 0;
      if (value(i,j) > 0) {
	oldg = grid_tot_j/value(i,j);
	/* damped change */
	_weight[i] = TMath::Power(((oldg-1.0)/oldg/log(oldg)), alpha);
      }
      tot_weight += _weight[i];
    }

    Double_t pts_per_bin = tot_weight / _bins;
    
    Double_t xold;
    Double_t xnew = 0;
    Double_t dw = 0;    

    UInt_t k;
    i = 1;
    for (k = 0; k < _bins; k++) {
      dw += _weight[k];
      xold = xnew;
      xnew = coord(k+1,j);
      
      while(dw > pts_per_bin) {
	dw -= pts_per_bin;
	newCoord(i++) = xnew - (xnew - xold) * dw / _weight[k];
      }
    }
    
    for (k = 1 ; k < _bins ; k++) {
      coord( k, j) = newCoord(k);
    }
    
    coord(_bins, j) = 1;
  }
}
 RooGrid.cxx:1
 RooGrid.cxx:2
 RooGrid.cxx:3
 RooGrid.cxx:4
 RooGrid.cxx:5
 RooGrid.cxx:6
 RooGrid.cxx:7
 RooGrid.cxx:8
 RooGrid.cxx:9
 RooGrid.cxx:10
 RooGrid.cxx:11
 RooGrid.cxx:12
 RooGrid.cxx:13
 RooGrid.cxx:14
 RooGrid.cxx:15
 RooGrid.cxx:16
 RooGrid.cxx:17
 RooGrid.cxx:18
 RooGrid.cxx:19
 RooGrid.cxx:20
 RooGrid.cxx:21
 RooGrid.cxx:22
 RooGrid.cxx:23
 RooGrid.cxx:24
 RooGrid.cxx:25
 RooGrid.cxx:26
 RooGrid.cxx:27
 RooGrid.cxx:28
 RooGrid.cxx:29
 RooGrid.cxx:30
 RooGrid.cxx:31
 RooGrid.cxx:32
 RooGrid.cxx:33
 RooGrid.cxx:34
 RooGrid.cxx:35
 RooGrid.cxx:36
 RooGrid.cxx:37
 RooGrid.cxx:38
 RooGrid.cxx:39
 RooGrid.cxx:40
 RooGrid.cxx:41
 RooGrid.cxx:42
 RooGrid.cxx:43
 RooGrid.cxx:44
 RooGrid.cxx:45
 RooGrid.cxx:46
 RooGrid.cxx:47
 RooGrid.cxx:48
 RooGrid.cxx:49
 RooGrid.cxx:50
 RooGrid.cxx:51
 RooGrid.cxx:52
 RooGrid.cxx:53
 RooGrid.cxx:54
 RooGrid.cxx:55
 RooGrid.cxx:56
 RooGrid.cxx:57
 RooGrid.cxx:58
 RooGrid.cxx:59
 RooGrid.cxx:60
 RooGrid.cxx:61
 RooGrid.cxx:62
 RooGrid.cxx:63
 RooGrid.cxx:64
 RooGrid.cxx:65
 RooGrid.cxx:66
 RooGrid.cxx:67
 RooGrid.cxx:68
 RooGrid.cxx:69
 RooGrid.cxx:70
 RooGrid.cxx:71
 RooGrid.cxx:72
 RooGrid.cxx:73
 RooGrid.cxx:74
 RooGrid.cxx:75
 RooGrid.cxx:76
 RooGrid.cxx:77
 RooGrid.cxx:78
 RooGrid.cxx:79
 RooGrid.cxx:80
 RooGrid.cxx:81
 RooGrid.cxx:82
 RooGrid.cxx:83
 RooGrid.cxx:84
 RooGrid.cxx:85
 RooGrid.cxx:86
 RooGrid.cxx:87
 RooGrid.cxx:88
 RooGrid.cxx:89
 RooGrid.cxx:90
 RooGrid.cxx:91
 RooGrid.cxx:92
 RooGrid.cxx:93
 RooGrid.cxx:94
 RooGrid.cxx:95
 RooGrid.cxx:96
 RooGrid.cxx:97
 RooGrid.cxx:98
 RooGrid.cxx:99
 RooGrid.cxx:100
 RooGrid.cxx:101
 RooGrid.cxx:102
 RooGrid.cxx:103
 RooGrid.cxx:104
 RooGrid.cxx:105
 RooGrid.cxx:106
 RooGrid.cxx:107
 RooGrid.cxx:108
 RooGrid.cxx:109
 RooGrid.cxx:110
 RooGrid.cxx:111
 RooGrid.cxx:112
 RooGrid.cxx:113
 RooGrid.cxx:114
 RooGrid.cxx:115
 RooGrid.cxx:116
 RooGrid.cxx:117
 RooGrid.cxx:118
 RooGrid.cxx:119
 RooGrid.cxx:120
 RooGrid.cxx:121
 RooGrid.cxx:122
 RooGrid.cxx:123
 RooGrid.cxx:124
 RooGrid.cxx:125
 RooGrid.cxx:126
 RooGrid.cxx:127
 RooGrid.cxx:128
 RooGrid.cxx:129
 RooGrid.cxx:130
 RooGrid.cxx:131
 RooGrid.cxx:132
 RooGrid.cxx:133
 RooGrid.cxx:134
 RooGrid.cxx:135
 RooGrid.cxx:136
 RooGrid.cxx:137
 RooGrid.cxx:138
 RooGrid.cxx:139
 RooGrid.cxx:140
 RooGrid.cxx:141
 RooGrid.cxx:142
 RooGrid.cxx:143
 RooGrid.cxx:144
 RooGrid.cxx:145
 RooGrid.cxx:146
 RooGrid.cxx:147
 RooGrid.cxx:148
 RooGrid.cxx:149
 RooGrid.cxx:150
 RooGrid.cxx:151
 RooGrid.cxx:152
 RooGrid.cxx:153
 RooGrid.cxx:154
 RooGrid.cxx:155
 RooGrid.cxx:156
 RooGrid.cxx:157
 RooGrid.cxx:158
 RooGrid.cxx:159
 RooGrid.cxx:160
 RooGrid.cxx:161
 RooGrid.cxx:162
 RooGrid.cxx:163
 RooGrid.cxx:164
 RooGrid.cxx:165
 RooGrid.cxx:166
 RooGrid.cxx:167
 RooGrid.cxx:168
 RooGrid.cxx:169
 RooGrid.cxx:170
 RooGrid.cxx:171
 RooGrid.cxx:172
 RooGrid.cxx:173
 RooGrid.cxx:174
 RooGrid.cxx:175
 RooGrid.cxx:176
 RooGrid.cxx:177
 RooGrid.cxx:178
 RooGrid.cxx:179
 RooGrid.cxx:180
 RooGrid.cxx:181
 RooGrid.cxx:182
 RooGrid.cxx:183
 RooGrid.cxx:184
 RooGrid.cxx:185
 RooGrid.cxx:186
 RooGrid.cxx:187
 RooGrid.cxx:188
 RooGrid.cxx:189
 RooGrid.cxx:190
 RooGrid.cxx:191
 RooGrid.cxx:192
 RooGrid.cxx:193
 RooGrid.cxx:194
 RooGrid.cxx:195
 RooGrid.cxx:196
 RooGrid.cxx:197
 RooGrid.cxx:198
 RooGrid.cxx:199
 RooGrid.cxx:200
 RooGrid.cxx:201
 RooGrid.cxx:202
 RooGrid.cxx:203
 RooGrid.cxx:204
 RooGrid.cxx:205
 RooGrid.cxx:206
 RooGrid.cxx:207
 RooGrid.cxx:208
 RooGrid.cxx:209
 RooGrid.cxx:210
 RooGrid.cxx:211
 RooGrid.cxx:212
 RooGrid.cxx:213
 RooGrid.cxx:214
 RooGrid.cxx:215
 RooGrid.cxx:216
 RooGrid.cxx:217
 RooGrid.cxx:218
 RooGrid.cxx:219
 RooGrid.cxx:220
 RooGrid.cxx:221
 RooGrid.cxx:222
 RooGrid.cxx:223
 RooGrid.cxx:224
 RooGrid.cxx:225
 RooGrid.cxx:226
 RooGrid.cxx:227
 RooGrid.cxx:228
 RooGrid.cxx:229
 RooGrid.cxx:230
 RooGrid.cxx:231
 RooGrid.cxx:232
 RooGrid.cxx:233
 RooGrid.cxx:234
 RooGrid.cxx:235
 RooGrid.cxx:236
 RooGrid.cxx:237
 RooGrid.cxx:238
 RooGrid.cxx:239
 RooGrid.cxx:240
 RooGrid.cxx:241
 RooGrid.cxx:242
 RooGrid.cxx:243
 RooGrid.cxx:244
 RooGrid.cxx:245
 RooGrid.cxx:246
 RooGrid.cxx:247
 RooGrid.cxx:248
 RooGrid.cxx:249
 RooGrid.cxx:250
 RooGrid.cxx:251
 RooGrid.cxx:252
 RooGrid.cxx:253
 RooGrid.cxx:254
 RooGrid.cxx:255
 RooGrid.cxx:256
 RooGrid.cxx:257
 RooGrid.cxx:258
 RooGrid.cxx:259
 RooGrid.cxx:260
 RooGrid.cxx:261
 RooGrid.cxx:262
 RooGrid.cxx:263
 RooGrid.cxx:264
 RooGrid.cxx:265
 RooGrid.cxx:266
 RooGrid.cxx:267
 RooGrid.cxx:268
 RooGrid.cxx:269
 RooGrid.cxx:270
 RooGrid.cxx:271
 RooGrid.cxx:272
 RooGrid.cxx:273
 RooGrid.cxx:274
 RooGrid.cxx:275
 RooGrid.cxx:276
 RooGrid.cxx:277
 RooGrid.cxx:278
 RooGrid.cxx:279
 RooGrid.cxx:280
 RooGrid.cxx:281
 RooGrid.cxx:282
 RooGrid.cxx:283
 RooGrid.cxx:284
 RooGrid.cxx:285
 RooGrid.cxx:286
 RooGrid.cxx:287
 RooGrid.cxx:288
 RooGrid.cxx:289
 RooGrid.cxx:290
 RooGrid.cxx:291
 RooGrid.cxx:292
 RooGrid.cxx:293
 RooGrid.cxx:294
 RooGrid.cxx:295
 RooGrid.cxx:296
 RooGrid.cxx:297
 RooGrid.cxx:298
 RooGrid.cxx:299
 RooGrid.cxx:300
 RooGrid.cxx:301
 RooGrid.cxx:302
 RooGrid.cxx:303
 RooGrid.cxx:304
 RooGrid.cxx:305
 RooGrid.cxx:306
 RooGrid.cxx:307
 RooGrid.cxx:308
 RooGrid.cxx:309
 RooGrid.cxx:310
 RooGrid.cxx:311
 RooGrid.cxx:312
 RooGrid.cxx:313
 RooGrid.cxx:314
 RooGrid.cxx:315
 RooGrid.cxx:316
 RooGrid.cxx:317
 RooGrid.cxx:318
 RooGrid.cxx:319
 RooGrid.cxx:320
 RooGrid.cxx:321
 RooGrid.cxx:322
 RooGrid.cxx:323
 RooGrid.cxx:324
 RooGrid.cxx:325
 RooGrid.cxx:326
 RooGrid.cxx:327
 RooGrid.cxx:328
 RooGrid.cxx:329
 RooGrid.cxx:330
 RooGrid.cxx:331
 RooGrid.cxx:332
 RooGrid.cxx:333
 RooGrid.cxx:334
 RooGrid.cxx:335
 RooGrid.cxx:336
 RooGrid.cxx:337
 RooGrid.cxx:338
 RooGrid.cxx:339
 RooGrid.cxx:340
 RooGrid.cxx:341
 RooGrid.cxx:342
 RooGrid.cxx:343
 RooGrid.cxx:344
 RooGrid.cxx:345
 RooGrid.cxx:346
 RooGrid.cxx:347
 RooGrid.cxx:348
 RooGrid.cxx:349
 RooGrid.cxx:350
 RooGrid.cxx:351
 RooGrid.cxx:352
 RooGrid.cxx:353
 RooGrid.cxx:354
 RooGrid.cxx:355
 RooGrid.cxx:356
 RooGrid.cxx:357
 RooGrid.cxx:358
 RooGrid.cxx:359
 RooGrid.cxx:360
 RooGrid.cxx:361
 RooGrid.cxx:362
 RooGrid.cxx:363
 RooGrid.cxx:364
 RooGrid.cxx:365
 RooGrid.cxx:366
 RooGrid.cxx:367
 RooGrid.cxx:368
 RooGrid.cxx:369
 RooGrid.cxx:370
 RooGrid.cxx:371
 RooGrid.cxx:372
 RooGrid.cxx:373
 RooGrid.cxx:374
 RooGrid.cxx:375
 RooGrid.cxx:376
 RooGrid.cxx:377
 RooGrid.cxx:378
 RooGrid.cxx:379
 RooGrid.cxx:380
 RooGrid.cxx:381
 RooGrid.cxx:382
 RooGrid.cxx:383
 RooGrid.cxx:384
 RooGrid.cxx:385
 RooGrid.cxx:386
 RooGrid.cxx:387
 RooGrid.cxx:388
 RooGrid.cxx:389
 RooGrid.cxx:390
 RooGrid.cxx:391
 RooGrid.cxx:392