RooFit accessing p.d.f's

Roofit and RooStats
1

Hi All,

I’m trying out the new RooFit, and have a question about separating the pdfs.
Here is my situation. I fit a distribution with a background plus gaussian function. The background is modeled after the simulation histogram.

RooRealVar x(“x”,“x”,1.4,2.);
RooRealVar mean(“mean”,“mean”,1.4,1.6);
RooRealVar sigma(“sigma”,“sigma”,0.0005,0.018);

RooGaussian gauss(“gauss”,“Theta Peak”,x,mean,sigma);

RooDataHist simulation(“simulation”,"",x,HSIM);
RooDataHist data(“data”,"",x,HDATA);

RooRealVar nsig(“nsig”,“signal Fracktion”,500,0.,10000.);
RooRealVar nbkg(“nbkg”,“backgroun Fracktion”,500,0.,10000.);

RooHistPdf BackGround(“BackGround”,“Phi Simulation”,x,simulation,7);
RooAddPdf fModel(“Model”,“Model”,RooArgList(BackGround,gauss),RooArgList(nbkg,nsig));

RooFitResult* res = fModel.fitTo(data,Range(1.40,1.80),Minos(kTRUE),PrintLevel(3),Save());

Everything works well, and plots well. Now I want to get a function or a PDF that can give me the value of the background at any point along X for fitted value of course. The final goal is comparing the log-likelihood ration for NULL hypothesis with the FULL hypothesis.

How can I extract the background and the gauss functions from my fitted model ??

_g

2

Hi,

If you have access to the BackGround and gauss pointers, those are the fitted
signal and background component. So with your example macro it would be
sufficient to do this, after you have completed the fit

RooAbsReal* nll_full = fModel.createNLL(data,Range(1.40,1.80))
RooAbsReal* nll_null = BackGround.createNLL(data,Range(1.40,1.80))

and then complare nll_full->getVal() with nll_null->getVal() ;

Another way to do this would be to create a profile likelihood in nsig, e.g.

RooAbsReal* profile_nsig = nll_full->createProfile(nsig) ;

and plot the profile likelihood

RooPlot* frame = nsig.frame(0,2000) ;
profile_nsig->plotOn(frame) ;

which will plot the -log(L) minimized w.r.t all parameters except nsig, minus
the -log(L) minimized w.r.t. all parameters. The value of the profile likelihood
at Nsig=0 is in effect equal to the ‘nll_full->getVal() - nll_null->getVal()’ expression
above. [ Note that a profile likelihood object is expensive: each evaluation requires
a minuit minimization cycle ]

See also tutorial macro root.cern.ch/root/html/tutorials … ell.C.html
on profile likelihood.

A separate question is how to obtain the Background and gauss pointers given
only a pointer to the top-level fModel object. [ NB: In the above profile likelihood approach
you don’t need this info ]. You can do this by calling calling RooAddPdf::pdfList()
which returns the list of pdfs added in your model. More generally you can also call
getComponents() method which will return a RooArgSet all component pdfs
(at any level in the expression tree), which you can use to e.g. look up a component by name. This implies you need to think of some naming scheme for your pdfs.

Wouter

3

Thanks,
It worked, but raise more questions. The nlls give me:

nll_full = -27315.195818
nll_null = 0.000000

I’m not sure this is right. The thing is that I have a peak on top of a background and trying to estimate the significance of the peak, and I would be happy to get the functional form of the background and the peak to also be able to calculate the significance.

When I plot the BackGround on the frame it looks nice, but the gauss function the amplitude comes out very large, larger than it’s fitted value, I can’t seem to make the gaussian plot normally scaled. If I understand it correctly the fitter fits the gaussian parameters and also the nsig number, now when I plot the gauss function, don’t I have to take into account the nsig ? This is what confuses me.

Second part of my question is: is it possible to evaluate a PDF at given x-values ?

4

Hi,

The nll_null seems wrong. If you send me a bit more details, or a tarball with
a macro and the neded input histograms so that I can reproduce your problem
i might be able to say more.

Concerning your nsig question: when you make an extended ML fit of

Nsig* Fsig + Nbkg*Fbkg

you make a plain likelihood fit with Nsig/(Nsig+Nbkg)*Fsig + Nbkg/(Nsig+Nkg)*Fbkg
as probability density function and add as extended term a Poisson with Nobs = Ndata
and Nexp = Nsig+Nbkg.

Whenever you plot a pdf over the data, it will plot the pdf shape defined by the above
fractions and multiply this unit-normalized function with the number of events in the
data that was previously plotted on the frame.

Answer to second question: Yes, that is very easy; you do x.setVal(something)
and then do pdf.geVal(x), where the argument in geVal() specified that the return
value should be normalized w.r.t. observable x.

Wouter

5

Dear Wouter,

I have composed a script that essentially reproduces my problem. I’m able to obtain very nice fits to my data, which are similar to what I programmed into the example, but the problem comes when I try to analyze obtained results.

Please find the script in the Attachment, I would really appreciate help with this, I think I’m making a very simple mistake.

_g
rooFit_Test.C (2.75 KB)