// 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"
#include "RooWorkspace.h"
#include "RooGlobalFunc.h"
#include "TError.h"
using namespace std;
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() ;
RooFIter oiter = _histObsList.fwdIterator() ;
RooAbsArg* hobs ;
while ((hobs = oiter.next())) {
RooAbsArg* dhobs = dhist.get()->find(hobs->GetName()) ;
RooRealVar* dhreal = dynamic_cast<RooRealVar*>(dhobs) ;
if (dhreal){
((RooRealVar*)hobs)->setRange(dhreal->getMin(),dhreal->getMax()) ;
}
}
}
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() ;
RooFIter oiter = _histObsList.fwdIterator() ;
RooAbsArg* hobs ;
while ((hobs = oiter.next())) {
RooAbsArg* dhobs = dhist.get()->find(hobs->GetName()) ;
RooRealVar* dhreal = dynamic_cast<RooRealVar*>(dhobs) ;
if (dhreal){
((RooRealVar*)hobs)->setRange(dhreal->getMin(),dhreal->getMax()) ;
}
}
}
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) ;
if (!harg->inRange(0)) {
return 0 ;
}
}
}
}
Double_t ret = _dataHist->weight(_histObsList,_intOrder,_unitNorm?kFALSE:kTRUE,_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 ;
}
namespace {
bool fullRange(const RooAbsArg& x, const RooAbsArg& y ,const char* range)
{
const RooAbsRealLValue *_x = dynamic_cast<const RooAbsRealLValue*>(&x);
const RooAbsRealLValue *_y = dynamic_cast<const RooAbsRealLValue*>(&y);
if (!_x || !_y) return false;
if (!range || !strlen(range) || !_x->hasRange(range) ||
_x->getBinningPtr(range)->isParameterized()) {
if (range && strlen(range) && _x->getBinningPtr(range)->isParameterized())
return false;
return (_x->getMin() == _y->getMin() && _x->getMax() == _y->getMax());
}
return (_x->getMin(range) == _y->getMin() && _x->getMax(range) == _y->getMax());
}
}
Int_t RooHistPdf::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName) const
{
RooFIter it = _pdfObsList.fwdIterator();
RooFIter jt = _histObsList.fwdIterator();
Int_t code = 0, frcode = 0, n = 0;
for (RooAbsArg *pa = 0, *ha = 0; (pa = it.next()) && (ha = jt.next()); ++n) {
if (allVars.find(*pa)) {
code |= 2 << n;
analVars.add(*pa);
if (fullRange(*pa, *ha, rangeName)) {
frcode |= 2 << n;
}
}
}
if (code == frcode) {
code |= 1;
}
if (_intOrder > 1 && !(code & 1)) {
analVars.removeAll();
return 0;
}
return (code >= 2) ? code : 0;
}
Double_t RooHistPdf::analyticalIntegral(Int_t code, const char* rangeName) const
{
if (((2 << _histObsList.getSize()) - 1) == code) {
return _dataHist->sum(kFALSE);
}
RooArgSet intSet;
std::map<const RooAbsArg*, std::pair<Double_t, Double_t> > ranges;
RooFIter it = _pdfObsList.fwdIterator();
RooFIter jt = _histObsList.fwdIterator();
Int_t n(0);
for (RooAbsArg *pa = 0, *ha = 0; (pa = it.next()) && (ha = jt.next()); ++n) {
if (code & (2 << n)) {
intSet.add(*ha);
}
if (!(code & 1)) {
RooAbsRealLValue* rlv = dynamic_cast<RooAbsRealLValue*>(pa);
if (rlv) {
const RooAbsBinning* binning = rlv->getBinningPtr(rangeName);
if (rangeName && rlv->hasRange(rangeName)) {
ranges[ha] = std::make_pair(
rlv->getMin(rangeName), rlv->getMax(rangeName));
} else if (binning) {
if (!binning->isParameterized()) {
ranges[ha] = std::make_pair(
binning->lowBound(), binning->highBound());
} else {
ranges[ha] = std::make_pair(
binning->lowBoundFunc()->getVal(), binning->highBoundFunc()->getVal());
}
}
}
}
if (ha != pa) {
pa->syncCache();
ha->copyCache(pa,kTRUE);
}
}
Double_t ret = (code & 1) ?
_dataHist->sum(intSet,_histObsList,kTRUE,kTRUE) :
_dataHist->sum(intSet,_histObsList,kFALSE,kTRUE, ranges);
return ret ;
}
list<Double_t>* RooHistPdf::plotSamplingHint(RooAbsRealLValue& obs, Double_t xlo, Double_t xhi) const
{
if (_intOrder>0) {
return 0 ;
}
_histObsIter->Reset() ;
_pdfObsIter->Reset() ;
RooAbsArg *pdfObs, *histObs, *dhObs(0) ;
while ((pdfObs = (RooAbsArg*)_pdfObsIter->Next()) && !dhObs) {
histObs = (RooAbsArg*) _histObsIter->Next() ;
if (TString(obs.GetName())==pdfObs->GetName()) {
dhObs = _dataHist->get()->find(histObs->GetName()) ;
}
}
if (!dhObs) {
return 0 ;
}
RooAbsLValue* lval = dynamic_cast<RooAbsLValue*>(dhObs) ;
if (!lval) {
return 0 ;
}
const RooAbsBinning* binning = lval->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 ;
}
std::list<Double_t>* RooHistPdf::binBoundaries(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> ;
for (Int_t i=0 ; i<binning->numBoundaries() ; i++) {
if (boundaries[i]>=xlo && boundaries[i]<=xhi) {
hint->push_back(boundaries[i]) ;
}
}
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
{
R__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 ;
}
Bool_t RooHistPdf::areIdentical(const RooDataHist& dh1, const RooDataHist& dh2)
{
if (fabs(dh1.sumEntries()-dh2.sumEntries())>1e-8) return kFALSE ;
if (dh1.numEntries() != dh2.numEntries()) return kFALSE ;
for (int i=0 ; i < dh1.numEntries() ; i++) {
dh1.get(i) ;
dh2.get(i) ;
if (fabs(dh1.weight()-dh2.weight())>1e-8) return kFALSE ;
}
return kTRUE ;
}
Bool_t RooHistPdf::importWorkspaceHook(RooWorkspace& ws)
{
std::list<RooAbsData*> allData = ws.allData() ;
std::list<RooAbsData*>::const_iterator iter ;
for (iter = allData.begin() ; iter != allData.end() ; ++iter) {
if (*iter == _dataHist) {
return kFALSE ;
}
}
RooAbsData* wsdata = ws.embeddedData(_dataHist->GetName()) ;
if (wsdata) {
if (wsdata->InheritsFrom(RooDataHist::Class())) {
if (areIdentical((RooDataHist&)*wsdata,*_dataHist)) {
_dataHist = (RooDataHist*) wsdata ;
} else {
TString uniqueName = Form("%s_%s",_dataHist->GetName(),GetName()) ;
Bool_t flag = ws.import(*_dataHist,RooFit::Rename(uniqueName.Data()),RooFit::Embedded()) ;
if (flag) {
coutE(ObjectHandling) << " RooHistPdf::importWorkspaceHook(" << GetName() << ") unable to import clone of underlying RooDataHist with unique name " << uniqueName << ", abort" << endl ;
return kTRUE ;
}
_dataHist = (RooDataHist*) ws.embeddedData(uniqueName.Data()) ;
}
} else {
TString uniqueName = Form("%s_%s",_dataHist->GetName(),GetName()) ;
Bool_t flag = ws.import(*_dataHist,RooFit::Rename(uniqueName.Data()),RooFit::Embedded()) ;
if (flag) {
coutE(ObjectHandling) << " RooHistPdf::importWorkspaceHook(" << GetName() << ") unable to import clone of underlying RooDataHist with unique name " << uniqueName << ", abort" << endl ;
return kTRUE ;
}
_dataHist = (RooDataHist*) ws.embeddedData(uniqueName.Data()) ;
}
return kFALSE ;
}
ws.import(*_dataHist,RooFit::Embedded()) ;
_dataHist = (RooDataHist*) ws.embeddedData(_dataHist->GetName()) ;
return kFALSE ;
}
void RooHistPdf::Streamer(TBuffer &R__b)
{
if (R__b.IsReading()) {
R__b.ReadClassBuffer(RooHistPdf::Class(),this);
} else {
R__b.WriteClassBuffer(RooHistPdf::Class(),this);
}
}