/*****************************************************************************
 * Project: RooFit                                                           *
 * Package: RooFitCore                                                       *
 * @(#)root/roofitcore:$Id$
 * 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
// RooBinIntegrator implements an adaptive one-dimensional 
// numerical integration algorithm. 
// END_HTML
//


#include "RooFit.h"
#include "Riostream.h"

#include "TClass.h"
#include "RooBinIntegrator.h"
#include "RooArgSet.h"
#include "RooRealVar.h"
#include "RooNumber.h"
#include "RooIntegratorBinding.h"
#include "RooNumIntConfig.h"
#include "RooNumIntFactory.h"
#include "RooMsgService.h"

#include <assert.h>



using namespace std;

ClassImp(RooBinIntegrator)
;

// Register this class with RooNumIntConfig

//_____________________________________________________________________________
void RooBinIntegrator::registerIntegrator(RooNumIntFactory& fact)
{
  // Register RooBinIntegrator, is parameters and capabilities with RooNumIntFactory

  RooRealVar numBins("numBins","Number of bins in range",100) ;
  RooBinIntegrator* proto = new RooBinIntegrator() ;
  fact.storeProtoIntegrator(proto,RooArgSet(numBins)) ;
  RooNumIntConfig::defaultConfig().method1D().setLabel(proto->IsA()->GetName()) ;
}



//_____________________________________________________________________________
RooBinIntegrator::RooBinIntegrator() : _numBins(0), _useIntegrandLimits(kFALSE), _x(0)
{
  // Default constructor
}


//_____________________________________________________________________________
RooBinIntegrator::RooBinIntegrator(const RooAbsFunc& function) : 
  RooAbsIntegrator(function)
{
  // Construct integrator on given function binding binding

  _useIntegrandLimits= kTRUE;
  assert(0 != integrand() && integrand()->isValid());

  // Allocate coordinate buffer size after number of function dimensions
  _x = new Double_t[_function->getDimension()] ;
  _numBins = 100 ;

  _xmin.resize(_function->getDimension()) ;
  _xmax.resize(_function->getDimension()) ;

  for (UInt_t i=0 ; i<_function->getDimension() ; i++) {
    _xmin[i]= integrand()->getMinLimit(i);
    _xmax[i]= integrand()->getMaxLimit(i);

    // Retrieve bin configuration from integrand
    list<Double_t>* tmp = integrand()->binBoundaries(i) ;
    if (!tmp) {
      oocoutW((TObject*)0,Integration) << "RooBinIntegrator::RooBinIntegrator WARNING: integrand provide no binning definition observable #" 
				     << i << " substituting default binning of " << _numBins << " bins" << endl ;
      tmp = new list<Double_t> ;
      for (Int_t j=0 ; j<=_numBins ; j++) {
	tmp->push_back(_xmin[i]+j*(_xmax[i]-_xmin[i])/_numBins) ;
      }
    }
    _binb.push_back(tmp) ;
  }
  checkLimits();

} 


//_____________________________________________________________________________
RooBinIntegrator::RooBinIntegrator(const RooAbsFunc& function, const RooNumIntConfig& config) : 
  RooAbsIntegrator(function), _binb(0)
{
  // Construct integrator on given function binding binding
  
  const RooArgSet& configSet = config.getConfigSection(IsA()->GetName()) ;  
  _useIntegrandLimits= kTRUE;
  _numBins = (Int_t) configSet.getRealValue("numBins") ;
  assert(0 != integrand() && integrand()->isValid());
  
  // Allocate coordinate buffer size after number of function dimensions
  _x = new Double_t[_function->getDimension()] ;

  for (UInt_t i=0 ; i<_function->getDimension() ; i++) {
    _xmin.push_back(integrand()->getMinLimit(i));
    _xmax.push_back(integrand()->getMaxLimit(i));
    
    // Retrieve bin configuration from integrand
    list<Double_t>* tmp = integrand()->binBoundaries(i) ;
    if (!tmp) {
      oocoutW((TObject*)0,Integration) << "RooBinIntegrator::RooBinIntegrator WARNING: integrand provide no binning definition observable #" 
				     << i << " substituting default binning of " << _numBins << " bins" << endl ;
      tmp = new list<Double_t> ;
      for (Int_t j=0 ; j<=_numBins ; j++) {
	tmp->push_back(_xmin[i]+j*(_xmax[i]-_xmin[i])/_numBins) ;
      }
    }
    _binb.push_back(tmp) ;
  }

  checkLimits();
} 


//_____________________________________________________________________________
RooAbsIntegrator* RooBinIntegrator::clone(const RooAbsFunc& function, const RooNumIntConfig& config) const
{
  // Clone integrator with new function binding and configuration. Needed by RooNumIntFactory
  return new RooBinIntegrator(function,config) ;
}





//_____________________________________________________________________________
RooBinIntegrator::~RooBinIntegrator()
{
  // Destructor
  if(_x) delete[] _x;
  for (vector<list<Double_t>*>::iterator iter = _binb.begin() ; iter!=_binb.end() ; ++iter) {
    delete (*iter) ;
  }

}


//_____________________________________________________________________________
Bool_t RooBinIntegrator::setLimits(Double_t *xmin, Double_t *xmax) 
{
  // Change our integration limits. Return kTRUE if the new limits are
  // ok, or otherwise kFALSE. Always returns kFALSE and does nothing
  // if this object was constructed to always use our integrand's limits.

  if(_useIntegrandLimits) {
    oocoutE((TObject*)0,Integration) << "RooBinIntegrator::setLimits: cannot override integrand's limits" << endl;
    return kFALSE;
  }
  _xmin[0]= *xmin;
  _xmax[0]= *xmax;
  return checkLimits();
}


//_____________________________________________________________________________
Bool_t RooBinIntegrator::checkLimits() const 
{
  // Check that our integration range is finite and otherwise return kFALSE.
  // Update the limits from the integrand if requested.

  if(_useIntegrandLimits) {
    assert(0 != integrand() && integrand()->isValid());
    _xmin.resize(_function->getDimension()) ;
    _xmax.resize(_function->getDimension()) ;
    for (UInt_t i=0 ; i<_function->getDimension() ; i++) {
      _xmin[i]= integrand()->getMinLimit(i);
      _xmax[i]= integrand()->getMaxLimit(i);
    }
  }
  for (UInt_t i=0 ; i<_function->getDimension() ; i++) {
    if (_xmax[i]<=_xmin[i]) {
      oocoutE((TObject*)0,Integration) << "RooBinIntegrator::checkLimits: bad range with min >= max (_xmin = " << _xmin[i] << " _xmax = " << _xmax[i] << ")" << endl;
      return kFALSE;
    }
    if (RooNumber::isInfinite(_xmin[i]) || RooNumber::isInfinite(_xmax[i])) {
      return kFALSE ;
    }
  }
  
  return kTRUE;
}


//_____________________________________________________________________________
Double_t RooBinIntegrator::integral(const Double_t *) 
{
  // Calculate numeric integral at given set of function binding parameters

  assert(isValid());

  double sum = 0. ;

  if (_function->getDimension()==1) {
    list<Double_t>::iterator iter = _binb[0]->begin() ;
    Double_t xlo = *iter ; iter++ ;
    for (; iter!=_binb[0]->end() ; ++iter) {
      Double_t xhi = *iter ;
      Double_t xcenter = (xhi+xlo)/2 ;
      Double_t binInt = integrand(xvec(xcenter))*(xhi-xlo) ;
      sum += binInt ;
      //cout << "RBI::integral over " << _function->getName() << " 1D binInt[" << xcenter << "] = " << binInt << " running sum = " << sum << endl ;
      xlo=xhi ;
    }
  }

  if (_function->getDimension()==2) {

    list<Double_t>::iterator iter1 = _binb[0]->begin() ;

    Double_t x1lo = *iter1 ; iter1++ ;
    for (; iter1!=_binb[0]->end() ; ++iter1) {

      Double_t x1hi = *iter1 ;
      Double_t x1center = (x1hi+x1lo)/2 ;
      
      list<Double_t>::iterator iter2 = _binb[1]->begin() ;
      Double_t x2lo = *iter2 ; iter2++ ;
      for (; iter2!=_binb[1]->end() ; ++iter2) {

	Double_t x2hi = *iter2 ;
	Double_t x2center = (x2hi+x2lo)/2 ;
      	
	Double_t binInt = integrand(xvec(x1center,x2center))*(x1hi-x1lo)*(x2hi-x2lo) ;
	//cout << "RBI::integral 2D binInt[" << x1center << "," << x2center << "] = " << binInt << " binv = " << (x1hi-x1lo) << "*" << (x2hi-x2lo) << endl ;
	sum += binInt ;
	x2lo=x2hi ;
      }
      x1lo=x1hi ;
    }    
  }

  if (_function->getDimension()==3) {

    list<Double_t>::iterator iter1 = _binb[0]->begin() ;

    Double_t x1lo = *iter1 ; iter1++ ;
    for (; iter1!=_binb[0]->end() ; ++iter1) {

      Double_t x1hi = *iter1 ;
      Double_t x1center = (x1hi+x1lo)/2 ;
      
      list<Double_t>::iterator iter2 = _binb[1]->begin() ;
      Double_t x2lo = *iter2 ; iter2++ ;
      for (; iter2!=_binb[1]->end() ; ++iter2) {

	Double_t x2hi = *iter2 ;
	Double_t x2center = (x2hi+x2lo)/2 ;

	list<Double_t>::iterator iter3 = _binb[2]->begin() ;
	Double_t x3lo = *iter3 ; iter3++ ;
	for (; iter3!=_binb[2]->end() ; ++iter3) {

	  Double_t x3hi = *iter3 ;
	  Double_t x3center = (x3hi+x3lo)/2 ;
	  
	  Double_t binInt = integrand(xvec(x1center,x2center,x3center))*(x1hi-x1lo)*(x2hi-x2lo)*(x3hi-x3lo) ;
	  //cout << "RBI::integral 3D binInt[" << x1center << "," << x2center << "," << x3center << "] = " << binInt << endl ;
	  sum += binInt ;
	  
	  x3lo=x3hi ;
	}
	x2lo=x2hi ;
      }
      x1lo=x1hi ;
    }    
  }

  return sum;
}


 RooBinIntegrator.cxx:1
 RooBinIntegrator.cxx:2
 RooBinIntegrator.cxx:3
 RooBinIntegrator.cxx:4
 RooBinIntegrator.cxx:5
 RooBinIntegrator.cxx:6
 RooBinIntegrator.cxx:7
 RooBinIntegrator.cxx:8
 RooBinIntegrator.cxx:9
 RooBinIntegrator.cxx:10
 RooBinIntegrator.cxx:11
 RooBinIntegrator.cxx:12
 RooBinIntegrator.cxx:13
 RooBinIntegrator.cxx:14
 RooBinIntegrator.cxx:15
 RooBinIntegrator.cxx:16
 RooBinIntegrator.cxx:17
 RooBinIntegrator.cxx:18
 RooBinIntegrator.cxx:19
 RooBinIntegrator.cxx:20
 RooBinIntegrator.cxx:21
 RooBinIntegrator.cxx:22
 RooBinIntegrator.cxx:23
 RooBinIntegrator.cxx:24
 RooBinIntegrator.cxx:25
 RooBinIntegrator.cxx:26
 RooBinIntegrator.cxx:27
 RooBinIntegrator.cxx:28
 RooBinIntegrator.cxx:29
 RooBinIntegrator.cxx:30
 RooBinIntegrator.cxx:31
 RooBinIntegrator.cxx:32
 RooBinIntegrator.cxx:33
 RooBinIntegrator.cxx:34
 RooBinIntegrator.cxx:35
 RooBinIntegrator.cxx:36
 RooBinIntegrator.cxx:37
 RooBinIntegrator.cxx:38
 RooBinIntegrator.cxx:39
 RooBinIntegrator.cxx:40
 RooBinIntegrator.cxx:41
 RooBinIntegrator.cxx:42
 RooBinIntegrator.cxx:43
 RooBinIntegrator.cxx:44
 RooBinIntegrator.cxx:45
 RooBinIntegrator.cxx:46
 RooBinIntegrator.cxx:47
 RooBinIntegrator.cxx:48
 RooBinIntegrator.cxx:49
 RooBinIntegrator.cxx:50
 RooBinIntegrator.cxx:51
 RooBinIntegrator.cxx:52
 RooBinIntegrator.cxx:53
 RooBinIntegrator.cxx:54
 RooBinIntegrator.cxx:55
 RooBinIntegrator.cxx:56
 RooBinIntegrator.cxx:57
 RooBinIntegrator.cxx:58
 RooBinIntegrator.cxx:59
 RooBinIntegrator.cxx:60
 RooBinIntegrator.cxx:61
 RooBinIntegrator.cxx:62
 RooBinIntegrator.cxx:63
 RooBinIntegrator.cxx:64
 RooBinIntegrator.cxx:65
 RooBinIntegrator.cxx:66
 RooBinIntegrator.cxx:67
 RooBinIntegrator.cxx:68
 RooBinIntegrator.cxx:69
 RooBinIntegrator.cxx:70
 RooBinIntegrator.cxx:71
 RooBinIntegrator.cxx:72
 RooBinIntegrator.cxx:73
 RooBinIntegrator.cxx:74
 RooBinIntegrator.cxx:75
 RooBinIntegrator.cxx:76
 RooBinIntegrator.cxx:77
 RooBinIntegrator.cxx:78
 RooBinIntegrator.cxx:79
 RooBinIntegrator.cxx:80
 RooBinIntegrator.cxx:81
 RooBinIntegrator.cxx:82
 RooBinIntegrator.cxx:83
 RooBinIntegrator.cxx:84
 RooBinIntegrator.cxx:85
 RooBinIntegrator.cxx:86
 RooBinIntegrator.cxx:87
 RooBinIntegrator.cxx:88
 RooBinIntegrator.cxx:89
 RooBinIntegrator.cxx:90
 RooBinIntegrator.cxx:91
 RooBinIntegrator.cxx:92
 RooBinIntegrator.cxx:93
 RooBinIntegrator.cxx:94
 RooBinIntegrator.cxx:95
 RooBinIntegrator.cxx:96
 RooBinIntegrator.cxx:97
 RooBinIntegrator.cxx:98
 RooBinIntegrator.cxx:99
 RooBinIntegrator.cxx:100
 RooBinIntegrator.cxx:101
 RooBinIntegrator.cxx:102
 RooBinIntegrator.cxx:103
 RooBinIntegrator.cxx:104
 RooBinIntegrator.cxx:105
 RooBinIntegrator.cxx:106
 RooBinIntegrator.cxx:107
 RooBinIntegrator.cxx:108
 RooBinIntegrator.cxx:109
 RooBinIntegrator.cxx:110
 RooBinIntegrator.cxx:111
 RooBinIntegrator.cxx:112
 RooBinIntegrator.cxx:113
 RooBinIntegrator.cxx:114
 RooBinIntegrator.cxx:115
 RooBinIntegrator.cxx:116
 RooBinIntegrator.cxx:117
 RooBinIntegrator.cxx:118
 RooBinIntegrator.cxx:119
 RooBinIntegrator.cxx:120
 RooBinIntegrator.cxx:121
 RooBinIntegrator.cxx:122
 RooBinIntegrator.cxx:123
 RooBinIntegrator.cxx:124
 RooBinIntegrator.cxx:125
 RooBinIntegrator.cxx:126
 RooBinIntegrator.cxx:127
 RooBinIntegrator.cxx:128
 RooBinIntegrator.cxx:129
 RooBinIntegrator.cxx:130
 RooBinIntegrator.cxx:131
 RooBinIntegrator.cxx:132
 RooBinIntegrator.cxx:133
 RooBinIntegrator.cxx:134
 RooBinIntegrator.cxx:135
 RooBinIntegrator.cxx:136
 RooBinIntegrator.cxx:137
 RooBinIntegrator.cxx:138
 RooBinIntegrator.cxx:139
 RooBinIntegrator.cxx:140
 RooBinIntegrator.cxx:141
 RooBinIntegrator.cxx:142
 RooBinIntegrator.cxx:143
 RooBinIntegrator.cxx:144
 RooBinIntegrator.cxx:145
 RooBinIntegrator.cxx:146
 RooBinIntegrator.cxx:147
 RooBinIntegrator.cxx:148
 RooBinIntegrator.cxx:149
 RooBinIntegrator.cxx:150
 RooBinIntegrator.cxx:151
 RooBinIntegrator.cxx:152
 RooBinIntegrator.cxx:153
 RooBinIntegrator.cxx:154
 RooBinIntegrator.cxx:155
 RooBinIntegrator.cxx:156
 RooBinIntegrator.cxx:157
 RooBinIntegrator.cxx:158
 RooBinIntegrator.cxx:159
 RooBinIntegrator.cxx:160
 RooBinIntegrator.cxx:161
 RooBinIntegrator.cxx:162
 RooBinIntegrator.cxx:163
 RooBinIntegrator.cxx:164
 RooBinIntegrator.cxx:165
 RooBinIntegrator.cxx:166
 RooBinIntegrator.cxx:167
 RooBinIntegrator.cxx:168
 RooBinIntegrator.cxx:169
 RooBinIntegrator.cxx:170
 RooBinIntegrator.cxx:171
 RooBinIntegrator.cxx:172
 RooBinIntegrator.cxx:173
 RooBinIntegrator.cxx:174
 RooBinIntegrator.cxx:175
 RooBinIntegrator.cxx:176
 RooBinIntegrator.cxx:177
 RooBinIntegrator.cxx:178
 RooBinIntegrator.cxx:179
 RooBinIntegrator.cxx:180
 RooBinIntegrator.cxx:181
 RooBinIntegrator.cxx:182
 RooBinIntegrator.cxx:183
 RooBinIntegrator.cxx:184
 RooBinIntegrator.cxx:185
 RooBinIntegrator.cxx:186
 RooBinIntegrator.cxx:187
 RooBinIntegrator.cxx:188
 RooBinIntegrator.cxx:189
 RooBinIntegrator.cxx:190
 RooBinIntegrator.cxx:191
 RooBinIntegrator.cxx:192
 RooBinIntegrator.cxx:193
 RooBinIntegrator.cxx:194
 RooBinIntegrator.cxx:195
 RooBinIntegrator.cxx:196
 RooBinIntegrator.cxx:197
 RooBinIntegrator.cxx:198
 RooBinIntegrator.cxx:199
 RooBinIntegrator.cxx:200
 RooBinIntegrator.cxx:201
 RooBinIntegrator.cxx:202
 RooBinIntegrator.cxx:203
 RooBinIntegrator.cxx:204
 RooBinIntegrator.cxx:205
 RooBinIntegrator.cxx:206
 RooBinIntegrator.cxx:207
 RooBinIntegrator.cxx:208
 RooBinIntegrator.cxx:209
 RooBinIntegrator.cxx:210
 RooBinIntegrator.cxx:211
 RooBinIntegrator.cxx:212
 RooBinIntegrator.cxx:213
 RooBinIntegrator.cxx:214
 RooBinIntegrator.cxx:215
 RooBinIntegrator.cxx:216
 RooBinIntegrator.cxx:217
 RooBinIntegrator.cxx:218
 RooBinIntegrator.cxx:219
 RooBinIntegrator.cxx:220
 RooBinIntegrator.cxx:221
 RooBinIntegrator.cxx:222
 RooBinIntegrator.cxx:223
 RooBinIntegrator.cxx:224
 RooBinIntegrator.cxx:225
 RooBinIntegrator.cxx:226
 RooBinIntegrator.cxx:227
 RooBinIntegrator.cxx:228
 RooBinIntegrator.cxx:229
 RooBinIntegrator.cxx:230
 RooBinIntegrator.cxx:231
 RooBinIntegrator.cxx:232
 RooBinIntegrator.cxx:233
 RooBinIntegrator.cxx:234
 RooBinIntegrator.cxx:235
 RooBinIntegrator.cxx:236
 RooBinIntegrator.cxx:237
 RooBinIntegrator.cxx:238
 RooBinIntegrator.cxx:239
 RooBinIntegrator.cxx:240
 RooBinIntegrator.cxx:241
 RooBinIntegrator.cxx:242
 RooBinIntegrator.cxx:243
 RooBinIntegrator.cxx:244
 RooBinIntegrator.cxx:245
 RooBinIntegrator.cxx:246
 RooBinIntegrator.cxx:247
 RooBinIntegrator.cxx:248
 RooBinIntegrator.cxx:249
 RooBinIntegrator.cxx:250
 RooBinIntegrator.cxx:251
 RooBinIntegrator.cxx:252
 RooBinIntegrator.cxx:253
 RooBinIntegrator.cxx:254
 RooBinIntegrator.cxx:255
 RooBinIntegrator.cxx:256
 RooBinIntegrator.cxx:257
 RooBinIntegrator.cxx:258
 RooBinIntegrator.cxx:259
 RooBinIntegrator.cxx:260
 RooBinIntegrator.cxx:261
 RooBinIntegrator.cxx:262
 RooBinIntegrator.cxx:263
 RooBinIntegrator.cxx:264
 RooBinIntegrator.cxx:265
 RooBinIntegrator.cxx:266
 RooBinIntegrator.cxx:267
 RooBinIntegrator.cxx:268
 RooBinIntegrator.cxx:269
 RooBinIntegrator.cxx:270
 RooBinIntegrator.cxx:271
 RooBinIntegrator.cxx:272
 RooBinIntegrator.cxx:273
 RooBinIntegrator.cxx:274
 RooBinIntegrator.cxx:275
 RooBinIntegrator.cxx:276
 RooBinIntegrator.cxx:277
 RooBinIntegrator.cxx:278
 RooBinIntegrator.cxx:279
 RooBinIntegrator.cxx:280
 RooBinIntegrator.cxx:281
 RooBinIntegrator.cxx:282
 RooBinIntegrator.cxx:283
 RooBinIntegrator.cxx:284
 RooBinIntegrator.cxx:285
 RooBinIntegrator.cxx:286
 RooBinIntegrator.cxx:287
 RooBinIntegrator.cxx:288
 RooBinIntegrator.cxx:289
 RooBinIntegrator.cxx:290
 RooBinIntegrator.cxx:291
 RooBinIntegrator.cxx:292
 RooBinIntegrator.cxx:293
 RooBinIntegrator.cxx:294
 RooBinIntegrator.cxx:295
 RooBinIntegrator.cxx:296
 RooBinIntegrator.cxx:297
 RooBinIntegrator.cxx:298