// 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), _unitNorm(kFALSE)
{
_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),
_pdfObsList("pdfObs","List of p.d.f. observables",this),
_dataHist((RooDataHist*)&dhist),
_codeReg(10),
_intOrder(intOrder),
_cdfBoundaries(kFALSE),
_totVolume(0),
_unitNorm(kFALSE)
{
_histObsList.addClone(vars) ;
_pdfObsList.add(vars) ;
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),
_pdfObsList("pdfObs","List of p.d.f. observables",this),
_dataHist((RooDataHist*)&dhist),
_codeReg(10),
_intOrder(intOrder),
_cdfBoundaries(kFALSE),
_totVolume(0),
_unitNorm(kFALSE)
{
_histObsList.addClone(histObs) ;
_pdfObsList.add(pdfObs) ;
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),
_pdfObsList("pdfObs",this,other._pdfObsList),
_dataHist(other._dataHist),
_codeReg(other._codeReg),
_intOrder(other._intOrder),
_cdfBoundaries(other._cdfBoundaries),
_totVolume(other._totVolume),
_unitNorm(other._unitNorm)
{
_histObsList.addClone(other._histObsList) ;
_histObsIter = _histObsList.createIterator() ;
_pdfObsIter = _pdfObsList.createIterator() ;
}
RooHistPdf::~RooHistPdf()
{
delete _histObsIter ;
delete _pdfObsIter ;
}
Double_t RooHistPdf::evaluate() const
{
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
{
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
{
if (rangeName!=0) {
return 0 ;
}
RooArgList hobsl(_histObsList),pobsl(_pdfObsList) ;
RooArgSet allVarsHist ;
TIterator* iter = allVars.createIterator() ;
RooAbsArg* pdfobs ;
while((pdfobs=(RooAbsArg*)iter->Next())) {
Int_t idx = pobsl.index(pdfobs) ;
if (idx>=0) {
RooAbsArg* hobs = hobsl.at(idx) ;
if (hobs) {
allVarsHist.add(*hobs) ;
}
}
}
delete iter ;
RooAbsCollection *allVarsCommon = allVarsHist.selectCommon(_histObsList) ;
Bool_t intAllObs = (allVarsCommon->getSize()==_histObsList.getSize()) ;
delete allVarsCommon ;
if (intAllObs) {
analVars.add(allVars) ;
return 1000 ;
}
RooArgSet* allVarsSel = (RooArgSet*) allVarsHist.selectCommon(_histObsList) ;
if (allVarsSel->getSize()==0) {
delete allVarsSel ;
return 0 ;
}
Int_t code(0),n(0) ;
iter = _histObsList.createIterator() ;
RooAbsArg* arg ;
while((arg=(RooAbsArg*)iter->Next())) {
if (allVars.find(arg->GetName())) {
code |= (1<<n) ;
analVars.add(*pobsl.at(n)) ;
}
n++ ;
}
delete iter ;
return code ;
}
Double_t RooHistPdf::analyticalIntegral(Int_t code, const char* ) const
{
if (code==1000) {
return _dataHist->sum(kTRUE) ;
}
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 ;
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->sum(intSet,_histObsList,kTRUE) ;
return ret ;
}
list<Double_t>* RooHistPdf::plotSamplingHint(RooAbsRealLValue& obs, Double_t xlo, Double_t xhi) const
{
if (_intOrder>0) {
return 0 ;
}
RooAbsLValue* lvarg = dynamic_cast<RooAbsLValue*>(_dataHist->get()->find(obs.GetName())) ;
if (!lvarg) {
return 0 ;
}
const RooAbsBinning* binning = lvarg->getBinningPtr(0) ;
Double_t* boundaries = binning->array() ;
list<Double_t>* hint = new list<Double_t> ;
xlo = xlo - 0.01*(xhi-xlo) ;
xhi = xhi + 0.01*(xhi-xlo) ;
Double_t delta = (xhi-xlo)*1e-8 ;
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
{
RooAbsCollection* common = _pdfObsList.selectCommon(vars) ;
if (common->getSize()==_pdfObsList.getSize()) {
delete common ;
return 1;
}
delete common ;
return 0 ;
}
Double_t RooHistPdf::maxVal(Int_t code) const
{
assert(code==1) ;
Double_t max(-1) ;
for (Int_t i=0 ; i<_dataHist->numEntries() ; i++) {
_dataHist->get(i) ;
Double_t wgt = _dataHist->weight() ;
if (wgt>max) max=wgt ;
}
return max*1.05 ;
}