/*****************************************************************************
 * Project: RooFit                                                           *
 * Package: RooFitCore                                                       *
 * @(#)root/roofit:$Id: RooHistPdf.cxx 28259 2009-04-16 16: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
// RooHistPdf implements a probablity density function sampled from a 
// multidimensional histogram. The histogram distribution is explicitly
// normalized by RooHistPdf and can have an arbitrary number of real or 
// discrete dimensions.
// END_HTML
//

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

#include "RooHistPdf.h"
#include "RooDataHist.h"
#include "RooMsgService.h"
#include "RooRealVar.h"
#include "RooCategory.h"



ClassImp(RooHistPdf)
;



//_____________________________________________________________________________
RooHistPdf::RooHistPdf() : _dataHist(0), _totVolume(0)
{
  // Default constructor
  _histObsIter = _histObsList.createIterator() ;
  _pdfObsIter = _pdfObsList.createIterator() ;
}


//_____________________________________________________________________________
RooHistPdf::RooHistPdf(const char *name, const char *title, const RooArgSet& vars, 
		       const RooDataHist& dhist, Int_t intOrder) :
  RooAbsPdf(name,title), 
  _histObsList("histObs","List of histogram observables",this),
  _pdfObsList("pdfObs","List of p.d.f. observables",this),
  _dataHist((RooDataHist*)&dhist), 
  _codeReg(10),
  _intOrder(intOrder),
  _cdfBoundaries(kFALSE),
  _totVolume(0),
  _unitNorm(kFALSE)
{
  // Constructor from a RooDataHist. RooDataHist dimensions
  // can be either real or discrete. See RooDataHist::RooDataHist for details on the binning.
  // RooHistPdf neither owns or clone 'dhist' and the user must ensure the input histogram exists
  // for the entire life span of this PDF.

  _histObsList.add(vars) ;

  // Verify that vars and dhist.get() have identical contents
  const RooArgSet* dvars = dhist.get() ;
  if (vars.getSize()!=dvars->getSize()) {
    coutE(InputArguments) << "RooHistPdf::ctor(" << GetName() 
			  << ") ERROR variable list and RooDataHist must contain the same variables." << endl ;
    assert(0) ;
  }
  TIterator* iter = vars.createIterator() ;
  RooAbsArg* arg ;
  while((arg=(RooAbsArg*)iter->Next())) {
    if (!dvars->find(arg->GetName())) {
      coutE(InputArguments) << "RooHistPdf::ctor(" << GetName() 
			    << ") ERROR variable list and RooDataHist must contain the same variables." << endl ;
      assert(0) ;
    }
  }
  delete iter ;

  _histObsIter = _histObsList.createIterator() ;
  _pdfObsIter = _pdfObsList.createIterator() ;
}




//_____________________________________________________________________________
RooHistPdf::RooHistPdf(const char *name, const char *title, const RooArgList& pdfObs, 
		       const RooArgList& histObs, const RooDataHist& dhist, Int_t intOrder) :
  RooAbsPdf(name,title), 
  _histObsList("histObs","List of histogram observables",this,kFALSE,kFALSE),
  _pdfObsList("pdfObs","List of p.d.f. observables",this),
  _dataHist((RooDataHist*)&dhist), 
  _codeReg(10),
  _intOrder(intOrder),
  _cdfBoundaries(kFALSE),
  _totVolume(0),
  _unitNorm(kFALSE)
{
  // Constructor from a RooDataHist. The first list of observables are the p.d.f.
  // observables, which may any RooAbsReal (function or variable). The second list
  // are the corresponding observables in the RooDataHist which must be of type
  // RooRealVar or RooCategory This constructor thus allows to apply a coordinate transformation
  // on the histogram data to be applied.

  _histObsList.add(histObs) ;
  _pdfObsList.add(pdfObs) ;

  // Verify that vars and dhist.get() have identical contents
  const RooArgSet* dvars = dhist.get() ;
  if (histObs.getSize()!=dvars->getSize()) {
    coutE(InputArguments) << "RooHistPdf::ctor(" << GetName() 
			  << ") ERROR histogram variable list and RooDataHist must contain the same variables." << endl ;
    throw(string("RooHistPdf::ctor() ERROR: histogram variable list and RooDataHist must contain the same variables")) ;
  }
  TIterator* iter = histObs.createIterator() ;
  RooAbsArg* arg ;
  while((arg=(RooAbsArg*)iter->Next())) {
    if (!dvars->find(arg->GetName())) {
      coutE(InputArguments) << "RooHistPdf::ctor(" << GetName() 
			    << ") ERROR variable list and RooDataHist must contain the same variables." << endl ;
      throw(string("RooHistPdf::ctor() ERROR: histogram variable list and RooDataHist must contain the same variables")) ;
    }
    if (!arg->isFundamental()) {
      coutE(InputArguments) << "RooHistPdf::ctor(" << GetName() 
			    << ") ERROR all elements of histogram observables set must be of type RooRealVar or RooCategory." << endl ;
      throw(string("RooHistPdf::ctor() ERROR all elements of histogram observables set must be of type RooRealVar or RooCategory.")) ;
    }
  }
  delete iter ;

  _histObsIter = _histObsList.createIterator() ;
  _pdfObsIter = _pdfObsList.createIterator() ;
}



//_____________________________________________________________________________
RooHistPdf::RooHistPdf(const RooHistPdf& other, const char* name) :
  RooAbsPdf(other,name), 
  _histObsList("histObs",this,other._histObsList),
  _pdfObsList("pdfObs",this,other._pdfObsList),
  _dataHist(other._dataHist),
  _codeReg(other._codeReg),
  _intOrder(other._intOrder),
  _cdfBoundaries(other._cdfBoundaries),
  _totVolume(other._totVolume),
  _unitNorm(other._unitNorm)
{
  // Copy constructor

  _histObsIter = _histObsList.createIterator() ;
  _pdfObsIter = _pdfObsList.createIterator() ;
}




//_____________________________________________________________________________
RooHistPdf::~RooHistPdf()
{
  // Destructor

  delete _histObsIter ;
  delete _pdfObsIter ;
}





//_____________________________________________________________________________
Double_t RooHistPdf::evaluate() const
{
  // Return the current value: The value of the bin enclosing the current coordinates
  // of the observables, normalized by the histograms contents. Interpolation
  // is applied if the RooHistPdf is configured to do that
  
  // Transfer values from   
  if (_pdfObsList.getSize()>0) {
    _histObsIter->Reset() ;
    _pdfObsIter->Reset() ;
    RooAbsArg* harg, *parg ;
    while((harg=(RooAbsArg*)_histObsIter->Next())) {
      parg = (RooAbsArg*)_pdfObsIter->Next() ;
      if (harg != parg) {
	parg->syncCache() ;
	harg->copyCache(parg,kTRUE) ;
      }
    }
  }

  Double_t ret =  _dataHist->weight(_histObsList,_intOrder,kFALSE,_cdfBoundaries) ;  
  if (ret<0) {
    ret=0 ;
  }
  return ret ;
}


//_____________________________________________________________________________
Double_t RooHistPdf::totVolume() const
{
  // Return the total volume spanned by the observables of the RooHistPdf

  // Return previously calculated value, if any
  if (_totVolume>0) {
    return _totVolume ;
  }
  _totVolume = 1. ;
  TIterator* iter = _histObsList.createIterator() ;
  RooAbsArg* arg ;
  while((arg=(RooAbsArg*)iter->Next())) {
    RooRealVar* real = dynamic_cast<RooRealVar*>(arg) ;
    if (real) {
      _totVolume *= (real->getMax()-real->getMin()) ;
    } else {
      RooCategory* cat = dynamic_cast<RooCategory*>(arg) ;
      if (cat) {
	_totVolume *= cat->numTypes() ;
      }
    }
  }
  delete iter ;
  return _totVolume ;
}



//_____________________________________________________________________________
Int_t RooHistPdf::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName) const 
{
  // Determine integration scenario. If no interpolation is used,
  // RooHistPdf can perform all integrals over its dependents
  // analytically via partial or complete summation of the input
  // histogram. If interpolation is used on the integral over
  // all histogram observables is supported

  // Only analytical integrals over the full range are defined
  if (rangeName!=0) {
    return 0 ;
  }

  // Simplest scenario, integrate over all dependents
  RooAbsCollection *allVarsCommon = allVars.selectCommon(_histObsList) ;  
  Bool_t intAllObs = (allVarsCommon->getSize()==_histObsList.getSize()) ;
  delete allVarsCommon ;
  if (intAllObs && matchArgs(allVars,analVars,_histObsList)) {
    return 1000 ;
  }

  // Disable partial analytical integrals if interpolation is used
  if (_intOrder>0) {
    return 0 ;
  }

  // Find subset of _histObsList that integration is requested over
  RooArgSet* allVarsSel = (RooArgSet*) allVars.selectCommon(_histObsList) ;
  if (allVarsSel->getSize()==0) {
    delete allVarsSel ;
    return 0 ;
  }

  // Partial integration scenarios.
  // Build unique code from bit mask of integrated variables in depList
  Int_t code(0),n(0) ;
  TIterator* iter = _histObsList.createIterator() ;
  RooAbsArg* arg ;
  while((arg=(RooAbsArg*)iter->Next())) {
    if (allVars.find(arg->GetName())) code |= (1<<n) ;
    n++ ;
  }
  delete iter ;
  analVars.add(*allVarsSel) ;

  return code ;

}



//_____________________________________________________________________________
Double_t RooHistPdf::analyticalIntegral(Int_t code, const char* /*rangeName*/) const 
{
  // Return integral identified by 'code'. The actual integration
  // is deferred to RooDataHist::sum() which implements partial
  // or complete summation over the histograms contents

  // WVE needs adaptation for rangeName feature
  // Simplest scenario, integration over all dependents
  if (code==1000) {
    return _dataHist->sum(kTRUE) ;
  }

  // Partial integration scenario, retrieve set of variables, calculate partial sum

  RooArgSet intSet ;
  TIterator* iter = _histObsList.createIterator() ;
  RooAbsArg* arg ;
  Int_t n(0) ;
  while((arg=(RooAbsArg*)iter->Next())) {
    if (code & (1<<n)) {
      intSet.add(*arg) ;
    }
    n++ ;
  }
  delete iter ;

  Double_t ret =  _dataHist->sum(intSet,_histObsList,kTRUE) ;
  return ret ;
}



