// A RooHist is a graphical representation of binned data based on the
// TGraphAsymmErrors class. Error bars are calculated using either Poisson
// or Binomial statistics. A RooHist is used to represent histograms in
// a RooPlot.
// END_HTML
#include "RooFit.h"
#include "RooHist.h"
#include "RooHist.h"
#include "RooHistError.h"
#include "RooCurve.h"
#include "RooMsgService.h"
#include "TH1.h"
#include "TClass.h"
#include "Riostream.h"
#include <iomanip>
#include <math.h>
ClassImp(RooHist)
;
RooHist::RooHist() :
_nominalBinWidth(1),
_nSigma(1),
_entries(0),
_rawEntries(0)
{
}
RooHist::RooHist(Double_t nominalBinWidth, Double_t nSigma, Double_t , Double_t ) :
TGraphAsymmErrors(), _nominalBinWidth(nominalBinWidth), _nSigma(nSigma), _rawEntries(-1)
{
initialize();
}
RooHist::RooHist(const TH1 &data, Double_t nominalBinWidth, Double_t nSigma, RooAbsData::ErrorType etype, Double_t xErrorFrac,
Bool_t correctForBinWidth, Double_t scaleFactor) :
TGraphAsymmErrors(), _nominalBinWidth(nominalBinWidth), _nSigma(nSigma), _rawEntries(-1)
{
initialize();
SetName(data.GetName());
SetTitle(data.GetTitle());
if(_nominalBinWidth == 0) {
const TAxis *axis= ((TH1&)data).GetXaxis();
if(axis->GetNbins() > 0) _nominalBinWidth= (axis->GetXmax() - axis->GetXmin())/axis->GetNbins();
}
setYAxisLabel(const_cast<TH1&>(data).GetYaxis()->GetTitle());
Int_t nbin= data.GetNbinsX();
for(Int_t bin= 1; bin <= nbin; bin++) {
Axis_t x= data.GetBinCenter(bin);
Stat_t y= data.GetBinContent(bin);
Stat_t dy = data.GetBinError(bin) ;
if (etype==RooAbsData::Poisson) {
addBin(x,y,data.GetBinWidth(bin),xErrorFrac,scaleFactor);
} else if (etype==RooAbsData::SumW2) {
addBinWithError(x,y,dy,dy,data.GetBinWidth(bin),xErrorFrac,correctForBinWidth,scaleFactor);
} else {
addBinWithError(x,y,0,0,data.GetBinWidth(bin),xErrorFrac,correctForBinWidth,scaleFactor);
}
}
_entries+= data.GetBinContent(0) + data.GetBinContent(nbin+1);
}
RooHist::RooHist(const TH1 &data1, const TH1 &data2, Double_t nominalBinWidth, Double_t nSigma,
RooAbsData::ErrorType etype, Double_t xErrorFrac, Bool_t efficiency, Double_t scaleFactor) :
TGraphAsymmErrors(), _nominalBinWidth(nominalBinWidth), _nSigma(nSigma), _rawEntries(-1)
{
initialize();
SetName(data1.GetName());
SetTitle(data1.GetTitle());
if(_nominalBinWidth == 0) {
const TAxis *axis= ((TH1&)data1).GetXaxis();
if(axis->GetNbins() > 0) _nominalBinWidth= (axis->GetXmax() - axis->GetXmin())/axis->GetNbins();
}
if (!efficiency) {
setYAxisLabel(Form("Asymmetry (%s - %s)/(%s + %s)",
data1.GetName(),data2.GetName(),data1.GetName(),data2.GetName()));
} else {
setYAxisLabel(Form("Efficiency (%s)/(%s + %s)",
data1.GetName(),data1.GetName(),data2.GetName()));
}
Int_t nbin= data1.GetNbinsX();
if(data2.GetNbinsX() != nbin) {
coutE(InputArguments) << "RooHist::RooHist: histograms have different number of bins" << endl;
return;
}
for(Int_t bin= 1; bin <= nbin; bin++) {
Axis_t x= data1.GetBinCenter(bin);
if(fabs(data2.GetBinCenter(bin)-x)>1e-10) {
coutW(InputArguments) << "RooHist::RooHist: histograms have different centers for bin " << bin << endl;
}
Stat_t y1= data1.GetBinContent(bin);
Stat_t y2= data2.GetBinContent(bin);
if (!efficiency) {
if (etype==RooAbsData::Poisson) {
addAsymmetryBin(x,roundBin(y1),roundBin(y2),data1.GetBinWidth(bin),xErrorFrac,scaleFactor);
} else if (etype==RooAbsData::SumW2) {
Stat_t dy1= data1.GetBinError(bin);
Stat_t dy2= data2.GetBinError(bin);
addAsymmetryBinWithError(x,y1,y2,dy1,dy2,data1.GetBinWidth(bin),xErrorFrac,scaleFactor);
} else {
addAsymmetryBinWithError(x,y1,y2,0,0,data1.GetBinWidth(bin),xErrorFrac,scaleFactor);
}
} else {
if (etype==RooAbsData::Poisson) {
addEfficiencyBin(x,roundBin(y1),roundBin(y2),data1.GetBinWidth(bin),xErrorFrac,scaleFactor);
} else if (etype==RooAbsData::SumW2) {
Stat_t dy1= data1.GetBinError(bin);
Stat_t dy2= data2.GetBinError(bin);
addEfficiencyBinWithError(x,y1,y2,dy1,dy2,data1.GetBinWidth(bin),xErrorFrac,scaleFactor);
} else {
addEfficiencyBinWithError(x,y1,y2,0,0,data1.GetBinWidth(bin),xErrorFrac,scaleFactor);
}
}
}
_entries= -1;
}
RooHist::RooHist(const RooHist& hist1, const RooHist& hist2, Double_t wgt1, Double_t wgt2,
RooAbsData::ErrorType etype, Double_t xErrorFrac) : _rawEntries(-1)
{
initialize() ;
SetName(hist1.GetName()) ;
SetTitle(hist1.GetTitle()) ;
_nominalBinWidth=hist1._nominalBinWidth ;
_nSigma=hist1._nSigma ;
setYAxisLabel(hist1.getYAxisLabel()) ;
if (!hist1.hasIdenticalBinning(hist2)) {
coutE(InputArguments) << "RooHist::RooHist input histograms have incompatible binning, combined histogram will remain empty" << endl ;
return ;
}
if (etype==RooAbsData::Poisson) {
if (wgt1!=1.0 || wgt2 != 1.0) {
coutW(InputArguments) << "RooHist::RooHist: WARNING: Poisson errors of weighted sum of two histograms is not well defined! " << endl
<< " Summed histogram bins will rounded to nearest integer for Poisson confidence interval calculation" << endl ;
}
Int_t i,n=hist1.GetN() ;
for(i=0 ; i<n ; i++) {
Double_t x1,y1,x2,y2,dx1 ;
#if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,1)
hist1.GetPoint(i,x1,y1) ;
#else
const_cast<RooHist&>(hist1).GetPoint(i,x1,y1) ;
#endif
dx1 = hist1.GetErrorX(i) ;
#if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,1)
hist2.GetPoint(i,x2,y2) ;
#else
const_cast<RooHist&>(hist2).GetPoint(i,x2,y2) ;
#endif
addBin(x1,roundBin(wgt1*y1+wgt2*y2),2*dx1/xErrorFrac,xErrorFrac) ;
}
} else {
Int_t i,n=hist1.GetN() ;
for(i=0 ; i<n ; i++) {
Double_t x1,y1,x2,y2,dx1,dy1,dy2 ;
#if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,1)
hist1.GetPoint(i,x1,y1) ;
#else
const_cast<RooHist&>(hist1).GetPoint(i,x1,y1) ;
#endif
dx1 = hist1.GetErrorX(i) ;
dy1 = hist1.GetErrorY(i) ;
dy2 = hist2.GetErrorY(i) ;
#if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,1)
hist2.GetPoint(i,x2,y2) ;
#else
const_cast<RooHist&>(hist2).GetPoint(i,x2,y2) ;
#endif
Double_t dy = sqrt(wgt1*wgt1*dy1*dy1+wgt2*wgt2*dy2*dy2) ;
addBinWithError(x1,wgt1*y1+wgt2*y2,dy,dy,2*dx1/xErrorFrac,xErrorFrac) ;
}
}
}
void RooHist::initialize()
{
SetMarkerStyle(8);
_entries= 0;
}
Double_t RooHist::getFitRangeNEvt() const
{
return (_rawEntries==-1 ? _entries : _rawEntries) ;
}
Double_t RooHist::getFitRangeNEvt(Double_t xlo, Double_t xhi) const
{
Double_t sum(0) ;
for (int i=0 ; i<GetN() ; i++) {
Double_t x,y ;
#if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,1)
GetPoint(i,x,y) ;
#else
const_cast<RooHist*>(this)->GetPoint(i,x,y) ;
#endif
if (x>=xlo && x<=xhi) {
sum += y ;
}
}
if (_rawEntries!=-1) {
coutW(Plotting) << "RooHist::getFitRangeNEvt() WARNING: Number of normalization events associated to histogram is not equal to number of events in histogram" << endl
<< " due cut made in RooAbsData::plotOn() call. Automatic normalization over sub-range of plot variable assumes" << endl
<< " that the effect of that cut is uniform across the plot, which may be an incorrect assumption. To be sure of" << endl
<< " correct normalization explicit pass normalization information to RooAbsPdf::plotOn() call using Normalization()" << endl ;
sum *= _rawEntries / _entries ;
}
return sum ;
}
Double_t RooHist::getFitRangeBinW() const
{
return _nominalBinWidth ;
}
Int_t RooHist::roundBin(Double_t y)
{
if(y < 0) {
coutW(Plotting) << fName << "::roundBin: rounding negative bin contents to zero: " << y << endl;
return 0;
}
Int_t n= (Int_t)(y+0.5);
if(fabs(y-n)>1e-6) {
coutW(Plotting) << fName << "::roundBin: rounding non-integer bin contents: " << y << endl;
}
return n;
}
void RooHist::addBin(Axis_t binCenter, Double_t n, Double_t binWidth, Double_t xErrorFrac, Double_t scaleFactor)
{
if (n<0) {
coutW(Plotting) << "RooHist::addBin(" << GetName() << ") WARNING: negative entry set to zero when Poisson error bars are requested" << endl ;
}
Double_t scale= 1;
if(binWidth > 0) {
scale= _nominalBinWidth/binWidth;
}
_entries+= n;
Int_t index= GetN();
Double_t ym,yp,dx(0.5*binWidth);
if (fabs((double)((n-Int_t(n))>1e-5))) {
Double_t ym1,yp1,ym2,yp2 ;
Int_t n1 = Int_t(n) ;
Int_t n2 = n1+1 ;
if(!RooHistError::instance().getPoissonInterval(n1,ym1,yp1,_nSigma) ||
!RooHistError::instance().getPoissonInterval(n2,ym2,yp2,_nSigma)) {
coutE(Plotting) << "RooHist::addBin: unable to add bin with " << n << " events" << endl;
}
ym = ym1 + (n-n1)*(ym2-ym1) ;
yp = yp1 + (n-n1)*(yp2-yp1) ;
coutW(Plotting) << "RooHist::addBin(" << GetName()
<< ") WARNING: non-integer bin entry " << n << " with Poisson errors, interpolating between Poisson errors of adjacent integer" << endl ;
} else {
if(!RooHistError::instance().getPoissonInterval(Int_t(n),ym,yp,_nSigma)) {
coutE(Plotting) << "RooHist::addBin: unable to add bin with " << n << " events" << endl;
return;
}
}
SetPoint(index,binCenter,n*scale*scaleFactor);
SetPointError(index,dx*xErrorFrac,dx*xErrorFrac,scale*(n-ym)*scaleFactor,scale*(yp-n)*scaleFactor);
updateYAxisLimits(scale*yp);
updateYAxisLimits(scale*ym);
}
void RooHist::addBinWithError(Axis_t binCenter, Double_t n, Double_t elow, Double_t ehigh, Double_t binWidth,
Double_t xErrorFrac, Bool_t correctForBinWidth, Double_t scaleFactor)
{
Double_t scale= 1;
if(binWidth > 0 && correctForBinWidth) {
scale= _nominalBinWidth/binWidth;
}
_entries+= n;
Int_t index= GetN();
Double_t dx(0.5*binWidth) ;
SetPoint(index,binCenter,n*scale*scaleFactor);
SetPointError(index,dx*xErrorFrac,dx*xErrorFrac,elow*scale*scaleFactor,ehigh*scale*scaleFactor);
updateYAxisLimits(scale*(n-elow));
updateYAxisLimits(scale*(n+ehigh));
}
void RooHist::addBinWithXYError(Axis_t binCenter, Double_t n, Double_t exlow, Double_t exhigh, Double_t eylow, Double_t eyhigh,
Double_t scaleFactor)
{
_entries+= n;
Int_t index= GetN();
SetPoint(index,binCenter,n*scaleFactor);
SetPointError(index,exlow,exhigh,eylow*scaleFactor,eyhigh*scaleFactor);
updateYAxisLimits(scaleFactor*(n-eylow));
updateYAxisLimits(scaleFactor*(n+eyhigh));
}
void RooHist::addAsymmetryBin(Axis_t binCenter, Int_t n1, Int_t n2, Double_t binWidth, Double_t xErrorFrac, Double_t scaleFactor)
{
Double_t scale= 1;
if(binWidth > 0) scale= _nominalBinWidth/binWidth;
Int_t index= GetN();
Double_t ym,yp,dx(0.5*binWidth);
if(!RooHistError::instance().getBinomialIntervalAsym(n1,n2,ym,yp,_nSigma)) {
coutE(Plotting) << "RooHist::addAsymmetryBin: unable to calculate binomial error for bin with " << n1 << "," << n2 << " events" << endl;
return;
}
Double_t a= (Double_t)(n1-n2)/(n1+n2);
SetPoint(index,binCenter,a*scaleFactor);
SetPointError(index,dx*xErrorFrac,dx*xErrorFrac,(a-ym)*scaleFactor,(yp-a)*scaleFactor);
updateYAxisLimits(scale*yp);
updateYAxisLimits(scale*ym);
}
void RooHist::addAsymmetryBinWithError(Axis_t binCenter, Double_t n1, Double_t n2, Double_t en1, Double_t en2, Double_t binWidth, Double_t xErrorFrac, Double_t scaleFactor)
{
Double_t scale= 1;
if(binWidth > 0) scale= _nominalBinWidth/binWidth;
Int_t index= GetN();
Double_t ym,yp,dx(0.5*binWidth);
Double_t a= (Double_t)(n1-n2)/(n1+n2);
Double_t error = 2*sqrt( pow(en1,2)*pow(n2,2) + pow(en2,2)*pow(n1,2) ) / pow(n1+n2,2) ;
ym=a-error ;
yp=a+error ;
SetPoint(index,binCenter,a*scaleFactor);
SetPointError(index,dx*xErrorFrac,dx*xErrorFrac,(a-ym)*scaleFactor,(yp-a)*scaleFactor);
updateYAxisLimits(scale*yp);
updateYAxisLimits(scale*ym);
}
void RooHist::addEfficiencyBin(Axis_t binCenter, Int_t n1, Int_t n2, Double_t binWidth, Double_t xErrorFrac, Double_t scaleFactor)
{
Double_t scale= 1;
if(binWidth > 0) scale= _nominalBinWidth/binWidth;
Int_t index= GetN();
Double_t a= (Double_t)(n1)/(n1+n2);
Double_t ym,yp,dx(0.5*binWidth);
if(!RooHistError::instance().getBinomialIntervalEff(n1,n2,ym,yp,_nSigma)) {
coutE(Plotting) << "RooHist::addEfficiencyBin: unable to calculate binomial error for bin with " << n1 << "," << n2 << " events" << endl;
return;
}
SetPoint(index,binCenter,a*scaleFactor);
SetPointError(index,dx*xErrorFrac,dx*xErrorFrac,(a-ym)*scaleFactor,(yp-a)*scaleFactor);
updateYAxisLimits(scale*yp);
updateYAxisLimits(scale*ym);
}
void RooHist::addEfficiencyBinWithError(Axis_t binCenter, Double_t n1, Double_t n2, Double_t en1, Double_t en2, Double_t binWidth, Double_t xErrorFrac, Double_t scaleFactor)
{
Double_t scale= 1;
if(binWidth > 0) scale= _nominalBinWidth/binWidth;
Int_t index= GetN();
Double_t a= (Double_t)(n1)/(n1+n2);
Double_t error = sqrt( pow(en1,2)*pow(n2,2) + pow(en2,2)*pow(n1,2) ) / pow(n1+n2,2) ;
Double_t ym,yp,dx(0.5*binWidth);
ym=a-error ;
yp=a+error ;
SetPoint(index,binCenter,a*scaleFactor);
SetPointError(index,dx*xErrorFrac,dx*xErrorFrac,(a-ym)*scaleFactor,(yp-a)*scaleFactor);
updateYAxisLimits(scale*yp);
updateYAxisLimits(scale*ym);
}
RooHist::~RooHist()
{
}
Bool_t RooHist::hasIdenticalBinning(const RooHist& other) const
{
if (GetN() != other.GetN()) {
return kFALSE ;
}
Int_t i ;
for (i=0 ; i<GetN() ; i++) {
Double_t x1,x2,y1,y2 ;
#if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,1)
GetPoint(i,x1,y1) ;
other.GetPoint(i,x2,y2) ;
#else
const_cast<RooHist&>(*this).GetPoint(i,x1,y1) ;
const_cast<RooHist&>(other).GetPoint(i,x2,y2) ;
#endif
if (fabs(x1-x2)>1e-10) {
return kFALSE ;
}
}
return kTRUE ;
}
Bool_t RooHist::isIdentical(const RooHist& other, Double_t tol) const
{
TH1::AddDirectory(kFALSE) ;
TH1F h_self("h_self","h_self",GetN(),0,1) ;
TH1F h_other("h_other","h_other",GetN(),0,1) ;
TH1::AddDirectory(kTRUE) ;
for (Int_t i=0 ; i<GetN() ; i++) {
h_self.SetBinContent(i+1,GetY()[i]) ;
h_other.SetBinContent(i+1,other.GetY()[i]) ;
}
Double_t M = h_self.KolmogorovTest(&h_other,"M") ;
if (M>tol) {
Double_t kprob = h_self.KolmogorovTest(&h_other) ;
cout << "RooHist::isIdentical() tolerance exceeded M=" << M << " (tol=" << tol << "), corresponding prob = " << kprob << endl ;
return kFALSE ;
}
return kTRUE ;
}
void RooHist::printMultiline(ostream& os, Int_t contents, Bool_t verbose, TString indent) const
{
RooPlotable::printMultiline(os,contents,verbose,indent);
os << indent << "--- RooHist ---" << endl;
Int_t n= GetN();
os << indent << " Contains " << n << " bins" << endl;
if(verbose) {
os << indent << " Errors calculated at" << _nSigma << "-sigma CL" << endl;
os << indent << " Bin Contents:" << endl;
for(Int_t i= 0; i < n; i++) {
os << indent << setw(3) << i << ") x= " << fX[i];
if(fEXhigh[i] > 0 || fEXlow[i] > 0) {
os << " +" << fEXhigh[i] << " -" << fEXlow[i];
}
os << " , y = " << fY[i] << " +" << fEYhigh[i] << " -" << fEYlow[i] << endl;
}
}
}
void RooHist::printName(ostream& os) const
{
os << GetName() ;
}
void RooHist::printTitle(ostream& os) const
{
os << GetTitle() ;
}
void RooHist::printClassName(ostream& os) const
{
os << IsA()->GetName() ;
}
RooHist* RooHist::makeResidHist(const RooCurve& curve,bool normalize) const
{
RooHist* hist = new RooHist(_nominalBinWidth) ;
hist->SetName(Form(normalize?"pull_%s_s":"resid_%s_s",GetName(),curve.GetName())) ;
hist->SetTitle(Form(normalize?"Pull of %s and %s":"Residual of %s and %s",GetTitle(),curve.GetTitle())) ;
Double_t xstart,xstop,y ;
#if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,1)
curve.GetPoint(0,xstart,y) ;
curve.GetPoint(curve.GetN()-1,xstop,y) ;
#else
const_cast<RooCurve&>(curve).GetPoint(0,xstart,y) ;
const_cast<RooCurve&>(curve).GetPoint(curve.GetN()-1,xstop,y) ;
#endif
for(Int_t i=0 ; i<GetN() ; i++) {
Double_t x,point;
#if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,1)
GetPoint(i,x,point) ;
#else
const_cast<RooHist&>(*this).GetPoint(i,x,point) ;
#endif
if (x<xstart || x>xstop) continue ;
Double_t yy = point - curve.interpolate(x) ;
Double_t dyl = GetErrorYlow(i) ;
Double_t dyh = GetErrorYhigh(i) ;
if (normalize) {
Double_t norm = (yy>0?dyl:dyh);
if (norm==0.) {
coutW(Plotting) << "RooHist::makeResisHist(" << GetName() << ") WARNING: point " << i << " has zero error, setting residual to zero" << endl ;
yy=0 ;
dyh=0 ;
dyl=0 ;
} else {
yy /= norm;
dyh /= norm;
dyl /= norm;
}
}
hist->addBinWithError(x,yy,dyl,dyh);
}
return hist ;
}