Chi^2 of 3D RooFit PDF fit to TH3D

Hello,

I want to get the chi^2 of a 3D fit. I have done this with a 1D fit:

``````        RooRealVar n("n","n",histmin,histmax);
RooDataHist ndata("ndata","ndata",n,test_hist);
RooRealVar fake("fake","fake",fakefraction,0,1);
RooDataHist real_temp("real_temp","real_temp",n,real_hist);
RooHistPdf pdf_real("pdf_real","pdf_real",n,real_temp,0);

RooDataHist fake_temp("fake_temp","fake_temp",n,fake_hist);
RooHistPdf pdf_fake("pdf_fake","pdf_fake",n,fake_temp,0);
RooFitResult *result = pdf.fitTo(ndata,Range(fitmin,fitmax),Save(kTRUE),PrintLevel(1),SumW2Error(kFALSE));
RooPlot *nframe = n.frame();
ndata.plotOn(nframe);
pdf.plotOn(nframe);
int n_param = result->floatParsFinal().getSize();
chi2ndof_out=nframe->chiSquare(n_param);

``````

but this uses a RooPlot frame which is only 1D (as far as I know).
I want to do a 3D fit:

``````        double_t nevents = test_hist->Integral(binx1,binx2,biny1,biny2,binz1,binz2);

RooRealVar nx("nx","nx",histmin[0],histmax[0]);
RooRealVar ny("ny","ny",histmin[1],histmax[1]);
RooRealVar nz("nz","nz",histmin[2],histmax[2]);
RooDataHist ndata("ndata","ndata",RooArgList(nx,ny,nz),test_hist);

RooRealVar fake("fake","fake",fakefraction,0,1);

RooDataHist real_temp("real_temp","real_temp",RooArgList(nx,ny,nz),real_hist);
RooHistPdf pdf_real("pdf_real","pdf_real",RooArgList(nx,ny,nz),RooArgList(nx,ny,nz),real_temp,0);

RooDataHist fake_temp("fake_temp","fake_temp",RooArgList(nx,ny,nz),fake_hist);
RooHistPdf pdf_fake("pdf_fake","pdf_fake",RooArgList(nx,ny,nz),RooArgList(nx,ny,nz),fake_temp,0);

nx.setRange("fitRange",fitmin[0],fitmax[0]);
ny.setRange("fitRange",fitmin[1],fitmax[1]);
nz.setRange("fitRange",fitmin[2],fitmax[2]);

RooFitResult *result = pdf3D.fitTo(ndata,Range("fitRange"),Save(kTRUE),PrintLevel(1),SumW2Error(kFALSE));

``````

And then calculate the chi^2/ndof for the fit. Both the 1D and 3D fits are PDFs from histogram templates.

If there is a roofit function that will do a chi^2 for a 3D fit, that is what I would want. But there doesn’t seem to be one so I am trying to code my own chi^2 and get the same answer for the 1D fit as the `nframe->chiSquare(n_param);` gives (as a sanity check, then apply the same method to the 3D fit).

This has led me to this bit of code: https://root.cern.ch/doc/master/RooCurve_8cxx_source.html#l00545 and I am trying to recreate it:

``````        TH1D *fit_hist = new TH1D("fit_hist","fit_hist",test_hist->GetNbinsX(),histmin,histmax);
RooDataHist *fit_data = pdf.generateBinned( n, Name("fit_data"), Verbose(0), NumEvents(nevents), ExpectedData() );
fit_data->fillHistogram(fit_hist,RooArgList(n),"",0);
chi2ndof_out =0;
int nbins = 0;
double_t x,y,eyl,eyh,avg;
for (int i = test_hist->FindBin(fitmin); i <= (test_hist->FindBin(fitmax)); i++) {
x = test_hist->GetXaxis()->GetBinCenter(i);
y = test_hist->GetBinContent(i);
eyl = test_hist->GetBinErrorLow(i);
eyh = test_hist->GetBinErrorUp(i);
avg = fit_hist->GetBinContent(i);
if(y!=0){
Double_t pull = (y>avg) ? ((y-avg)/eyl) : ((y-avg)/eyh) ;
chi2ndof_out  += pull*pull;
nbins++;
}
}
int ndof = nbins - n_param;
chi2ndof_out = chi2ndof_out/ndof;
``````

But there must be something that I am doing wrong. I get a chi^/ndof of 61.9025 compared to the `nframe->chiSquare(n_param);` gives 3.6192e-13.

I am wondering if it has something to do with the errors of the pdf vs the errors of the generated histogram.
I hope someone can help me see what I am missing.

Best,
Sarah

Hi Sarah,

for plots, you can indeed only compare the Chi^2 of the stuff that’s in the plot. That will give very different Chi^2 values because the PDF is projected onto the plot variable.
You could try to construct a RooChi2Var:
https://root.cern.ch/doc/master/classRooChi2Var.html

That’s what the Chi^2 fit uses.

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.