Helllo,
I have a distribution that I fit with a function in some sub-range (“sideband”) of a histogrammed distribution. I then use the fitting function to extrapolate to a region of interest (“extrapolate”). I would like to evaluate the integral of the fitting function with its uncertainty, as determined by the parameter errors returned by the fit.
I know how to do this in ROOT (TF1::IntegralError), but not in RooFit. I naively tried using an extended-PDF and taking the fraction of the distribution in the region of interest times the uncertainty on the normalization constant. However, this clearly does not take into account the uncertainty in the shape of the fitting function, and underestimates the overall integral uncertainty.
I am attaching:
o a macro that demonstrates the problem.
o a picture showing the output of the macro (ROOT fit on the left, RooFit plot on the right, yellow box with numerical results at the bottom).
Any help would be very much appreciated.
–Christos
PS An aside: is there a simple way to get the RooFit stats box to look like the ROOT one?
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooExtendPdf.h"
#include "RooFitResult.h"
#include "RooDataHist.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "TH1F.h"
#include "TF1.h"
#include "TRandom.h"
#include "TFitResult.h"
#include "TMatrixDSym.h"
#include "RooPlot.h"
using namespace RooFit ;
void test()
{
const float XMIN = -10; const float XMAX = 10;
const float sideband_min = -2.0;
const float sideband_max = 2.0;
const float signal_min = 2.0;
const float signal_max = XMAX;
const float MEAN_INPUT = 0.0;
const float SIGMA_INPUT = 1.0;
// Declare observable x
RooRealVar x("x","x",XMIN,XMAX) ;
// Create Gaussian PDFs g(x,mean,sigma)
RooRealVar mean("mean","mean of gaussian",MEAN_INPUT, -1000, 1000) ;
RooRealVar sigma("sigma","width of gaussians",SIGMA_INPUT, 0, 1000) ;
RooGaussian sig("sig","gaussian",x,mean,sigma);
TH1F * h = new TH1F("h", "ROOT fit", 100, XMIN, XMAX);
for(int i = 0; i != 1000; ++i)
h->Fill(gRandom->Gaus(MEAN_INPUT, SIGMA_INPUT));
// Define signal range in which events counts are to be defined
x.setRange("sideband", sideband_min, sideband_max);
x.setRange("extrapolate", signal_min, signal_max);
x.setRange("fullRange", XMIN, XMAX);
// Associated nsig as expected number of events
RooRealVar nsig("nsig","number of signal events",50,0.,10000) ;
RooExtendPdf esig("esig","extended signal p.d.f",sig,nsig,"fullRange") ;
RooDataHist h2("h2", "h2", x, Import(*h));
TF1 *f1 = new TF1("f1", "gaus", sideband_min, sideband_max);
TFitResultPtr r = h->Fit("f1", "LRS");
Double_t integ = f1->Integral(signal_min, signal_max, r->GetParams())
/ h->GetBinWidth(0);
Double_t dinteg = f1->IntegralError(signal_min, signal_max, r->GetParams(),
r->GetCovarianceMatrix().GetMatrixArray())/ h->GetBinWidth(0);
// Perform unbinned extended ML fit to data
RooFitResult * rf = esig.fitTo(h2,Extended(kTRUE), Range("sideband"), Save()) ;
RooAbsReal* xtra_int = esig.createIntegral(x,NormSet(x), Range("extrapolate")) ;
RooAbsReal* full_int = esig.createIntegral(x,NormSet(x), Range("fullRange")) ;
// cout << " xtra_int = " << xtra_int->getVal() << endl;
// cout << " full_int = " << full_int->getVal() << endl;
RooArgList pars = rf->floatParsFinal();
float n = ((RooRealVar*) pars.find("nsig"))->getVal();
float dn = ((RooRealVar*) pars.find("nsig"))->getError();
RooPlot * xframe = x.frame(Title("RooFit fit"));
h2.plotOn(xframe);
esig.plotOn(xframe);
esig.paramOn(xframe,Layout(0.55));
TCanvas * c1 = new TCanvas();
xframe->Draw();
cout << ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << endl;
int min_bin = h->FindBin(signal_min);
int max_bin = h->FindBin(signal_max);
cout << " Actual # of entries between [" << signal_min << ", " << signal_max
<< "] = " << h->Integral(min_bin, max_bin) << endl;
cout << " Extrapolation between [" << signal_min << ", " << signal_max
<< "] with ROOT, integral = " << integ << " +- " << dinteg << endl;
cout << " Extrapolation between [" << signal_min << ", " << signal_max
<< "] with RooFit, integral = "
<< xtra_int->getVal() * n/full_int->getVal() << " +- "
<< xtra_int->getVal() * dn/full_int->getVal() << endl;
cout << ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << endl;
}