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);
RooAddPdf pdf("pdf","pdf",RooArgList(pdf_fake,pdf_real),fake);
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);
RooAddPdf pdf3D("pdf3D","pdf3D",RooArgList(pdf_fake,pdf_real),fake);
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