//_____________________________________________________________________________
list<Double_t>* RooHistPdf::plotSamplingHint(RooAbsRealLValue& obs, Double_t xlo, Double_t xhi) const
{
  // Return sampling hint for making curves of (projections) of this function
  // as the recursive division strategy of RooCurve cannot deal efficiently
  // with the vertical lines that occur in a non-interpolated histogram

  // No hints are required when interpolation is used
  if (_intOrder>0) {
    return 0 ;
  }

  // Check that observable is in dataset, if not no hint is generated
  RooAbsLValue* lvarg = dynamic_cast<RooAbsLValue*>(_dataHist->get()->find(obs.GetName())) ;
  if (!lvarg) {
    return 0 ;
  }

  // Retrieve position of all bin boundaries
  const RooAbsBinning* binning = lvarg->getBinningPtr(0) ;
  Double_t* boundaries = binning->array() ;

  list<Double_t>* hint = new list<Double_t> ;

  // Widen range slighty
  xlo = xlo - 0.01*(xhi-xlo) ;
  xhi = xhi + 0.01*(xhi-xlo) ;

  Double_t delta = (xhi-xlo)*1e-8 ;
 
  // Construct array with pairs of points positioned epsilon to the left and
  // right of the bin boundaries
  for (Int_t i=0 ; i<binning->numBoundaries() ; i++) {
    if (boundaries[i]>=xlo && boundaries[i]<=xhi) {
      hint->push_back(boundaries[i]-delta) ;
      hint->push_back(boundaries[i]+delta) ;
    }
  }

  return hint ;
}



//_____________________________________________________________________________
Int_t RooHistPdf::getMaxVal(const RooArgSet& /*vars*/) const 
{
  // Not implemented yet
  return 0 ;
}


//_____________________________________________________________________________
Double_t RooHistPdf::maxVal(Int_t /*code*/) const 
{
  // Not implemented yet
  return 0 ;
}