#include "RooFit.h"
#include "TH1.h"
#include "TH1.h"
#include "TMath.h"
#include "RooDataHist.h"
#include "RooAbsLValue.h"
#include "RooArgList.h"
#include "RooRealVar.h"
#include "RooMath.h"
#include "RooBinning.h"
#include "RooPlot.h"
#include "RooHistError.h"
ClassImp(RooDataHist)
;
RooDataHist::RooDataHist()
{
_arrSize = 0 ;
_wgt = 0 ;
_errLo = 0 ;
_errHi = 0 ;
_sumw2 = 0 ;
_binv = 0 ;
_pbinv = 0 ;
_idxMult = 0 ;
_curWeight = 0 ;
_curIndex = -1 ;
_realIter = 0 ;
}
RooDataHist::RooDataHist(const char *name, const char *title, const RooArgSet& vars) :
RooTreeData(name,title,vars), _curWeight(0), _curVolume(1), _pbinv(0)
{
initialize() ;
appendToDir(this,kTRUE) ;
}
RooDataHist::RooDataHist(const char *name, const char *title, const RooArgSet& vars, const RooAbsData& data, Double_t weight) :
RooTreeData(name,title,vars), _curWeight(0), _curVolume(1), _pbinv(0)
{
initialize() ;
add(data,(const RooFormulaVar*)0,weight) ;
appendToDir(this,kTRUE) ;
}
RooDataHist::RooDataHist(const char *name, const char *title, const RooArgList& vars, const TH1* hist, Double_t weight) :
RooTreeData(name,title,vars), _curWeight(0), _curVolume(1), _pbinv(0)
{
if (vars.getSize() != hist->GetDimension()) {
cout << "RooDataHist::ctor(" << GetName() << ") ERROR: dimension of input histogram must match "
<< "number of dimension variables" << endl ;
assert(0) ;
}
TH1* histo = const_cast<TH1*>(hist) ;
RooRealVar* xvar = (RooRealVar*) _vars.find(vars.at(0)->GetName()) ;
if (!dynamic_cast<RooRealVar*>(xvar)) {
cout << "RooDataHist::ctor(" << GetName() << ") ERROR: dimension " << xvar->GetName() << " must be real" << endl ;
assert(0) ;
}
Double_t xlo = ((RooRealVar*)vars.at(0))->getMin() ;
Double_t xhi = ((RooRealVar*)vars.at(0))->getMax() ;
Int_t xmin(0) ;
if (histo->GetXaxis()->GetXbins()->GetArray()) {
RooBinning xbins(histo->GetNbinsX(),histo->GetXaxis()->GetXbins()->GetArray()) ;
Double_t xloAdj = xbins.binLow(xbins.binNumber(xlo+1e-6)) ;
Double_t xhiAdj = xbins.binHigh(xbins.binNumber(xhi-1e-6)) ;
xbins.setRange(xloAdj,xhiAdj) ;
((RooRealVar*)vars.at(0))->setBinning(xbins) ;
if (fabs(xloAdj-xlo)>1e-6||fabs(xhiAdj-xhi)) {
cout << "RooDataHist::ctor(" << GetName() << "): fit range of variable " << xvar->GetName() << " expanded to nearest bin boundaries: ["
<< xlo << "," << xhi << "] --> [" << xloAdj << "," << xhiAdj << "]" << endl ;
}
xvar->setBinning(xbins) ;
xmin = xbins.rawBinNumber(xloAdj+1e-6) ;
} else {
RooBinning xbins(histo->GetXaxis()->GetXmin(),histo->GetXaxis()->GetXmax()) ;
xbins.addUniform(histo->GetNbinsX(),histo->GetXaxis()->GetXmin(),histo->GetXaxis()->GetXmax()) ;
Double_t xloAdj = xbins.binLow(xbins.binNumber(xlo+1e-6)) ;
Double_t xhiAdj = xbins.binHigh(xbins.binNumber(xhi-1e-6)) ;
xbins.setRange(xloAdj,xhiAdj) ;
((RooRealVar*)vars.at(0))->setRange(xloAdj,xhiAdj) ;
if (fabs(xloAdj-xlo)>1e-6||fabs(xhiAdj-xhi)) {
cout << "RooDataHist::ctor(" << GetName() << "): fit range of variable " << xvar->GetName() << " expanded to nearest bin boundaries: ["
<< xlo << "," << xhi << "] --> [" << xloAdj << "," << xhiAdj << "]" << endl ;
}
RooUniformBinning xbins2(xloAdj,xhiAdj,xbins.numBins()) ;
xvar->setBinning(xbins2) ;
xmin = xbins.rawBinNumber(xloAdj+1e-6) ;
}
RooRealVar* yvar = (RooRealVar*) (vars.at(1) ? _vars.find(vars.at(1)->GetName()) : 0 ) ;
Int_t ymin(0) ;
if (yvar) {
Double_t ylo = ((RooRealVar*)vars.at(1))->getMin() ;
Double_t yhi = ((RooRealVar*)vars.at(1))->getMax() ;
if (!dynamic_cast<RooRealVar*>(yvar)) {
cout << "RooDataHist::ctor(" << GetName() << ") ERROR: dimension " << yvar->GetName() << " must be real" << endl ;
assert(0) ;
}
if (histo->GetYaxis()->GetXbins()->GetArray()) {
RooBinning ybins(histo->GetNbinsY(),histo->GetYaxis()->GetXbins()->GetArray()) ;
Double_t yloAdj = ybins.binLow(ybins.binNumber(ylo+1e-6)) ;
Double_t yhiAdj = ybins.binHigh(ybins.binNumber(yhi-1e-6)) ;
ybins.setRange(yloAdj,yhiAdj) ;
((RooRealVar*)vars.at(1))->setBinning(ybins) ;
if (fabs(yloAdj-ylo)>1e-6||fabs(yhiAdj-yhi)) {
cout << "RooDataHist::ctor(" << GetName() << "): fit range of variable " << yvar->GetName() << " expanded to nearest bin boundaries: ["
<< ylo << "," << yhi << "] --> [" << yloAdj << "," << yhiAdj << "]" << endl ;
}
yvar->setBinning(ybins) ;
ymin = ybins.rawBinNumber(yloAdj+1e-6) ;
} else {
RooBinning ybins(histo->GetYaxis()->GetXmin(),histo->GetYaxis()->GetXmax()) ;
ybins.addUniform(histo->GetNbinsY(),histo->GetYaxis()->GetXmin(),histo->GetYaxis()->GetXmax()) ;
Double_t yloAdj = ybins.binLow(ybins.binNumber(ylo+1e-6)) ;
Double_t yhiAdj = ybins.binHigh(ybins.binNumber(yhi-1e-6)) ;
ybins.setRange(yloAdj,yhiAdj) ;
((RooRealVar*)vars.at(1))->setRange(yloAdj,yhiAdj) ;
if (fabs(yloAdj-ylo)>1e-6||fabs(yhiAdj-yhi)) {
cout << "RooDataHist::ctor(" << GetName() << "): fit range of variable " << yvar->GetName() << " expanded to nearest bin boundaries: ["
<< ylo << "," << yhi << "] --> [" << yloAdj << "," << yhiAdj << "]" << endl ;
}
RooUniformBinning ybins2(yloAdj,yhiAdj,ybins.numBins()) ;
yvar->setBinning(ybins2) ;
ymin = ybins.rawBinNumber(yloAdj+1e-6) ;
}
}
RooRealVar* zvar = (RooRealVar*) (vars.at(2) ? _vars.find(vars.at(2)->GetName()) : 0 ) ;
Int_t zmin(0) ;
if (zvar) {
Double_t zlo = ((RooRealVar*)vars.at(2))->getMin() ;
Double_t zhi = ((RooRealVar*)vars.at(2))->getMax() ;
if (!dynamic_cast<RooRealVar*>(zvar)) {
cout << "RooDataHist::ctor(" << GetName() << ") ERROR: dimension " << zvar->GetName() << " must be real" << endl ;
assert(0) ;
}
if (histo->GetZaxis()->GetXbins()->GetArray()) {
RooBinning zbins(histo->GetNbinsZ(),histo->GetZaxis()->GetXbins()->GetArray()) ;
Double_t zloAdj = zbins.binLow(zbins.binNumber(zlo+1e-6)) ;
Double_t zhiAdj = zbins.binHigh(zbins.binNumber(zhi-1e-6)) ;
zbins.setRange(zloAdj,zhiAdj) ;
((RooRealVar*)vars.at(2))->setBinning(zbins) ;
if (fabs(zloAdj-zlo)>1e-6||fabs(zhiAdj-zhi)) {
cout << "RooDataHist::ctor(" << GetName() << "): fit range of variable " << zvar->GetName() << " expanded to nearest bin boundaries: ["
<< zlo << "," << zhi << "] --> [" << zloAdj << "," << zhiAdj << "]" << endl ;
}
zvar->setBinning(zbins) ;
zmin = zbins.rawBinNumber(zloAdj+1e-6) ;
} else {
RooBinning zbins(histo->GetZaxis()->GetXmin(),histo->GetZaxis()->GetXmax()) ;
zbins.addUniform(histo->GetNbinsZ(),histo->GetZaxis()->GetXmin(),histo->GetZaxis()->GetXmax()) ;
Double_t zloAdj = zbins.binLow(zbins.binNumber(zlo+1e-6)) ;
Double_t zhiAdj = zbins.binHigh(zbins.binNumber(zhi-1e-6)) ;
zbins.setRange(zloAdj,zhiAdj) ;
((RooRealVar*)vars.at(2))->setRange(zloAdj,zhiAdj) ;
if (fabs(zloAdj-zlo)>1e-6||fabs(zhiAdj-zhi)) {
cout << "RooDataHist::ctor(" << GetName() << "): fit range of variable " << zvar->GetName() << " expanded to nearest bin boundaries: ["
<< zlo << "," << zhi << "] --> [" << zloAdj << "," << zhiAdj << "]" << endl ;
}
RooUniformBinning zbins2(zloAdj,zhiAdj,zbins.numBins()) ;
zvar->setBinning(zbins2) ;
zmin = zbins.rawBinNumber(zloAdj+1e-6) ;
}
}
initialize() ;
appendToDir(this,kTRUE) ;
RooArgSet set(*xvar) ;
if (yvar) set.add(*yvar) ;
if (zvar) set.add(*zvar) ;
Int_t ix(0),iy(0),iz(0) ;
for (ix=0 ; ix < xvar->getBins() ; ix++) {
xvar->setBin(ix) ;
if (yvar) {
for (iy=0 ; iy < yvar->getBins() ; iy++) {
yvar->setBin(iy) ;
if (zvar) {
for (iz=0 ; iz < zvar->getBins() ; iz++) {
zvar->setBin(iz) ;
add(set,histo->GetBinContent(ix+1+xmin,iy+1+ymin,iz+1+zmin)*weight,TMath::Power(histo->GetBinError(ix+1+xmin,iy+1+ymin,iz+1+zmin)*weight,2)) ;
}
} else {
add(set,histo->GetBinContent(ix+1+xmin,iy+1+ymin)*weight,TMath::Power(histo->GetBinError(ix+1+xmin,iy+1+ymin)*weight,2)) ;
}
}
} else {
add(set,histo->GetBinContent(ix+1+xmin)*weight,TMath::Power(histo->GetBinError(ix+1+xmin)*weight,2)) ;
}
}
}
void RooDataHist::initialize(Bool_t fillTree)
{
_idxMult = new Int_t[_vars.getSize()] ;
_arrSize = 1 ;
_iterator->Reset() ;
RooAbsLValue* arg ;
Int_t n(0), i ;
while((arg=dynamic_cast<RooAbsLValue*>(_iterator->Next()))) {
for (i=0 ; i<n ; i++) {
_idxMult[i] *= arg->numBins() ;
}
_idxMult[n++] = 1 ;
_arrSize *= arg->numBins() ;
}
_wgt = new Double_t[_arrSize] ;
_errLo = new Double_t[_arrSize] ;
_errHi = new Double_t[_arrSize] ;
_sumw2 = new Double_t[_arrSize] ;
_binv = new Double_t[_arrSize] ;
for (i=0 ; i<_arrSize ; i++) {
_wgt[i] = 0 ;
_errLo[i] = -1 ;
_errHi[i] = -1 ;
_sumw2[i] = 0 ;
}
RooAbsArg* real ;
_iterator->Reset() ;
while((real=(RooAbsArg*)_iterator->Next())) {
if (dynamic_cast<RooAbsReal*>(real)) _realVars.add(*real) ;
}
_realIter = _realVars.createIterator() ;
if (!fillTree) return ;
Int_t ibin ;
for (ibin=0 ; ibin<_arrSize ; ibin++) {
_iterator->Reset() ;
RooAbsLValue* arg ;
Int_t i(0), idx(0), tmp(ibin) ;
Double_t binVolume(1) ;
while((arg=dynamic_cast<RooAbsLValue*>(_iterator->Next()))) {
idx = tmp / _idxMult[i] ;
tmp -= idx*_idxMult[i++] ;
RooAbsLValue* arglv = dynamic_cast<RooAbsLValue*>(arg) ;
arglv->setBin(idx) ;
binVolume *= arglv->getBinWidth(idx) ;
}
_binv[ibin] = binVolume ;
Fill() ;
}
}
RooDataHist::RooDataHist(const RooDataHist& other, const char* newname) :
RooTreeData(other,newname), RooDirItem(), _curWeight(0), _curVolume(1), _pbinv(0)
{
Int_t i ;
Int_t nVar = _vars.getSize() ;
_idxMult = new Int_t[nVar] ;
for (i=0 ; i<nVar ; i++) {
_idxMult[i] = other._idxMult[i] ;
}
_arrSize = other._arrSize ;
_wgt = new Double_t[_arrSize] ;
_errLo = new Double_t[_arrSize] ;
_errHi = new Double_t[_arrSize] ;
_binv = new Double_t[_arrSize] ;
_sumw2 = new Double_t[_arrSize] ;
for (i=0 ; i<_arrSize ; i++) {
_wgt[i] = other._wgt[i] ;
_errLo[i] = other._errLo[i] ;
_errHi[i] = other._errHi[i] ;
_sumw2[i] = other._sumw2[i] ;
_binv[i] = other._binv[i] ;
}
RooAbsArg* arg ;
_iterator->Reset() ;
while((arg=(RooAbsArg*)_iterator->Next())) {
if (dynamic_cast<RooAbsReal*>(arg)) _realVars.add(*arg) ;
}
_realIter = _realVars.createIterator() ;
appendToDir(this,kTRUE) ;
}
RooDataHist::RooDataHist(const char* name, const char* title, RooDataHist* h, const RooArgSet& varSubset,
const RooFormulaVar* cutVar, const char* cutRange, Int_t nStart, Int_t nStop, Bool_t copyCache) :
RooTreeData(name,title,h,varSubset,cutVar,cutRange,nStart,nStop,copyCache), _curWeight(0), _curVolume(1), _pbinv(0)
{
initialize(kFALSE) ;
Int_t i ;
for (i=0 ; i<_arrSize ; i++) {
_wgt[i] = h->_wgt[i] ;
_errLo[i] = h->_errLo[i] ;
_errHi[i] = h->_errHi[i] ;
_sumw2[i] = h->_sumw2[i] ;
_binv[i] = h->_binv[i] ;
}
appendToDir(this,kTRUE) ;
}
RooAbsData* RooDataHist::cacheClone(const RooArgSet* newCacheVars, const char* newName)
{
RooDataHist* dhist = new RooDataHist(newName?newName:GetName(),GetTitle(),this,*get(),0,0,0,2000000000,kTRUE) ;
RooArgSet* selCacheVars = (RooArgSet*) newCacheVars->selectCommon(dhist->_cachedVars) ;
dhist->initCache(*selCacheVars) ;
delete selCacheVars ;
return dhist ;
}
RooAbsData* RooDataHist::reduceEng(const RooArgSet& varSubset, const RooFormulaVar* cutVar, const char* cutRange,
Int_t nStart, Int_t nStop, Bool_t )
{
checkInit() ;
RooArgSet* myVarSubset = (RooArgSet*) _vars.selectCommon(varSubset) ;
RooDataHist *rdh = new RooDataHist(GetName(), GetTitle(), *myVarSubset) ;
RooFormulaVar* cloneVar = 0;
RooArgSet* tmp(0) ;
if (cutVar) {
tmp = (RooArgSet*) RooArgSet(*cutVar).snapshot() ;
if (!tmp) {
cout << "RooDataHist::reduceEng(" << GetName() << ") Couldn't deep-clone cut variable, abort," << endl ;
return 0 ;
}
cloneVar = (RooFormulaVar*) tmp->find(cutVar->GetName()) ;
cloneVar->attachDataSet(*this) ;
}
Int_t i ;
Double_t lo,hi ;
Int_t nevt = nStop < numEntries() ? nStop : numEntries() ;
TIterator* vIter = get()->createIterator() ;
for (i=nStart ; i<nevt ; i++) {
const RooArgSet* row = get(i) ;
Bool_t doSelect(kTRUE) ;
if (cutRange) {
RooAbsArg* arg ;
vIter->Reset() ;
while((arg=(RooAbsArg*)vIter->Next())) {
if (!arg->inRange(cutRange)) {
doSelect = kFALSE ;
break ;
}
}
}
if (!doSelect) continue ;
if (!cloneVar || cloneVar->getVal()) {
weightError(lo,hi,SumW2) ;
rdh->add(*row,weight(),lo*lo) ;
}
}
delete vIter ;
if (cloneVar) {
delete tmp ;
}
return rdh ;
}
RooDataHist::~RooDataHist()
{
if (_wgt) delete[] _wgt ;
if (_errLo) delete[] _errLo ;
if (_errHi) delete[] _errHi ;
if (_sumw2) delete[] _sumw2 ;
if (_binv) delete[] _binv ;
if (_pbinv) delete[] _pbinv ;
if (_idxMult) delete[] _idxMult ;
if (_realIter) delete _realIter ;
removeFromDir(this) ;
}
Int_t RooDataHist::calcTreeIndex() const
{
_iterator->Reset() ;
RooAbsLValue* arg ;
Int_t masterIdx(0), i(0) ;
while((arg=dynamic_cast<RooAbsLValue*>(_iterator->Next()))) {
masterIdx += _idxMult[i++]*arg->getBin() ;
}
return masterIdx ;
}
void RooDataHist::dump2()
{
Int_t i ;
for (i=0 ; i<_arrSize ; i++) {
if (_wgt[i]!=0) {
cout << "wgt[" << i << "] = " << _wgt[i] << "err[" << i << "] = " << _errLo[i] << " vol[" << i << "] = " << _binv[i] << endl ;
}
}
}
RooPlot *RooDataHist::plotOn(RooPlot *frame, PlotOpt o) const
{
if (o.bins) return RooTreeData::plotOn(frame,o) ;
if(0 == frame) {
cout << ClassName() << "::" << GetName() << ":plotOn: frame is null" << endl;
return 0;
}
RooAbsRealLValue *var= (RooAbsRealLValue*) frame->getPlotVar();
if(0 == var) {
cout << ClassName() << "::" << GetName()
<< ":plotOn: frame does not specify a plot variable" << endl;
return 0;
}
RooRealVar* dataVar = (RooRealVar*) _vars.find(var->GetName()) ;
if (!dataVar) {
cout << ClassName() << "::" << GetName()
<< ":plotOn: dataset doesn't contain plot frame variable" << endl;
return 0;
}
o.bins = &dataVar->getBinning() ;
return RooTreeData::plotOn(frame,o) ;
}
Double_t RooDataHist::weight(const RooArgSet& bin, Int_t intOrder, Bool_t correctForBinSize)
{
if (intOrder<0) {
cout << "RooDataHist::weight(" << GetName() << ") ERROR: interpolation order must be positive" << endl ;
return 0 ;
}
if (intOrder==0) {
_vars = bin ;
Int_t idx = calcTreeIndex() ;
if (correctForBinSize) {
return _wgt[idx] / _binv[idx] ;
} else {
return _wgt[idx] ;
}
}
_vars = bin ;
Double_t wInt(0) ;
if (_realVars.getSize()==1) {
_realIter->Reset() ;
RooRealVar* real=(RooRealVar*)_realIter->Next() ;
wInt = interpolateDim(*real,((RooAbsReal*)bin.find(real->GetName()))->getVal(), intOrder, correctForBinSize) ;
} else if (_realVars.getSize()==2) {
_realIter->Reset() ;
RooRealVar* realX=(RooRealVar*)_realIter->Next() ;
RooRealVar* realY=(RooRealVar*)_realIter->Next() ;
Double_t xval = ((RooAbsReal*)bin.find(realX->GetName()))->getVal() ;
Double_t yval = ((RooAbsReal*)bin.find(realY->GetName()))->getVal() ;
Int_t ybinC = realY->getBin() ;
Int_t ybinLo = ybinC-intOrder/2 - ((yval<realY->getBinning().binCenter(ybinC))?1:0) ;
Int_t ybinM = realY->numBins() ;
Int_t i ;
Double_t yarr[10] ;
Double_t xarr[10] ;
for (i=ybinLo ; i<=intOrder+ybinLo ; i++) {
Int_t ibin ;
if (i>=0 && i<ybinM) {
ibin = i ;
realY->setBin(ibin) ;
xarr[i-ybinLo] = realY->getVal() ;
} else if (i>=ybinM) {
ibin = 2*ybinM-i-1 ;
realY->setBin(ibin) ;
xarr[i-ybinLo] = 2*realY->getMax()-realY->getVal() ;
} else {
ibin = -i ;
realY->setBin(ibin) ;
xarr[i-ybinLo] = 2*realY->getMin()-realY->getVal() ;
}
yarr[i-ybinLo] = interpolateDim(*realX,xval,intOrder,correctForBinSize) ;
}
wInt = RooMath::interpolate(xarr,yarr,intOrder+1,yval) ;
} else {
cout << "RooDataHist::weight(" << GetName() << ") interpolation in "
<< _realVars.getSize() << " dimensions not yet implemented" << endl ;
return weight(bin,0) ;
}
if (wInt<0) wInt=0. ;
return wInt ;
}
void RooDataHist::weightError(Double_t& lo, Double_t& hi, ErrorType etype) const
{
switch (etype) {
case Poisson:
if (_curWgtErrLo>=0) {
lo = _curWgtErrLo ;
hi = _curWgtErrHi ;
return ;
}
Double_t ym,yp ;
RooHistError::instance().getPoissonInterval(Int_t(weight()+0.5),ym,yp,1) ;
_curWgtErrLo = weight()-ym ;
_curWgtErrHi = yp-weight() ;
_errLo[_curIndex] = _curWgtErrLo ;
_errHi[_curIndex] = _curWgtErrHi ;
lo = _curWgtErrLo ;
hi = _curWgtErrHi ;
return ;
case SumW2:
lo = sqrt(_curSumW2) ;
hi = sqrt(_curSumW2) ;
return ;
}
}
Double_t RooDataHist::interpolateDim(RooRealVar& dim, Double_t xval, Int_t intOrder, Bool_t correctForBinSize)
{
Int_t fbinC = dim.getBin() ;
Int_t fbinLo = fbinC-intOrder/2 - ((xval<dim.getBinning().binCenter(fbinC))?1:0) ;
Int_t fbinM = dim.numBins() ;
Int_t i ;
Double_t yarr[10] ;
Double_t xarr[10] ;
for (i=fbinLo ; i<=intOrder+fbinLo ; i++) {
Int_t ibin ;
if (i>=0 && i<fbinM) {
ibin = i ;
dim.setBin(ibin) ;
xarr[i-fbinLo] = dim.getVal() ;
} else if (i>=fbinM) {
ibin = 2*fbinM-i-1 ;
dim.setBin(ibin) ;
xarr[i-fbinLo] = 2*dim.getMax()-dim.getVal() ;
} else {
ibin = -i ;
dim.setBin(ibin) ;
xarr[i-fbinLo] = 2*dim.getMin()-dim.getVal() ;
}
Int_t idx = calcTreeIndex() ;
yarr[i-fbinLo] = _wgt[idx] ;
if (correctForBinSize) yarr[i-fbinLo] /= _binv[idx] ;
}
dim.setBin(fbinC) ;
Double_t ret = RooMath::interpolate(xarr,yarr,intOrder+1,xval) ;
return ret ;
}
void RooDataHist::add(const RooArgSet& row, Double_t weight, Double_t sumw2)
{
checkInit() ;
_vars = row ;
Int_t idx = calcTreeIndex() ;
_wgt[idx] += weight ;
_sumw2[idx] += (sumw2>0?sumw2:weight*weight) ;
_errLo[idx] = -1 ;
_errHi[idx] = -1 ;
}
void RooDataHist::set(const RooArgSet& row, Double_t weight, Double_t wgtErrLo, Double_t wgtErrHi)
{
checkInit() ;
_vars = row ;
Int_t idx = calcTreeIndex() ;
_wgt[idx] = weight ;
_errLo[idx] = wgtErrLo ;
_errHi[idx] = wgtErrHi ;
}
void RooDataHist::set(const RooArgSet& row, Double_t weight, Double_t wgtErr)
{
checkInit() ;
_vars = row ;
Int_t idx = calcTreeIndex() ;
_wgt[idx] = weight ;
_errLo[idx] = wgtErr ;
_errHi[idx] = wgtErr ;
}
void RooDataHist::add(const RooAbsData& dset, const char* cut, Double_t weight)
{
RooFormulaVar cutVar("select",cut,*dset.get()) ;
add(dset,&cutVar,weight) ;
}
void RooDataHist::add(const RooAbsData& dset, const RooFormulaVar* cutVar, Double_t weight)
{
checkInit() ;
RooFormulaVar* cloneVar = 0;
RooArgSet* tmp(0) ;
if (cutVar) {
tmp = (RooArgSet*) RooArgSet(*cutVar).snapshot() ;
if (!tmp) {
cout << "RooDataHist::add(" << GetName() << ") Couldn't deep-clone cut variable, abort," << endl ;
return ;
}
cloneVar = (RooFormulaVar*) tmp->find(cutVar->GetName()) ;
cloneVar->attachDataSet(dset) ;
}
Int_t i ;
for (i=0 ; i<dset.numEntries() ; i++) {
const RooArgSet* row = dset.get(i) ;
if (!cloneVar || cloneVar->getVal()) {
add(*row,weight*dset.weight()) ;
}
}
if (cloneVar) {
delete tmp ;
}
}
Double_t RooDataHist::sum(Bool_t correctForBinSize) const
{
Int_t i ;
Double_t total(0) ;
for (i=0 ; i<_arrSize ; i++) {
Double_t binVolume = correctForBinSize ? _binv[i] : 1.0 ;
total += _wgt[i]/binVolume ;
}
return total ;
}
Double_t RooDataHist::sum(const RooArgSet& sumSet, const RooArgSet& sliceSet, Bool_t correctForBinSize)
{
RooArgSet sliceOnlySet(sliceSet) ;
sliceOnlySet.remove(sumSet,kTRUE,kTRUE) ;
_vars = sliceOnlySet ;
calculatePartialBinVolume(sliceOnlySet) ;
TIterator* ssIter = sumSet.createIterator() ;
RooAbsArg* arg ;
Bool_t* mask = new Bool_t[_vars.getSize()] ;
Int_t* refBin = new Int_t[_vars.getSize()] ;
Int_t i(0) ;
_iterator->Reset() ;
while((arg=(RooAbsArg*)_iterator->Next())) {
if (sumSet.find(arg->GetName())) {
mask[i] = kFALSE ;
} else {
mask[i] = kTRUE ;
refBin[i] = (dynamic_cast<RooAbsLValue*>(arg))->getBin() ;
}
i++ ;
}
Double_t total(0) ;
Int_t ibin ;
for (ibin=0 ; ibin<_arrSize ; ibin++) {
Int_t idx(0), tmp(ibin), ivar(0) ;
Bool_t skip(kFALSE) ;
_iterator->Reset() ;
while((!skip && (arg=(RooAbsArg*)_iterator->Next()))) {
idx = tmp / _idxMult[ivar] ;
tmp -= idx*_idxMult[ivar] ;
if (mask[ivar] && idx!=refBin[ivar]) skip=kTRUE ;
ivar++ ;
}
if (!skip) {
Double_t binVolume = correctForBinSize ? _pbinv[ibin] : 1.0 ;
total += _wgt[ibin]/binVolume ;
}
}
delete ssIter ;
delete[] mask ;
delete[] refBin ;
return total ;
}
void RooDataHist::calculatePartialBinVolume(const RooArgSet& dimSet) const
{
if (_pbinv==0) {
_pbinv = new Double_t[_arrSize] ;
} else {
if (RooNameSet(dimSet)==_pbinvCache) {
return ;
}
}
Bool_t* selDim = new Bool_t[_vars.getSize()] ;
_iterator->Reset() ;
RooAbsArg* v ;
Int_t i(0) ;
while((v=(RooAbsArg*)_iterator->Next())) {
selDim[i++] = dimSet.find(v->GetName()) ? kTRUE : kFALSE ;
}
Int_t ibin ;
for (ibin=0 ; ibin<_arrSize ; ibin++) {
_iterator->Reset() ;
RooAbsLValue* arg ;
Int_t i(0), idx(0), tmp(ibin) ;
Double_t binVolume(1) ;
while((arg=dynamic_cast<RooAbsLValue*>(_iterator->Next()))) {
idx = tmp / _idxMult[i] ;
tmp -= idx*_idxMult[i++] ;
if (selDim[i-1]) {
RooAbsLValue* arglv = dynamic_cast<RooAbsLValue*>(arg) ;
binVolume *= arglv->getBinWidth(idx) ;
}
}
_pbinv[ibin] = binVolume ;
}
delete[] selDim ;
_pbinvCache.refill(dimSet) ;
}
Int_t RooDataHist::numEntries(Bool_t useWeights) const
{
if (!useWeights) return RooTreeData::numEntries() ;
return Int_t(sumEntries()) ;
}
Double_t RooDataHist::sumEntries(const char* cutSpec, const char* cutRange) const
{
if (cutSpec==0 && cutRange==0) {
Int_t i ;
Double_t n(0) ;
for (i=0 ; i<_arrSize ; i++) {
n+= _wgt[i] ;
}
return n ;
} else {
RooFormula* select = 0 ;
if (cutSpec) {
select = new RooFormula("select",cutSpec,*get()) ;
}
Double_t sumw(0) ;
Int_t i ;
for (i=0 ; i<GetEntries() ; i++) {
get(i) ;
if (select && select->eval()==0.) continue ;
if (cutRange && !_vars.allInRange(cutRange)) continue ;
sumw += weight() ;
}
if (select) delete select ;
return sumw ;
}
}
void RooDataHist::reset()
{
RooTreeData::reset() ;
Int_t i ;
for (i=0 ; i<_arrSize ; i++) {
_wgt[i] = 0. ;
_errLo[i] = -1 ;
_errHi[i] = -1 ;
}
_curWeight = 0 ;
_curWgtErrLo = -1 ;
_curWgtErrHi = -1 ;
_curVolume = 1 ;
}
const RooArgSet* RooDataHist::get(Int_t masterIdx) const
{
_curWeight = _wgt[masterIdx] ;
_curWgtErrLo = _errLo[masterIdx] ;
_curWgtErrHi = _errHi[masterIdx] ;
_curSumW2 = _sumw2[masterIdx] ;
_curVolume = _binv[masterIdx] ;
_curIndex = masterIdx ;
return RooTreeData::get(masterIdx) ;
}
const RooArgSet* RooDataHist::get(const RooArgSet& coord) const
{
((RooDataHist*)this)->_vars = coord ;
return get(calcTreeIndex()) ;
}
Double_t RooDataHist::binVolume(const RooArgSet& coord)
{
((RooDataHist*)this)->_vars = coord ;
return _binv[calcTreeIndex()] ;
}
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.