RooFit: binned poisson likelihood with nuisance parameters

[stephan.zimmer]

Hi ROOTers & RooFitters,
I’m desperately trying to make RooFit do what I want but can’t seem to figure out how to do this. Here’s my problem:

GOAL:

• I have a measured spectrum, which I want to compare to a theoretical spectrum and determine the best fit values, so my likelihood is typically written like this: L( model | data ) = Prod_i ( Poisson(n_obs_i | n_pred_i) )
• where n_obs_i is the observed number of events in bin i, and n_pred is the number of predicted events. Part of the complication arises from the computation of n_pred_i: n_pred_i = acceptance(E_i) * int_e_low,i^e_high,i { flux(E) dE ) with flux(E) being my model under test, containing a set of parameters of interest (e.g. spectral index, normalization)
• Note however, that in order to compute the predicted events in bin i, I have to integrate over the bin edges
• The problem is further complicated by the fact that I have to introduce a nuisance parameter to account for the systematic uncertainty sigma_syst which varies from bin to bin, so the likelihood becomes: L( model | data ) = Prod_i ( Poisson(n_obs_i | n_pred_i * (1+a*sigma_syst,i) ) x Gaus(0|a,1), where Gaus(0|a,1) is the constraint term/ side-band measurement, but taken from the data directly

MY PROBLEMS:

1. how to construct NLL?

Given the structure of the problem, I’ve been wondering if I could construct the pdf iteratively, by bin-wise defining the Poisson(n_obs,i | n_pred_i * (1+a*sigma_syst,i) ) which allows me to properly account for the bin-wise change of sigma_syst, the pdf would then be

binList = RooArgList([... list of poisson terms ... + Gaussian constraint term])
pdf = RooProdPdf("pdf","myPdf",binList)

but since the poissonian would contain the number of observed events in bin i, I have no idea how to compute the NLL and initialize the minimizer

2. use of histograms:

An alternative way to deal with the above problem is to construct my pdf like this:

pdf = RooGenericPdf("pdf"," acceptance(E_i) * int_e_low,i^e_high,i { flux(E) dE ) * time * (1+a*sigma_syst))

with acceptance, sigma_syst & flux being RooFormulaVar objects, obtained from fitting the histograms of acceptance and syst. uncertainty; however, I have no idea how to construct the integrated flux, which would still have to obey the bin boundaries. This method has the advantage, that I can then create the NLL, but It’s absolutely unclear to me how to introduce the integral of the flux over the bin noting that the parameters of interest are implicitly inside the flux formula.

I’ve scoured the web and looked through numerous RooFit tutorials - but I can’t find any example problem to be even close to the one I sketched above.

Thank you so much for taking your time reading this and hopefully you’ll have some ideas.

[moneta]

Hi,

Have you looked at the HistFactory ? It helps in making model based on histograms and adding systematic errors.
See https://cdsweb.cern.ch/record/1456844

Best Regards

Lorenzo

[stephan.zimmer]

Ciao Lorenzo,
indeed I had a look at HistFactory, but I have to admit that I have yet to see how this is applicable to my problem. I’ve looked at the paper with great interest only to see that this tool is great for HEP analysis, but not so much for astro particle physics where we e.g. don’t have luminosities in our calculations. At least that’s what I thought when I read your paper the first time. Now, after re-reading, it’s certainly worth a try to see if I can use it to come up with a sensible result.
Thanks!
-Stephan