Skellam distribution maximum likelihood fitter?

Roofit and RooStats
1

Hello,

I have two invariant mass distributions: signal and (event mixed) background, both of which have the standard Poisson distribution errors. The difference of the two would have Skellam distribution errors and it is this difference that I need to fit to extract the number of resonance decays in my data.

Does ROOT have implementation of Skellam distribution and/or a “Fitter” class using it?

Thank you,
Siarhei.

2

Hi,

We don’t have this implementation in ROOT, however you can implement by using the formula in Wikipedia:

double skellam(double k, double mu1, double mu2) { 
   return  std::exp(-(mu1+mu2))*std::pow(mu1/mu2,0.5*k) * ROOT::Math::cyl_bessel_i(std::abs(k), 2.0 * std::sqrt(mu1*mu2);
}

In your case, I would not subtract the histogram and instead do a global fit to the S+B and B histograms using the Poisson statistics.

Lorenzo

3

Hello Lorenzo,

Could you please clarify what you meant by that? Are you saying there is a way to fit two histograms simultaneously effectively fitting their difference?

Thank you,
Siarhei.

4

Hi,

yes, if you have an histogram with S+B and an histogram with only B, it is much more effective to fit simultaneously S+B with a model for S + a model for B and B with only the B model, than to fit the difference with only the B model. In the first case you have more information.
If the B histogram contains much more statistics that the S+B, you can also fit first the B model in the B histogram and then using the fit result for the B model to fit for only the S model in the S+B histogram.

Best Regards

Lorenzo

5

Hi Lorenzo,

I am not sure what you mean.
[list=]
[] How do I fit two histograms simultaneously?
[] I do not have any model for B - my background is generated using event mixing.
[*] Are you suggesting that I fit my background first and then fix all the parameters of the background fit function (bFitFun) except overall scaling, and then use the function to fit my S+B (for example, BreitWigner + [0]*bFitFun) ?
[/list]

Also I want to at least be able to look at my (S+B)-B, however when I look at the code of void TH1::Add(TF1 *f1, Double_t c1, Option_t *option) , I do not see any recalculation of the histogram errors, taking into account errors on the “f1” parameters. Is it somehow done implicitly?

Thank you,
Siarhei.

6

Hi,

You have to build a global chi2 or likelihood function for the two histograms and then minimize it. A tutorial might exist for RooFit. If you want I could provide you an example using only ROOT.

if you don’t have a parametrized model for the background, you could use the B histogram directly.
I.e. the bFitFun is your B histogram except an overall scaling.
If you think you have a reasanble model for the background, then it is easier to fit it first and then use it for fitting the S+B histogram

[quote] I do not see any recalculation of the histogram errors, taking into account errors on the “f1” parameters. Is it somehow done implicitly?
[/quote]

You are right, in TH1::Add(TF1 *f1, Double_t c1, Option_t *option) the errors on the function are not considered. It is maybe something we could add in the future. Now, you will have to adjust the bin errors yourself using the function error from the fit. The function error from the fit you can get from FitResult::GetConfidenceInterval
See root.cern.ch/root/htmldoc/ROOT__ … eIntervals

Lorenzo