Hello,
Its been several hours I’m trying to access the chisquare of a fit i’ve done for 2D type of data with roofit. But i can’t find a good way to do it, the chiSquare() function seems to not work if I don’t project my fit, and therefore i’m not really accessing the whole informations.
Here is my code :
#include "TF2.h"
#include "TH2.h"
#include "TMath.h"
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooGaussian.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "TTree.h"
#include "TH1D.h"
#include "TRandom.h"
TH1 *makeTH1(TRandom &trnd);
TTree *makeTTree(TRandom &trnd);
TString DataFile = Form("Root_Files/AnalysisResults_symmetric.root");
TString Table = Form("hf-task-correlation-dplus-dplus-reduced");
//Choose the data we will look into : 1 for Dplus, 2 : for Dminus, 3 : for Dplusminus
int Choose_type = 2;
TString Data_type = Form("");
if (Choose_type == 1) {Data_type = Form("DplusPair");}
else if (Choose_type == 2) {Data_type = Form("DminusPair");}
else if (Choose_type == 3) {Data_type = Form("DplusminusPair");}
else {
cout<<"Please choose a correct choose_type, it should be 1, 2 or 3"<<endl;
exit(0);
}
TString Data_name = Form("hMass");
Data_name = Data_name + Data_type;
TString Save_name = Form("Root_Files/Results");
Save_name = Save_name + Data_type;
Save_name = Save_name + Form(".root");
void Fitting() {
// O p e n t h e d a t a f i l e
// -------------------------------------------------------------------
TFile *myFile = new TFile(DataFile);
TDirectoryFile *Folder = (TDirectoryFile*)myFile->Get(Table);
TH2F *Graph = (TH2F*)Folder->Get(Data_name.Data());
using namespace RooFit;
// P r e p a r a t i o n o f t h e f i t t i n g
// -------------------------------------------------------------------
//Creation of the variables
RooRealVar invmassx("invmassx", "inv_mass_first_candidate (Gev/c²)", 1.80, 1.95);
RooRealVar invmassy("invmassy", "inv_mass_second_candidate (Gev/c²)", 1.80, 1.95);
RooRealVar mean1("mean1", "mean", 1.87112, 1.85, 1.88);
RooRealVar sigma1("sigma1", "sigma", 0.00904307, 0.007, 0.015);
RooRealVar bkgfirst1("bkgfirst1", "background first parameter", 275.332, -100, 1000);
RooRealVar bkgsecond1("bkgsecond1", "background second parameter", 21.1663, -1, 300);
RooRealVar bkgthird1("bkgthird1", "background first parameter", -60.1751, -300, 2);
//Normalisation variable for the additional pdf
RooRealVar nsig("nsig", "#signal events", 800., 0., 3000);
RooRealVar nsigbkg("nsigbkg", "#signal events", 83598., 3000., 1000000);
RooRealVar nbkgsig("nbkgsig", "#signal events", 83598., 3000., 1000000);
RooRealVar nbkg("nbkg", "#signal events", 449101., 0., 1000000);
//Creation of the frame
RooDataHist dh("dh", "dh", {invmassx,invmassy}, Graph);
RooPlot *frame = new RooPlot(invmassx,invmassy,1.80,1.95,1.80,1.95);
dh.plotOn(frame,Project({invmassx,invmassy}), Name("dh"));
//Creation of the function
RooGaussian gaussX("gaussX", "gauss", invmassx, mean1, sigma1);
RooPolynomial backgroundX("backgroundX", "Polynomial", invmassx, {bkgfirst1,bkgsecond1,bkgthird1});
RooGaussian gaussY("gaussY", "gauss", invmassy, mean1, sigma1);
RooPolynomial backgroundY("backgroundY", "Polynomial", invmassy, {bkgfirst1,bkgsecond1,bkgthird1});
// C o n s t r u c t u n c o r r e l a t e d p r o d u c t p d f
// -------------------------------------------------------------------
//Multiply gaussx and gaussy into a two-dimensional p.d.f. gaussxy
RooProdPdf gaussxy("gaussxy","gaussX*gaussY",RooArgList(gaussX,gaussY));
RooProdPdf gaussxbkgy("gaussxbkgy","gaussX*backgroundY",RooArgList(gaussX,backgroundY));
RooProdPdf bkgxgaussy("bkgxgaussy","backgroundX*gaussY",RooArgList(backgroundX,gaussY));
RooProdPdf bkgxy("bkgxy","backgroundX*backgroundY",RooArgList(backgroundX,backgroundY));
//Addition of all the functions
RooAddPdf total("total", "g+a", {gaussxy,gaussxbkgy,bkgxgaussy,bkgxy}, {nsig, nsigbkg, nbkgsig, nbkg});
RooFitResult *fit_result_data = total.fitTo(dh, Save());
total.plotOn(frame,Project({invmassx,invmassy}), Name("model_curve"));
TCanvas *c = new TCanvas("rf304_uncorrprod","rf304_uncorrprod",800, 400);
Graph->GetXaxis()->SetRangeUser(1.8, 1.95);
Graph->GetYaxis()->SetRangeUser(1.8, 1.95);
// I n t e g r a t i o n o f t h e d o u b l e g a u s s i a n
// -------------------------------------------------------------------
invmassx.setRange("range",mean1.getVal()-sigma1.getVal()*3,mean1.getVal()+sigma1.getVal()*3);
invmassy.setRange("range",mean1.getVal()-sigma1.getVal()*3,mean1.getVal()+sigma1.getVal()*3);
RooAbsReal *Integral = gaussxy.createIntegral(RooArgSet(invmassx,invmassy),NormSet(RooArgSet(invmassx,invmassy)),Range("range"));
// G e t t h e n u m b e r o f m a t c h e d
// -------------------------------------------------------------------
TString Matched_name = Data_name + "sMatched";
TH2F *Matched_graph = (TH2F*)Folder->Get(Matched_name.Data());
int Matched_value = Matched_graph->GetEntries();
// C r e a t i o n o f p r o j e c t i o n t o a c c e s s t h e c h i s q u a r e
// -------------------------------------------------------------------
RooPlot *xframe = invmassx.frame(Title("Projection on x"));
dh.plotOn(xframe);
total.plotOn(xframe);
RooPlot *yframe = invmassy.frame(Title("Projection on y"));
dh.plotOn(yframe);
total.plotOn(yframe);
// P r e s e n t a t i o n o f r e s u l t s
// -------------------------------------------------------------------
TH2* pdf_2d = (TH2*) total.createHistogram("invmassx,invmassy",Graph->FindBin(1.95,1.80)-Graph->FindBin(1.80,1.80),Graph->FindBin(1.95,1.80)-Graph->FindBin(1.80,1.80)) ;
pdf_2d->SetLineColor(2);
Graph->Draw("lego2");
pdf_2d->Draw("surf same");
c->Print(Save_name);
cout<<endl<<endl<<endl<<endl<<"E n d r e s u l t s :"<<endl<<endl<<endl<<endl;
cout<<"Number of events :"<<pdf_2d->Integral(0,-1,0,-1)<<endl;
cout<<"Sigma :"<<sigma1.getVal()<<"+-"<<sigma1.getError()<<endl;
cout<<"Mean :"<<mean1.getVal()<<"+-"<<mean1.getError()<<endl;
cout<<"Normalisation factor :"<<nsig.getVal()<<"+-"<<nsig.getError()<<endl;
cout<<"MonteCarlo truth : "<<Matched_value<<endl;
cout<<"Raw yield :"<<nsig.getVal()*Integral->getVal()<<"+-"<<nsig.getError()*Integral->getVal()+nsig.getVal()*Integral->getPropagatedError(*fit_result_data,{invmassx,invmassy})+nsig.getError()*Integral->getPropagatedError(*fit_result_data,{invmassx,invmassy})<<endl;
cout<<"Efficiency : "<<nsig.getVal()*Integral->getVal()/Matched_value<<"+-"<<(nsig.getError()*Integral->getVal()+nsig.getVal()*Integral->getPropagatedError(*fit_result_data,{invmassx,invmassy})+nsig.getError()*Integral->getPropagatedError(*fit_result_data,{invmassx,invmassy}))/Matched_value<<endl;
cout<<"chisquare for x projection :"<<xframe->chiSquare()<<endl;
cout<<"chisquare for y projection :"<<yframe->chiSquare()<<endl;
}