// RooDataSet is a container class to hold N-dimensional binned data. Each bins central
// coordinates in N-dimensional space are represented by a RooArgSet of RooRealVar, RooCategory
// or RooStringVar objects, thus data can be binned in real and/or discrete dimensions
// END_HTML
#include "RooFit.h"
#include "Riostream.h"
#include "TH1.h"
#include "TH1.h"
#include "TDirectory.h"
#include "TMath.h"
#include "RooMsgService.h"
#include "RooDataHist.h"
#include "RooDataHistSliceIter.h"
#include "RooAbsLValue.h"
#include "RooArgList.h"
#include "RooRealVar.h"
#include "RooMath.h"
#include "RooBinning.h"
#include "RooPlot.h"
#include "RooHistError.h"
#include "RooCategory.h"
#include "RooCmdConfig.h"
#include "RooLinkedListIter.h"
#include "RooTreeDataStore.h"
#include "RooVectorDataStore.h"
#include "TTree.h"
#include "RooTrace.h"
#include "RooTreeData.h"
using namespace std ;
ClassImp(RooDataHist)
;
RooDataHist::RooDataHist() : _pbinvCacheMgr(0,10)
{
_arrSize = 0 ;
_wgt = 0 ;
_errLo = 0 ;
_errHi = 0 ;
_sumw2 = 0 ;
_binv = 0 ;
_pbinv = 0 ;
_curWeight = 0 ;
_curIndex = -1 ;
_realIter = _realVars.createIterator() ;
_binValid = 0 ;
_curSumW2 = 0 ;
_curVolume = 1 ;
_curWgtErrHi = 0 ;
_curWgtErrLo = 0 ;
_cache_sum_valid = 0 ;
TRACE_CREATE
}
RooDataHist::RooDataHist(const char *name, const char *title, const RooArgSet& vars, const char* binningName) :
RooAbsData(name,title,vars), _wgt(0), _binValid(0), _curWeight(0), _curVolume(1), _pbinv(0), _pbinvCacheMgr(0,10), _cache_sum_valid(0)
{
_dstore = (defaultStorageType==Tree) ? ((RooAbsDataStore*) new RooTreeDataStore(name,title,_vars)) :
((RooAbsDataStore*) new RooVectorDataStore(name,title,_vars)) ;
initialize(binningName) ;
_dstore->setExternalWeightArray(_wgt,_errLo,_errHi,_sumw2) ;
appendToDir(this,kTRUE) ;
TRACE_CREATE
}
RooDataHist::RooDataHist(const char *name, const char *title, const RooArgSet& vars, const RooAbsData& data, Double_t wgt) :
RooAbsData(name,title,vars), _wgt(0), _binValid(0), _curWeight(0), _curVolume(1), _pbinv(0), _pbinvCacheMgr(0,10), _cache_sum_valid(0)
{
_dstore = (defaultStorageType==Tree) ? ((RooAbsDataStore*) new RooTreeDataStore(name,title,_vars)) :
((RooAbsDataStore*) new RooVectorDataStore(name,title,_vars)) ;
initialize() ;
_dstore->setExternalWeightArray(_wgt,_errLo,_errHi,_sumw2) ;
add(data,(const RooFormulaVar*)0,wgt) ;
appendToDir(this,kTRUE) ;
TRACE_CREATE
}
RooDataHist::RooDataHist(const char *name, const char *title, const RooArgList& vars, RooCategory& indexCat,
map<string,TH1*> histMap, Double_t wgt) :
RooAbsData(name,title,RooArgSet(vars,&indexCat)),
_wgt(0), _binValid(0), _curWeight(0), _curVolume(1), _pbinv(0), _pbinvCacheMgr(0,10), _cache_sum_valid(0)
{
_dstore = (defaultStorageType==Tree) ? ((RooAbsDataStore*) new RooTreeDataStore(name,title,_vars)) :
((RooAbsDataStore*) new RooVectorDataStore(name,title,_vars)) ;
importTH1Set(vars, indexCat, histMap, wgt, kFALSE) ;
_dstore->setExternalWeightArray(_wgt,_errLo,_errHi,_sumw2) ;
TRACE_CREATE
}
RooDataHist::RooDataHist(const char *name, const char *title, const RooArgList& vars, RooCategory& indexCat,
map<string,RooDataHist*> dhistMap, Double_t wgt) :
RooAbsData(name,title,RooArgSet(vars,&indexCat)),
_wgt(0), _binValid(0), _curWeight(0), _curVolume(1), _pbinv(0), _pbinvCacheMgr(0,10), _cache_sum_valid(0)
{
_dstore = (defaultStorageType==Tree) ? ((RooAbsDataStore*) new RooTreeDataStore(name,title,_vars)) :
((RooAbsDataStore*) new RooVectorDataStore(name,title,_vars)) ;
importDHistSet(vars, indexCat, dhistMap, wgt) ;
_dstore->setExternalWeightArray(_wgt,_errLo,_errHi,_sumw2) ;
TRACE_CREATE
}
RooDataHist::RooDataHist(const char *name, const char *title, const RooArgList& vars, const TH1* hist, Double_t wgt) :
RooAbsData(name,title,vars), _wgt(0), _binValid(0), _curWeight(0), _curVolume(1), _pbinv(0), _pbinvCacheMgr(0,10), _cache_sum_valid(0)
{
_dstore = (defaultStorageType==Tree) ? ((RooAbsDataStore*) new RooTreeDataStore(name,title,_vars)) :
((RooAbsDataStore*) new RooVectorDataStore(name,title,_vars)) ;
if (vars.getSize() != hist->GetDimension()) {
coutE(InputArguments) << "RooDataHist::ctor(" << GetName() << ") ERROR: dimension of input histogram must match "
<< "number of dimension variables" << endl ;
assert(0) ;
}
importTH1(vars,*const_cast<TH1*>(hist),wgt, kFALSE) ;
_dstore->setExternalWeightArray(_wgt,_errLo,_errHi,_sumw2) ;
TRACE_CREATE
}
RooDataHist::RooDataHist(const char *name, const char *title, const RooArgList& vars, const RooCmdArg& arg1, const RooCmdArg& arg2, const RooCmdArg& arg3,
const RooCmdArg& arg4,const RooCmdArg& arg5,const RooCmdArg& arg6,const RooCmdArg& arg7,const RooCmdArg& arg8) :
RooAbsData(name,title,RooArgSet(vars,(RooAbsArg*)RooCmdConfig::decodeObjOnTheFly("RooDataHist::RooDataHist", "IndexCat",0,0,arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8))),
_wgt(0), _binValid(0), _curWeight(0), _curVolume(1), _pbinv(0), _pbinvCacheMgr(0,10), _cache_sum_valid(0)
{
_dstore = (defaultStorageType==Tree) ? ((RooAbsDataStore*) new RooTreeDataStore(name,title,_vars)) :
((RooAbsDataStore*) new RooVectorDataStore(name,title,_vars)) ;
RooCmdConfig pc(Form("RooDataHist::ctor(%s)",GetName())) ;
pc.defineObject("impHist","ImportHisto",0) ;
pc.defineInt("impDens","ImportHisto",0) ;
pc.defineObject("indexCat","IndexCat",0) ;
pc.defineObject("impSliceHist","ImportHistoSlice",0,0,kTRUE) ;
pc.defineString("impSliceState","ImportHistoSlice",0,"",kTRUE) ;
pc.defineObject("impSliceDHist","ImportDataHistSlice",0,0,kTRUE) ;
pc.defineString("impSliceDState","ImportDataHistSlice",0,"",kTRUE) ;
pc.defineDouble("weight","Weight",0,1) ;
pc.defineObject("dummy1","ImportDataHistSliceMany",0) ;
pc.defineObject("dummy2","ImportHistoSliceMany",0) ;
pc.defineMutex("ImportHisto","ImportHistoSlice","ImportDataHistSlice") ;
pc.defineDependency("ImportHistoSlice","IndexCat") ;
pc.defineDependency("ImportDataHistSlice","IndexCat") ;
RooLinkedList l ;
l.Add((TObject*)&arg1) ; l.Add((TObject*)&arg2) ;
l.Add((TObject*)&arg3) ; l.Add((TObject*)&arg4) ;
l.Add((TObject*)&arg5) ; l.Add((TObject*)&arg6) ;
l.Add((TObject*)&arg7) ; l.Add((TObject*)&arg8) ;
pc.process(l) ;
if (!pc.ok(kTRUE)) {
assert(0) ;
return ;
}
TH1* impHist = static_cast<TH1*>(pc.getObject("impHist")) ;
Bool_t impDens = pc.getInt("impDens") ;
Double_t initWgt = pc.getDouble("weight") ;
const char* impSliceNames = pc.getString("impSliceState","",kTRUE) ;
const RooLinkedList& impSliceHistos = pc.getObjectList("impSliceHist") ;
RooCategory* indexCat = static_cast<RooCategory*>(pc.getObject("indexCat")) ;
const char* impSliceDNames = pc.getString("impSliceDState","",kTRUE) ;
const RooLinkedList& impSliceDHistos = pc.getObjectList("impSliceDHist") ;
if (impHist) {
importTH1(vars,*impHist,initWgt, impDens) ;
} else if (indexCat) {
if (impSliceHistos.GetSize()>0) {
map<string,TH1*> hmap ;
char tmp[1024] ;
strlcpy(tmp,impSliceNames,1024) ;
char* token = strtok(tmp,",") ;
TIterator* hiter = impSliceHistos.MakeIterator() ;
while(token) {
hmap[token] = (TH1*) hiter->Next() ;
token = strtok(0,",") ;
}
importTH1Set(vars,*indexCat,hmap,initWgt,kFALSE) ;
} else {
map<string,RooDataHist*> dmap ;
char tmp[1024] ;
strlcpy(tmp,impSliceDNames,1024) ;
char* token = strtok(tmp,",") ;
TIterator* hiter = impSliceDHistos.MakeIterator() ;
while(token) {
dmap[token] = (RooDataHist*) hiter->Next() ;
token = strtok(0,",") ;
}
importDHistSet(vars,*indexCat,dmap,initWgt) ;
}
} else {
initialize() ;
appendToDir(this,kTRUE) ;
}
_dstore->setExternalWeightArray(_wgt,_errLo,_errHi,_sumw2) ;
TRACE_CREATE
}
void RooDataHist::importTH1(const RooArgList& vars, TH1& histo, Double_t wgt, Bool_t doDensityCorrection)
{
Int_t offset[3] ;
adjustBinning(vars,histo,offset) ;
initialize() ;
appendToDir(this,kTRUE) ;
RooRealVar* xvar = (RooRealVar*) _vars.find(vars.at(0)->GetName()) ;
RooRealVar* yvar = (RooRealVar*) (vars.at(1) ? _vars.find(vars.at(1)->GetName()) : 0 ) ;
RooRealVar* zvar = (RooRealVar*) (vars.at(2) ? _vars.find(vars.at(2)->GetName()) : 0 ) ;
Int_t xmin(0),ymin(0),zmin(0) ;
RooArgSet vset(*xvar) ;
Double_t volume = xvar->getMax()-xvar->getMin() ;
xmin = offset[0] ;
if (yvar) {
vset.add(*yvar) ;
ymin = offset[1] ;
volume *= (yvar->getMax()-yvar->getMin()) ;
}
if (zvar) {
vset.add(*zvar) ;
zmin = offset[2] ;
volume *= (zvar->getMax()-zvar->getMin()) ;
}
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) ;
Double_t bv = doDensityCorrection ? binVolume(vset) : 1;
add(vset,bv*histo.GetBinContent(ix+1+xmin,iy+1+ymin,iz+1+zmin)*wgt,bv*TMath::Power(histo.GetBinError(ix+1+xmin,iy+1+ymin,iz+1+zmin)*wgt,2)) ;
}
} else {
Double_t bv = doDensityCorrection ? binVolume(vset) : 1;
add(vset,bv*histo.GetBinContent(ix+1+xmin,iy+1+ymin)*wgt,bv*TMath::Power(histo.GetBinError(ix+1+xmin,iy+1+ymin)*wgt,2)) ;
}
}
} else {
Double_t bv = doDensityCorrection ? binVolume(vset) : 1 ;
add(vset,bv*histo.GetBinContent(ix+1+xmin)*wgt,bv*TMath::Power(histo.GetBinError(ix+1+xmin)*wgt,2)) ;
}
}
}
void RooDataHist::importTH1Set(const RooArgList& vars, RooCategory& indexCat, map<string,TH1*> hmap, Double_t wgt, Bool_t doDensityCorrection)
{
RooCategory* icat = (RooCategory*) _vars.find(indexCat.GetName()) ;
TH1* histo(0) ;
Bool_t init(kFALSE) ;
for (map<string,TH1*>::iterator hiter = hmap.begin() ; hiter!=hmap.end() ; ++hiter) {
if (!histo) {
histo = hiter->second ;
}
if (!indexCat.lookupType(hiter->first.c_str())) {
indexCat.defineType(hiter->first.c_str()) ;
coutI(InputArguments) << "RooDataHist::importTH1Set(" << GetName() << ") defining state \"" << hiter->first << "\" in index category " << indexCat.GetName() << endl ;
}
if (!icat->lookupType(hiter->first.c_str())) {
icat->defineType(hiter->first.c_str()) ;
}
}
if (histo && (vars.getSize() != histo->GetDimension())) {
coutE(InputArguments) << "RooDataHist::ctor(" << GetName() << ") ERROR: dimension of input histogram must match "
<< "number of continuous variables" << endl ;
assert(0) ;
}
Int_t offset[3] ;
adjustBinning(vars,*histo,offset) ;
if (!init) {
initialize() ;
appendToDir(this,kTRUE) ;
init = kTRUE ;
}
RooRealVar* xvar = (RooRealVar*) _vars.find(vars.at(0)->GetName()) ;
RooRealVar* yvar = (RooRealVar*) (vars.at(1) ? _vars.find(vars.at(1)->GetName()) : 0 ) ;
RooRealVar* zvar = (RooRealVar*) (vars.at(2) ? _vars.find(vars.at(2)->GetName()) : 0 ) ;
Int_t xmin(0),ymin(0),zmin(0) ;
RooArgSet vset(*xvar) ;
Double_t volume = xvar->getMax()-xvar->getMin() ;
xmin = offset[0] ;
if (yvar) {
vset.add(*yvar) ;
ymin = offset[1] ;
volume *= (yvar->getMax()-yvar->getMin()) ;
}
if (zvar) {
vset.add(*zvar) ;
zmin = offset[2] ;
volume *= (zvar->getMax()-zvar->getMin()) ;
}
Double_t avgBV = volume / numEntries() ;
Int_t ic(0),ix(0),iy(0),iz(0) ;
for (ic=0 ; ic < icat->numBins(0) ; ic++) {
icat->setBin(ic) ;
histo = hmap[icat->getLabel()] ;
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) ;
Double_t bv = doDensityCorrection ? binVolume(vset)/avgBV : 1;
add(vset,bv*histo->GetBinContent(ix+1+xmin,iy+1+ymin,iz+1+zmin)*wgt,bv*TMath::Power(histo->GetBinError(ix+1+xmin,iy+1+ymin,iz+1+zmin)*wgt,2)) ;
}
} else {
Double_t bv = doDensityCorrection ? binVolume(vset)/avgBV : 1;
add(vset,bv*histo->GetBinContent(ix+1+xmin,iy+1+ymin)*wgt,bv*TMath::Power(histo->GetBinError(ix+1+xmin,iy+1+ymin)*wgt,2)) ;
}
}
} else {
Double_t bv = doDensityCorrection ? binVolume(vset)/avgBV : 1;
add(vset,bv*histo->GetBinContent(ix+1+xmin)*wgt,bv*TMath::Power(histo->GetBinError(ix+1+xmin)*wgt,2)) ;
}
}
}
}
void RooDataHist::importDHistSet(const RooArgList& , RooCategory& indexCat, std::map<std::string,RooDataHist*> dmap, Double_t initWgt)
{
RooCategory* icat = (RooCategory*) _vars.find(indexCat.GetName()) ;
for (map<string,RooDataHist*>::iterator diter = dmap.begin() ; diter!=dmap.end() ; ++diter) {
if (!indexCat.lookupType(diter->first.c_str())) {
indexCat.defineType(diter->first.c_str()) ;
coutI(InputArguments) << "RooDataHist::importDHistSet(" << GetName() << ") defining state \"" << diter->first << "\" in index category " << indexCat.GetName() << endl ;
}
if (!icat->lookupType(diter->first.c_str())) {
icat->defineType(diter->first.c_str()) ;
}
}
initialize() ;
appendToDir(this,kTRUE) ;
for (map<string,RooDataHist*>::iterator diter = dmap.begin() ; diter!=dmap.end() ; ++diter) {
RooDataHist* dhist = diter->second ;
icat->setLabel(diter->first.c_str()) ;
for (Int_t i=0 ; i<dhist->numEntries() ; i++) {
_vars = *dhist->get(i) ;
add(_vars,dhist->weight()*initWgt, pow(dhist->weightError(SumW2),2) ) ;
}
}
}
void RooDataHist::adjustBinning(const RooArgList& vars, TH1& href, Int_t* offset)
{
RooRealVar* xvar = (RooRealVar*) _vars.find(*vars.at(0)) ;
if (!dynamic_cast<RooRealVar*>(xvar)) {
coutE(InputArguments) << "RooDataHist::adjustBinning(" << 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 (href.GetXaxis()->GetXbins()->GetArray()) {
RooBinning xbins(href.GetNbinsX(),href.GetXaxis()->GetXbins()->GetArray()) ;
Double_t tolerance = 1e-6*xbins.averageBinWidth() ;
Double_t xloAdj = xbins.binLow(xbins.binNumber(xlo+tolerance)) ;
Double_t xhiAdj = xbins.binHigh(xbins.binNumber(xhi-tolerance)) ;
xbins.setRange(xloAdj,xhiAdj) ;
((RooRealVar*)vars.at(0))->setBinning(xbins) ;
if (fabs(xloAdj-xlo)>tolerance||fabs(xhiAdj-xhi)<tolerance) {
coutI(DataHandling) << "RooDataHist::adjustBinning(" << GetName() << "): fit range of variable " << xvar->GetName() << " expanded to nearest bin boundaries: ["
<< xlo << "," << xhi << "] --> [" << xloAdj << "," << xhiAdj << "]" << endl ;
}
xvar->setBinning(xbins) ;
xmin = xbins.rawBinNumber(xloAdj+tolerance) ;
if (offset) {
offset[0] = xmin ;
}
} else {
RooBinning xbins(href.GetXaxis()->GetXmin(),href.GetXaxis()->GetXmax()) ;
xbins.addUniform(href.GetNbinsX(),href.GetXaxis()->GetXmin(),href.GetXaxis()->GetXmax()) ;
Double_t tolerance = 1e-6*xbins.averageBinWidth() ;
Double_t xloAdj = xbins.binLow(xbins.binNumber(xlo+tolerance)) ;
Double_t xhiAdj = xbins.binHigh(xbins.binNumber(xhi-tolerance)) ;
xbins.setRange(xloAdj,xhiAdj) ;
((RooRealVar*)vars.at(0))->setRange(xloAdj,xhiAdj) ;
if (fabs(xloAdj-xlo)>tolerance||fabs(xhiAdj-xhi)<tolerance) {
coutI(DataHandling) << "RooDataHist::adjustBinning(" << 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+tolerance) ;
if (offset) {
offset[0] = xmin ;
}
}
RooRealVar* yvar = (RooRealVar*) (vars.at(1) ? _vars.find(*vars.at(1)) : 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)) {
coutE(InputArguments) << "RooDataHist::adjustBinning(" << GetName() << ") ERROR: dimension " << yvar->GetName() << " must be real" << endl ;
assert(0) ;
}
if (href.GetYaxis()->GetXbins()->GetArray()) {
RooBinning ybins(href.GetNbinsY(),href.GetYaxis()->GetXbins()->GetArray()) ;
Double_t tolerance = 1e-6*ybins.averageBinWidth() ;
Double_t yloAdj = ybins.binLow(ybins.binNumber(ylo+tolerance)) ;
Double_t yhiAdj = ybins.binHigh(ybins.binNumber(yhi-tolerance)) ;
ybins.setRange(yloAdj,yhiAdj) ;
((RooRealVar*)vars.at(1))->setBinning(ybins) ;
if (fabs(yloAdj-ylo)>tolerance||fabs(yhiAdj-yhi)<tolerance) {
coutI(DataHandling) << "RooDataHist::adjustBinning(" << GetName() << "): fit range of variable " << yvar->GetName() << " expanded to nearest bin boundaries: ["
<< ylo << "," << yhi << "] --> [" << yloAdj << "," << yhiAdj << "]" << endl ;
}
yvar->setBinning(ybins) ;
ymin = ybins.rawBinNumber(yloAdj+tolerance) ;
if (offset) {
offset[1] = ymin ;
}
} else {
RooBinning ybins(href.GetYaxis()->GetXmin(),href.GetYaxis()->GetXmax()) ;
ybins.addUniform(href.GetNbinsY(),href.GetYaxis()->GetXmin(),href.GetYaxis()->GetXmax()) ;
Double_t tolerance = 1e-6*ybins.averageBinWidth() ;
Double_t yloAdj = ybins.binLow(ybins.binNumber(ylo+tolerance)) ;
Double_t yhiAdj = ybins.binHigh(ybins.binNumber(yhi-tolerance)) ;
ybins.setRange(yloAdj,yhiAdj) ;
((RooRealVar*)vars.at(1))->setRange(yloAdj,yhiAdj) ;
if (fabs(yloAdj-ylo)>tolerance||fabs(yhiAdj-yhi)<tolerance) {
coutI(DataHandling) << "RooDataHist::adjustBinning(" << 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+tolerance) ;
if (offset) {
offset[1] = ymin ;
}
}
}
RooRealVar* zvar = (RooRealVar*) (vars.at(2) ? _vars.find(*vars.at(2)) : 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)) {
coutE(InputArguments) << "RooDataHist::adjustBinning(" << GetName() << ") ERROR: dimension " << zvar->GetName() << " must be real" << endl ;
assert(0) ;
}
if (href.GetZaxis()->GetXbins()->GetArray()) {
RooBinning zbins(href.GetNbinsZ(),href.GetZaxis()->GetXbins()->GetArray()) ;
Double_t tolerance = 1e-6*zbins.averageBinWidth() ;
Double_t zloAdj = zbins.binLow(zbins.binNumber(zlo+tolerance)) ;
Double_t zhiAdj = zbins.binHigh(zbins.binNumber(zhi-tolerance)) ;
zbins.setRange(zloAdj,zhiAdj) ;
((RooRealVar*)vars.at(2))->setBinning(zbins) ;
if (fabs(zloAdj-zlo)>tolerance||fabs(zhiAdj-zhi)<tolerance) {
coutI(DataHandling) << "RooDataHist::adjustBinning(" << GetName() << "): fit range of variable " << zvar->GetName() << " expanded to nearest bin boundaries: ["
<< zlo << "," << zhi << "] --> [" << zloAdj << "," << zhiAdj << "]" << endl ;
}
zvar->setBinning(zbins) ;
zmin = zbins.rawBinNumber(zloAdj+tolerance) ;
if (offset) {
offset[2] = zmin ;
}
} else {
RooBinning zbins(href.GetZaxis()->GetXmin(),href.GetZaxis()->GetXmax()) ;
zbins.addUniform(href.GetNbinsZ(),href.GetZaxis()->GetXmin(),href.GetZaxis()->GetXmax()) ;
Double_t tolerance = 1e-6*zbins.averageBinWidth() ;
Double_t zloAdj = zbins.binLow(zbins.binNumber(zlo+tolerance)) ;
Double_t zhiAdj = zbins.binHigh(zbins.binNumber(zhi-tolerance)) ;
zbins.setRange(zloAdj,zhiAdj) ;
((RooRealVar*)vars.at(2))->setRange(zloAdj,zhiAdj) ;
if (fabs(zloAdj-zlo)>tolerance||fabs(zhiAdj-zhi)<tolerance) {
coutI(DataHandling) << "RooDataHist::adjustBinning(" << 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+tolerance) ;
if (offset) {
offset[2] = zmin ;
}
}
}
}
void RooDataHist::initialize(const char* binningName, Bool_t fillTree)
{
RooAbsArg* real ;
_iterator->Reset() ;
while((real=(RooAbsArg*)_iterator->Next())) {
if (dynamic_cast<RooAbsReal*>(real)) _realVars.add(*real);
}
_realIter = _realVars.createIterator() ;
_iterator->Reset();
RooAbsArg* rvarg;
while((rvarg=(RooAbsArg*)_iterator->Next())) {
if (binningName) {
RooRealVar* rrv = dynamic_cast<RooRealVar*>(rvarg);
if (rrv) {
rrv->setBinning(rrv->getBinning(binningName));
}
}
_lvvars.push_back(dynamic_cast<RooAbsLValue*>(rvarg));
const RooAbsBinning* binning = dynamic_cast<RooAbsLValue*>(rvarg)->getBinningPtr(0);
_lvbins.push_back(binning ? binning->clone() : 0);
}
_idxMult.resize(_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() ;
}
if (!_wgt) {
_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] ;
if (fillTree==kFALSE) {
_dstore->setExternalWeightArray(_wgt,_errLo,_errHi,_sumw2) ;
}
for (i=0 ; i<_arrSize ; i++) {
_wgt[i] = 0 ;
_errLo[i] = -1 ;
_errHi[i] = -1 ;
_sumw2[i] = 0 ;
}
}
if (!fillTree) return ;
Int_t ibin ;
for (ibin=0 ; ibin<_arrSize ; ibin++) {
_iterator->Reset() ;
RooAbsLValue* arg2 ;
Int_t j(0), idx(0), tmp(ibin) ;
Double_t theBinVolume(1) ;
while((arg2=dynamic_cast<RooAbsLValue*>(_iterator->Next()))) {
idx = tmp / _idxMult[j] ;
tmp -= idx*_idxMult[j++] ;
RooAbsLValue* arglv = dynamic_cast<RooAbsLValue*>(arg2) ;
arglv->setBin(idx) ;
theBinVolume *= arglv->getBinWidth(idx) ;
}
_binv[ibin] = theBinVolume ;
fill() ;
}
}
void RooDataHist::checkBinBounds() const
{
if (!_binbounds.empty()) return;
for (std::vector<const RooAbsBinning*>::const_iterator it = _lvbins.begin();
_lvbins.end() != it; ++it) {
_binbounds.push_back(std::vector<Double_t>());
if (*it) {
std::vector<Double_t>& bounds = _binbounds.back();
bounds.reserve(2 * (*it)->numBins());
for (Int_t i = 0; i < (*it)->numBins(); ++i) {
bounds.push_back((*it)->binLow(i));
bounds.push_back((*it)->binHigh(i));
}
}
}
}
RooDataHist::RooDataHist(const RooDataHist& other, const char* newname) :
RooAbsData(other,newname), RooDirItem(), _idxMult(other._idxMult), _binValid(0), _curWeight(0), _curVolume(1), _pbinv(0), _pbinvCacheMgr(other._pbinvCacheMgr,0), _cache_sum_valid(0)
{
Int_t 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() ;
_iterator->Reset() ;
RooAbsArg* rvarg ;
while((rvarg=(RooAbsArg*)_iterator->Next())) {
_lvvars.push_back(dynamic_cast<RooAbsLValue*>(rvarg)) ;
const RooAbsBinning* binning = dynamic_cast<RooAbsLValue*>(rvarg)->getBinningPtr(0) ;
_lvbins.push_back(binning ? binning->clone() : 0) ;
}
_dstore->setExternalWeightArray(_wgt,_errLo,_errHi,_sumw2) ;
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) :
RooAbsData(name,title,varSubset),
_wgt(0), _binValid(0), _curWeight(0), _curVolume(1), _pbinv(0), _pbinvCacheMgr(0,10), _cache_sum_valid(0)
{
_dstore = new RooTreeDataStore(name,title,*h->_dstore,_vars,cutVar,cutRange,nStart,nStop,copyCache) ;
initialize(0,kFALSE) ;
_dstore->setExternalWeightArray(_wgt,_errLo,_errHi,_sumw2) ;
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) ;
TRACE_CREATE
}
RooAbsData* RooDataHist::cacheClone(const RooAbsArg* newCacheOwner, const RooArgSet* newCacheVars, const char* newName)
{
checkInit() ;
RooDataHist* dhist = new RooDataHist(newName?newName:GetName(),GetTitle(),this,*get(),0,0,0,2000000000,kTRUE) ;
RooArgSet* selCacheVars = (RooArgSet*) newCacheVars->selectCommon(dhist->_cachedVars) ;
dhist->attachCache(newCacheOwner, *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) ;
delete myVarSubset ;
RooFormulaVar* cloneVar = 0;
RooArgSet* tmp(0) ;
if (cutVar) {
tmp = (RooArgSet*) RooArgSet(*cutVar).snapshot() ;
if (!tmp) {
coutE(DataHandling) << "RooDataHist::reduceEng(" << GetName() << ") Couldn't deep-clone cut variable, abort," << endl ;
return 0 ;
}
cloneVar = (RooFormulaVar*) tmp->find(*cutVar) ;
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 (_realIter) delete _realIter ;
if (_binValid) delete[] _binValid ;
vector<const RooAbsBinning*>::iterator iter = _lvbins.begin() ;
while(iter!=_lvbins.end()) {
delete *iter ;
iter++ ;
}
removeFromDir(this) ;
TRACE_DESTROY
}
Int_t RooDataHist::getIndex(const RooArgSet& coord, Bool_t fast)
{
checkInit() ;
if (fast) {
_vars.assignFast(coord,kFALSE) ;
} else {
_vars.assignValueOnly(coord) ;
}
return calcTreeIndex() ;
}
Int_t RooDataHist::calcTreeIndex() const
{
Int_t masterIdx(0), i(0) ;
vector<RooAbsLValue*>::const_iterator iter = _lvvars.begin() ;
vector<const RooAbsBinning*>::const_iterator biter = _lvbins.begin() ;
for (;iter!=_lvvars.end() ; ++iter) {
const RooAbsBinning* binning = (*biter) ;
masterIdx += _idxMult[i++]*(*iter)->getBin(binning) ;
biter++ ;
}
return masterIdx ;
}
void RooDataHist::dump2()
{
cout << "_arrSize = " << _arrSize << endl ;
for (Int_t i=0 ; i<_arrSize ; i++) {
cout << "wgt[" << i << "] = " << _wgt[i] << "sumw2[" << i << "] = " << _sumw2[i] << " vol[" << i << "] = " << _binv[i] << endl ;
}
}
RooPlot *RooDataHist::plotOn(RooPlot *frame, PlotOpt o) const
{
checkInit() ;
if (o.bins) return RooAbsData::plotOn(frame,o) ;
if(0 == frame) {
coutE(InputArguments) << ClassName() << "::" << GetName() << ":plotOn: frame is null" << endl;
return 0;
}
RooAbsRealLValue *var= (RooAbsRealLValue*) frame->getPlotVar();
if(0 == var) {
coutE(InputArguments) << ClassName() << "::" << GetName()
<< ":plotOn: frame does not specify a plot variable" << endl;
return 0;
}
RooRealVar* dataVar = (RooRealVar*) _vars.find(*var) ;
if (!dataVar) {
coutE(InputArguments) << ClassName() << "::" << GetName()
<< ":plotOn: dataset doesn't contain plot frame variable" << endl;
return 0;
}
o.bins = &dataVar->getBinning() ;
o.correctForBinWidth = kFALSE ;
return RooAbsData::plotOn(frame,o) ;
}
Double_t RooDataHist::weightSquared() const {
return _curSumW2 ;
}
Double_t RooDataHist::weight(const RooArgSet& bin, Int_t intOrder, Bool_t correctForBinSize, Bool_t cdfBoundaries, Bool_t oneSafe)
{
checkInit() ;
if (intOrder<0) {
coutE(InputArguments) << "RooDataHist::weight(" << GetName() << ") ERROR: interpolation order must be positive" << endl ;
return 0 ;
}
if (intOrder==0) {
_vars.assignValueOnly(bin,oneSafe) ;
Int_t idx = calcTreeIndex() ;
if (correctForBinSize) {
return _wgt[idx] / _binv[idx] ;
} else {
return _wgt[idx] ;
}
}
_vars.assignValueOnly(bin) ;
Double_t wInt(0) ;
if (_realVars.getSize()==1) {
RooFIter realIter = _realVars.fwdIterator() ;
RooRealVar* real=(RooRealVar*)realIter.next() ;
const RooAbsBinning* binning = real->getBinningPtr(0) ;
wInt = interpolateDim(*real,binning,((RooAbsReal*)bin.find(*real))->getVal(), intOrder, correctForBinSize, cdfBoundaries) ;
} else if (_realVars.getSize()==2) {
RooFIter realIter = _realVars.fwdIterator() ;
RooRealVar* realX=(RooRealVar*)realIter.next() ;
RooRealVar* realY=(RooRealVar*)realIter.next() ;
Double_t xval = ((RooAbsReal*)bin.find(*realX))->getVal() ;
Double_t yval = ((RooAbsReal*)bin.find(*realY))->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] ;
const RooAbsBinning* binning = realX->getBinningPtr(0) ;
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 -1;
realY->setBin(ibin) ;
xarr[i-ybinLo] = 2*realY->getMin()-realY->getVal() ;
}
yarr[i-ybinLo] = interpolateDim(*realX,binning,xval,intOrder,correctForBinSize,kFALSE) ;
}
if (gDebug>7) {
cout << "RooDataHist interpolating data is" << endl ;
cout << "xarr = " ;
for (int q=0; q<=intOrder ; q++) cout << xarr[q] << " " ;
cout << " yarr = " ;
for (int q=0; q<=intOrder ; q++) cout << yarr[q] << " " ;
cout << endl ;
}
wInt = RooMath::interpolate(xarr,yarr,intOrder+1,yval) ;
} else {
coutE(InputArguments) << "RooDataHist::weight(" << GetName() << ") interpolation in "
<< _realVars.getSize() << " dimensions not yet implemented" << endl ;
return weight(bin,0) ;
}
return wInt ;
}
void RooDataHist::weightError(Double_t& lo, Double_t& hi, ErrorType etype) const
{
checkInit() ;
switch (etype) {
case Auto:
throw string(Form("RooDataHist::weightError(%s) error type Auto not allowed here",GetName())) ;
break ;
case Expected:
throw string(Form("RooDataHist::weightError(%s) error type Expected not allowed here",GetName())) ;
break ;
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 ;
case None:
lo = 0 ;
hi = 0 ;
return ;
}
}
Double_t RooDataHist::interpolateDim(RooRealVar& dim, const RooAbsBinning* binning, Double_t xval, Int_t intOrder, Bool_t correctForBinSize, Bool_t cdfBoundaries)
{
Int_t fbinC = dim.getBin(*binning) ;
Int_t fbinLo = fbinC-intOrder/2 - ((xval<binning->binCenter(fbinC))?1:0) ;
Int_t fbinM = dim.numBins(*binning) ;
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.setBinFast(ibin,*binning) ;
xarr[i-fbinLo] = dim.getVal() ;
Int_t idx = calcTreeIndex() ;
yarr[i-fbinLo] = _wgt[idx] ;
if (correctForBinSize) yarr[i-fbinLo] /= _binv[idx] ;
} else if (i>=fbinM) {
ibin = 2*fbinM-i-1 ;
dim.setBinFast(ibin,*binning) ;
if (cdfBoundaries) {
xarr[i-fbinLo] = dim.getMax()+1e-10*(i-fbinM+1) ;
yarr[i-fbinLo] = 1.0 ;
} else {
Int_t idx = calcTreeIndex() ;
xarr[i-fbinLo] = 2*dim.getMax()-dim.getVal() ;
yarr[i-fbinLo] = _wgt[idx] ;
if (correctForBinSize) yarr[i-fbinLo] /= _binv[idx] ;
}
} else {
ibin = -i - 1 ;
dim.setBinFast(ibin,*binning) ;
if (cdfBoundaries) {
xarr[i-fbinLo] = dim.getMin()-ibin*(1e-10) ; ;
yarr[i-fbinLo] = 0.0 ;
} else {
Int_t idx = calcTreeIndex() ;
xarr[i-fbinLo] = 2*dim.getMin()-dim.getVal() ;
yarr[i-fbinLo] = _wgt[idx] ;
if (correctForBinSize) yarr[i-fbinLo] /= _binv[idx] ;
}
}
}
dim.setBinFast(fbinC,*binning) ;
Double_t ret = RooMath::interpolate(xarr,yarr,intOrder+1,xval) ;
return ret ;
}
void RooDataHist::add(const RooArgSet& row, Double_t wgt, Double_t sumw2)
{
checkInit() ;
_vars = row ;
Int_t idx = calcTreeIndex() ;
_wgt[idx] += wgt ;
_sumw2[idx] += (sumw2>0?sumw2:wgt*wgt) ;
_errLo[idx] = -1 ;
_errHi[idx] = -1 ;
_cache_sum_valid = kFALSE ;
}
void RooDataHist::set(const RooArgSet& row, Double_t wgt, Double_t wgtErrLo, Double_t wgtErrHi)
{
checkInit() ;
_vars = row ;
Int_t idx = calcTreeIndex() ;
_wgt[idx] = wgt ;
_errLo[idx] = wgtErrLo ;
_errHi[idx] = wgtErrHi ;
_cache_sum_valid = kFALSE ;
}
void RooDataHist::set(Double_t wgt, Double_t wgtErr)
{
checkInit() ;
if (_curIndex<0) {
_curIndex = calcTreeIndex() ;
}
_wgt[_curIndex] = wgt ;
_errLo[_curIndex] = wgtErr ;
_errHi[_curIndex] = wgtErr ;
_sumw2[_curIndex] = wgtErr*wgtErr ;
_cache_sum_valid = kFALSE ;
}
void RooDataHist::set(const RooArgSet& row, Double_t wgt, Double_t wgtErr)
{
checkInit() ;
_vars = row ;
Int_t idx = calcTreeIndex() ;
_wgt[idx] = wgt ;
_errLo[idx] = wgtErr ;
_errHi[idx] = wgtErr ;
_sumw2[idx] = wgtErr*wgtErr ;
_cache_sum_valid = kFALSE ;
}
void RooDataHist::add(const RooAbsData& dset, const char* cut, Double_t wgt)
{
RooFormulaVar cutVar("select",cut,*dset.get()) ;
add(dset,&cutVar,wgt) ;
}
void RooDataHist::add(const RooAbsData& dset, const RooFormulaVar* cutVar, Double_t wgt)
{
checkInit() ;
RooFormulaVar* cloneVar = 0;
RooArgSet* tmp(0) ;
if (cutVar) {
tmp = (RooArgSet*) RooArgSet(*cutVar).snapshot() ;
if (!tmp) {
coutE(DataHandling) << "RooDataHist::add(" << GetName() << ") Couldn't deep-clone cut variable, abort," << endl ;
return ;
}
cloneVar = (RooFormulaVar*) tmp->find(*cutVar) ;
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,wgt*dset.weight(), wgt*wgt*dset.weightSquared()) ;
}
}
if (cloneVar) {
delete tmp ;
}
_cache_sum_valid = kFALSE ;
}
Double_t RooDataHist::sum(Bool_t correctForBinSize, Bool_t inverseBinCor) const
{
checkInit() ;
Int_t cache_code = 1 + (correctForBinSize?1:0) + ((correctForBinSize&&inverseBinCor)?1:0) ;
if (_cache_sum_valid==cache_code) {
return _cache_sum ;
}
Int_t i ;
Double_t total(0), carry(0);
for (i=0 ; i<_arrSize ; i++) {
Double_t theBinVolume = correctForBinSize ? (inverseBinCor ? 1/_binv[i] : _binv[i]) : 1.0 ;
Double_t y = _wgt[i]*theBinVolume - carry;
Double_t t = total + y;
carry = (t - total) - y;
total = t;
}
_cache_sum_valid=cache_code ;
_cache_sum = total ;
return total ;
}
Double_t RooDataHist::sum(const RooArgSet& sumSet, const RooArgSet& sliceSet, Bool_t correctForBinSize, Bool_t inverseBinCor)
{
checkInit() ;
RooArgSet varSave ;
varSave.addClone(_vars) ;
RooArgSet* sliceOnlySet = new RooArgSet(sliceSet) ;
sliceOnlySet->remove(sumSet,kTRUE,kTRUE) ;
_vars = *sliceOnlySet ;
calculatePartialBinVolume(*sliceOnlySet) ;
delete 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)) {
mask[i] = kFALSE ;
} else {
mask[i] = kTRUE ;
refBin[i] = (dynamic_cast<RooAbsLValue*>(arg))->getBin() ;
}
i++ ;
}
Double_t total(0), carry(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 theBinVolume = correctForBinSize ? (inverseBinCor ? 1/(*_pbinv)[i] : (*_pbinv)[i] ) : 1.0 ;
Double_t y = _wgt[ibin]*theBinVolume - carry;
Double_t t = total + y;
carry = (t - total) - y;
total = t;
}
}
delete ssIter ;
delete[] mask ;
delete[] refBin ;
_vars = varSave ;
return total ;
}
Double_t RooDataHist::sum(const RooArgSet& sumSet, const RooArgSet& sliceSet,
Bool_t correctForBinSize, Bool_t inverseBinCor,
const std::map<const RooAbsArg*, std::pair<Double_t, Double_t> >& ranges)
{
checkInit();
checkBinBounds();
RooArgSet varSave;
varSave.addClone(_vars);
{
RooArgSet sliceOnlySet(sliceSet);
sliceOnlySet.remove(sumSet,kTRUE,kTRUE);
_vars = sliceOnlySet;
}
std::vector<bool> mask(_vars.getSize());
std::vector<Int_t> refBin(_vars.getSize());
std::vector<Double_t> rangeLo(_vars.getSize(), -std::numeric_limits<Double_t>::infinity());
std::vector<Double_t> rangeHi(_vars.getSize(), +std::numeric_limits<Double_t>::infinity());
_iterator->Reset();
RooAbsArg* arg;
for (Int_t i = 0; (arg=(RooAbsArg*)_iterator->Next()); ++i) {
RooAbsArg* sumsetv = sumSet.find(*arg);
RooAbsArg* slicesetv = sliceSet.find(*arg);
mask[i] = !sumsetv;
if (mask[i]) {
refBin[i] = (dynamic_cast<RooAbsLValue*>(arg))->getBin();
}
std::map<const RooAbsArg*, std::pair<Double_t, Double_t> >::const_iterator
it = ranges.find(sumsetv ? sumsetv : slicesetv);
if (ranges.end() != it) {
rangeLo[i] = it->second.first;
rangeHi[i] = it->second.second;
}
}
Double_t total(0), carry(0);
for (Int_t ibin = 0; ibin < _arrSize; ++ibin) {
_iterator->Reset();
Bool_t skip(kFALSE);
for (Int_t ivar = 0, tmp = ibin;
(!skip && (arg=(RooAbsArg*)_iterator->Next())); ++ivar) {
const Int_t idx = tmp / _idxMult[ivar];
tmp -= idx*_idxMult[ivar];
if (mask[ivar] && idx!=refBin[ivar]) skip=kTRUE;
}
if (skip) continue;
_iterator->Reset();
Double_t theBinVolume = 1.;
for (Int_t ivar = 0, tmp = ibin;
(arg=(RooAbsArg*)_iterator->Next()); ++ivar) {
const Int_t idx = tmp / _idxMult[ivar];
tmp -= idx*_idxMult[ivar];
if (_binbounds[ivar].empty()) continue;
const Double_t binLo = _binbounds[ivar][2 * idx];
const Double_t binHi = _binbounds[ivar][2 * idx + 1];
if (binHi < rangeLo[ivar] || binLo > rangeHi[ivar]) {
theBinVolume = 0.;
break;
}
theBinVolume *=
(std::min(rangeHi[ivar], binHi) - std::max(rangeLo[ivar], binLo));
}
const Double_t corrPartial = theBinVolume / _binv[ibin];
if (0. == corrPartial) continue;
const Double_t corr = correctForBinSize ? (inverseBinCor ? 1. / _binv[ibin] : _binv[ibin] ) : 1.0;
const Double_t y = _wgt[ibin] * corr * corrPartial - carry;
const Double_t t = total + y;
carry = (t - total) - y;
total = t;
}
_vars = varSave;
return total;
}
void RooDataHist::calculatePartialBinVolume(const RooArgSet& dimSet) const
{
vector<Double_t> *pbinv = _pbinvCacheMgr.getObj(&dimSet) ;
if (pbinv) {
_pbinv = pbinv ;
return ;
}
pbinv = new vector<Double_t>(_arrSize) ;
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) ? kTRUE : kFALSE ;
}
Int_t ibin ;
for (ibin=0 ; ibin<_arrSize ; ibin++) {
_iterator->Reset() ;
RooAbsLValue* arg ;
Int_t j(0), idx(0), tmp(ibin) ;
Double_t theBinVolume(1) ;
while((arg=dynamic_cast<RooAbsLValue*>(_iterator->Next()))) {
idx = tmp / _idxMult[j] ;
tmp -= idx*_idxMult[j++] ;
if (selDim[j-1]) {
RooAbsLValue* arglv = dynamic_cast<RooAbsLValue*>(arg) ;
theBinVolume *= arglv->getBinWidth(idx) ;
}
}
(*pbinv)[ibin] = theBinVolume ;
}
delete[] selDim ;
_pbinvCacheMgr.setObj(&dimSet,pbinv) ;
_pbinv = pbinv ;
}
Int_t RooDataHist::numEntries() const
{
return RooAbsData::numEntries() ;
}
Double_t RooDataHist::sumEntries() const
{
Int_t i ;
Double_t n(0), carry(0);
for (i=0 ; i<_arrSize ; i++) {
if (!_binValid || _binValid[i]) {
Double_t y = _wgt[i] - carry;
Double_t t = n + y;
carry = (t - n) - y;
n = t;
}
}
return n ;
}
Double_t RooDataHist::sumEntries(const char* cutSpec, const char* cutRange) const
{
checkInit() ;
if (cutSpec==0 && cutRange==0) {
return sumEntries();
} else {
RooFormula* select = 0 ;
if (cutSpec) {
select = new RooFormula("select",cutSpec,*get()) ;
}
Double_t sumw(0), carry(0);
Int_t i ;
for (i=0 ; i<numEntries() ; i++) {
get(i) ;
if (select && select->eval()==0.) continue ;
if (cutRange && !_vars.allInRange(cutRange)) continue ;
if (!_binValid || _binValid[i]) {
Double_t y = weight() - carry;
Double_t t = sumw + y;
carry = (t - sumw) - y;
sumw = t;
}
}
if (select) delete select ;
return sumw ;
}
}
void RooDataHist::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 ;
_cache_sum_valid = kFALSE ;
}
const RooArgSet* RooDataHist::get(Int_t masterIdx) const
{
checkInit() ;
_curWeight = _wgt[masterIdx] ;
_curWgtErrLo = _errLo[masterIdx] ;
_curWgtErrHi = _errHi[masterIdx] ;
_curSumW2 = _sumw2[masterIdx] ;
_curVolume = _binv[masterIdx] ;
_curIndex = masterIdx ;
return RooAbsData::get(masterIdx) ;
}
const RooArgSet* RooDataHist::get(const RooArgSet& coord) const
{
((RooDataHist*)this)->_vars = coord ;
return get(calcTreeIndex()) ;
}
Double_t RooDataHist::binVolume(const RooArgSet& coord)
{
checkInit() ;
((RooDataHist*)this)->_vars = coord ;
return _binv[calcTreeIndex()] ;
}
void RooDataHist::setAllWeights(Double_t value)
{
for (Int_t i=0 ; i<_arrSize ; i++) {
_wgt[i] = value ;
}
_cache_sum_valid = kFALSE ;
}
TIterator* RooDataHist::sliceIterator(RooAbsArg& sliceArg, const RooArgSet& otherArgs)
{
_vars = otherArgs ;
_curIndex = calcTreeIndex() ;
RooAbsArg* intArg = _vars.find(sliceArg) ;
if (!intArg) {
coutE(InputArguments) << "RooDataHist::sliceIterator() variable " << sliceArg.GetName() << " is not part of this RooDataHist" << endl ;
return 0 ;
}
return new RooDataHistSliceIter(*this,*intArg) ;
}
void RooDataHist::SetName(const char *name)
{
if (_dir) _dir->GetList()->Remove(this);
TNamed::SetName(name) ;
if (_dir) _dir->GetList()->Add(this);
}
void RooDataHist::SetNameTitle(const char *name, const char* title)
{
if (_dir) _dir->GetList()->Remove(this);
TNamed::SetNameTitle(name,title) ;
if (_dir) _dir->GetList()->Add(this);
}
void RooDataHist::printValue(ostream& os) const
{
os << numEntries() << " bins (" << sumEntries() << " weights)" ;
}
void RooDataHist::printArgs(ostream& os) const
{
os << "[" ;
_iterator->Reset() ;
RooAbsArg* arg ;
Bool_t first(kTRUE) ;
while((arg=(RooAbsArg*)_iterator->Next())) {
if (first) {
first=kFALSE ;
} else {
os << "," ;
}
os << arg->GetName() ;
}
os << "]" ;
}
void RooDataHist::cacheValidEntries()
{
checkInit() ;
if (!_binValid) {
_binValid = new Bool_t[_arrSize] ;
}
TIterator* iter = _vars.createIterator() ;
RooAbsArg* arg ;
for (Int_t i=0 ; i<_arrSize ; i++) {
get(i) ;
_binValid[i] = kTRUE ;
iter->Reset() ;
while((arg=(RooAbsArg*)iter->Next())) {
_binValid[i] &= arg->inRange(0) ;
}
}
delete iter ;
}
Bool_t RooDataHist::valid() const
{
if (_binValid) {
return _binValid[_curIndex] ;
}
return kTRUE ;
}
Bool_t RooDataHist::isNonPoissonWeighted() const
{
for (int i=0 ; i<numEntries() ; i++) {
if (fabs(_wgt[i]-Int_t(_wgt[i]))>1e-10) return kTRUE ;
}
return kFALSE ;
}
void RooDataHist::printMultiline(ostream& os, Int_t content, Bool_t verbose, TString indent) const
{
RooAbsData::printMultiline(os,content,verbose,indent) ;
os << indent << "Binned Dataset " << GetName() << " (" << GetTitle() << ")" << endl ;
os << indent << " Contains " << numEntries() << " bins with a total weight of " << sumEntries() << endl;
if (!verbose) {
os << indent << " Observables " << _vars << endl ;
} else {
os << indent << " Observables: " ;
_vars.printStream(os,kName|kValue|kExtras|kTitle,kVerbose,indent+" ") ;
}
if(verbose) {
if (_cachedVars.getSize()>0) {
os << indent << " Caches " << _cachedVars << endl ;
}
}
}
void RooDataHist::Streamer(TBuffer &R__b)
{
if (R__b.IsReading()) {
UInt_t R__s, R__c;
Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
if (R__v>2) {
R__b.ReadClassBuffer(RooDataHist::Class(),this,R__v,R__s,R__c);
initialize(0,kFALSE) ;
} else {
UInt_t R__s1, R__c1;
Version_t R__v1 = R__b.ReadVersion(&R__s1, &R__c1); if (R__v1) { }
RooAbsData::Streamer(R__b);
TTree* X_tree(0) ; R__b >> X_tree;
RooArgSet X_truth ; X_truth.Streamer(R__b);
TString X_blindString ; X_blindString.Streamer(R__b);
R__b.CheckByteCount(R__s1, R__c1, RooTreeData::Class());
_dstore = new RooTreeDataStore(X_tree,_vars) ;
_dstore->SetName(GetName()) ;
_dstore->SetTitle(GetTitle()) ;
_dstore->checkInit() ;
RooDirItem::Streamer(R__b);
R__b >> _arrSize;
delete [] _wgt;
_wgt = new Double_t[_arrSize];
R__b.ReadFastArray(_wgt,_arrSize);
delete [] _errLo;
_errLo = new Double_t[_arrSize];
R__b.ReadFastArray(_errLo,_arrSize);
delete [] _errHi;
_errHi = new Double_t[_arrSize];
R__b.ReadFastArray(_errHi,_arrSize);
delete [] _sumw2;
_sumw2 = new Double_t[_arrSize];
R__b.ReadFastArray(_sumw2,_arrSize);
delete [] _binv;
_binv = new Double_t[_arrSize];
R__b.ReadFastArray(_binv,_arrSize);
_realVars.Streamer(R__b);
R__b >> _curWeight;
R__b >> _curWgtErrLo;
R__b >> _curWgtErrHi;
R__b >> _curSumW2;
R__b >> _curVolume;
R__b >> _curIndex;
R__b.CheckByteCount(R__s, R__c, RooDataHist::IsA());
}
} else {
R__b.WriteClassBuffer(RooDataHist::Class(),this);
}
}