// 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() :
_valid(kFALSE), _dim(0), _bins(0), _boxes(0), _vol(0), _xl(0), _xu(0), _delx(0), _d(0), _xi(0), _xin(0), _weight(0)
{
}
RooGrid::RooGrid(const RooAbsFunc &function)
: _valid(kTRUE), _xl(0),_xu(0),_delx(0),_xi(0)
{
if(!(_valid= function.isValid())) {
oocoutE((TObject*)0,InputArguments) << ClassName() << ": cannot initialize using an invalid function" << endl;
return;
}
_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;
}
_valid= initialize(function);
}
RooGrid::~RooGrid()
{
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)
{
_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)
{
if(bins == _bins) return;
Double_t pts_per_bin = (Double_t) _bins / (Double_t) bins;
for (UInt_t j = 0; j < _dim; j++) {
Double_t xold,xnew(0),dw(0);
Int_t i = 1;
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;
}
}
for(k = 1 ; k < bins; k++) {
coord(k, j) = newCoord(k);
}
coord(bins, j) = 1;
}
_bins = bins;
}
void RooGrid::resetValues()
{
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
{
vol= 1;
if(useQuasiRandom) {
RooRandom::quasi(_dim,x);
}
else {
RooRandom::uniform(_dim,x);
}
for(UInt_t j= 0; j < _dim; ++j) {
Double_t z= ((box[j] + x[j])/_boxes)*_bins;
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;
}
x[j] = _xl[j] + y*_delx[j];
vol *= bin_width;
}
}
void RooGrid::firstBox(UInt_t box[]) const
{
for(UInt_t i= 0; i < _dim; i++) box[i]= 0;
}
Bool_t RooGrid::nextBox(UInt_t box[]) const
{
Int_t j(_dim-1);
while (j >= 0) {
box[j]= (box[j] + 1) % _boxes;
if (0 != box[j]) return kTRUE;
j--;
}
return kFALSE;
}
void RooGrid::printMultiline(ostream& os, Int_t , Bool_t verbose, TString indent) const
{
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
{
os << GetName() ;
}
void RooGrid::printTitle(ostream& os) const
{
os << GetTitle() ;
}
void RooGrid::printClassName(ostream& os) const
{
os << IsA()->GetName() ;
}
void RooGrid::accumulate(const UInt_t bin[], Double_t amount)
{
for(UInt_t j = 0; j < _dim; j++) value(bin[j],j) += amount;
}
void RooGrid::refine(Double_t alpha)
{
for (UInt_t j = 0; j < _dim; 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);
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);
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);
_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;
}
}