Fit a double Crystal Ball with RooFitRanges

Roofit and RooStats
,
1

Hello Experts,
I try to fit some data with a pdf consisting of two combined Crystal Ball functions for the signal and an exponential for the combinatorial background. In addition a data window is blinded. Therefore I use RooFit.Range to specify the range on which the pdf should be fitted.
Now the problem is that the fit only converges if I fit on the full range which is of course not right in the present of this blinded window.
I found out that it works when using a gaussian as a signal like in my previous question: roofitranges with a gaussian.
My Code is here:

import ROOT as R

# Create observables
x = R.RooRealVar("x", "x", 5100, 5700)
werte = 5000
bkg = 500

left = 5350
right = 5400

# Construct signal pdf
mean = R.RooRealVar('mean', '', 5280, 5250, 5300)
sigma = R.RooRealVar('sigma', '', 6, 0, 40)
alpha_1 = R.RooRealVar('alpha_1', '', 1, 0, 100)
alpha_2 = R.RooRealVar('alpha_2', '', 1, 0, 100)
n_1 = R.RooRealVar('n_1', '', 3, 1, 100)
n_2 = R.RooRealVar('n_2', '', 3, 1, 100)

cbs_1 = R.RooCBShape("CrystallBall_1", "CrystallBall_1", x, mean, sigma, alpha_1, n_1)
cbs_2 = R.RooCBShape("CrystallBall_2", "CrystallBall_2", x, mean, sigma, alpha_2, n_2)
mc_frac = R.RooRealVar('mc_frac', '', 0.5)

sig = R.RooAddPdf('sig', '', R.RooArgList(cbs_1, cbs_2), R.RooArgList(mc_frac))
sig_yield = R.RooRealVar('sig_yield', '', 1794.4, 0, 100000)



exp = R.RooRealVar('exp', '', -2e-3, -4e-3, -0.01e-3)
combinatorial = R.RooExponential('combinatorial', '', x, exp)
bkg_yield = R.RooRealVar('bkg_yield', '', bkg, 0, werte)


pdf_cbWithExp = R.RooAddPdf(
    'pdf_cbWithExp',
    '',
    R.RooArgList(sig, combinatorial),
    R.RooArgList(sig_yield, bkg_yield)
)

x.setRange("low", 5100, left)
x.setRange("high", right, 5700)
x.setRange("full", 5100, 5700)

fit_range = "low,high"

data = R.RooDataSet(pdf_cbWithExp.generate(R.RooArgSet(x), werte))
data_red = data.reduce("x<{}||x>{}".format(left, right))

pdf_cbWithExp.fitTo(data_red, R.RooFit.Range(fit_range), R.RooFit.NumCPU(4), R.RooFit.Save(True), R.RooFit.Extended())
my_plot = x.frame()

data_red.plotOn(my_plot)
pdf_cbWithExp.plotOn(my_plot, R.RooFit.Range(fit_range), R.RooFit.NormRange(fit_range))
pdf_cbWithExp.getParameters(data).writeToFile("results.txt")

can = R.TCanvas()
can.SetLogy()
my_plot.Draw()

can.SaveAs("blind_test_ranges.pdf")

The results are:

alpha_1 =  2.2282 +/- 21.297 L(0 - 100) 
alpha_2 =  1.0000 +/- 22.525 L(0 - 100) 
bkg_yield =  500.00 +/- 1262.2 L(0 - 5000) 
exp = -0.00200000 +/- 0.0016787 L(-0.004 - -1e-05) 
mc_frac =  0.50000 C L(-INF - +INF) 
mean =  5280.0 +/- 20.612 L(5250 - 5300) 
n_1 =  3.0000 +/- 21.836 L(1 - 100) 
n_2 =  3.0000 +/- 21.836 L(1 - 100) 
sig_yield =  1794.4 +/- 22160 L(0 - 100000) 
sigma =  6.0000 +/- 12.019 L(0 - 40)

and the terminal output is given here: terminal_output.txt (1.7 MB).
The plot looks like this: blind_test_ranges.pdf (23.9 KB).

My root version is 6.17/01 built for linuxx8664gcc on Feb 01 2019.

Thanks in advance,
Jasper

2 
     p.d.f value is Not-a-Number (-nan), forcing value to zero @ !refCoefNorm=(x = 5676.03), !pdfs=(CrystallBall_1 = 0/0,CrystallBall_2 = 0/0), !coefficients=(mc_frac = 0.5)

This looks like both crystal ball distributions are completely inside the blinded window. When it’s being evaluated, the likelihood is zero. Worse, though, is that the integrals over the crystal ball distributions are zero on the active ranges. In the end, RooFit divides 0/0 (function value / integral) to compute the likelihood, which yields NaN.
What would be needed to stabilise this is a Crystal ball that’s wide enough such that the integral yields a non-zero probability or the crystal balls need to be deactivated. Did you try to make the crystal balls very wide such that the integrals leak out of the blinded range?

3

Yeah you’re right. If I put the sigma on 30, the fit converges and delivers this plot blind_test_ranges.pdf (23.7 KB) and this result:

alpha_1 =  1.0460 +/- 0.25321 L(0 - 100) 
alpha_2 =  1.0460 +/- 0.25321 L(0 - 100) 
bkg_yield =  1103.7 +/- 127.24 L(0 - 5000) 
exp = -0.00212929 +/- 0.00052837 L(-0.004 - -1e-05) 
mc_frac =  0.50000 C L(-INF - +INF) 
mean =  5279.9 +/- 0.85934 L(5250 - 5300) 
n_1 =  2.6226 +/- 62.992 L(1 - 100) 
n_2 =  2.6226 +/- 62.992 L(1 - 100) 
sig_yield =  3913.2 +/- 132.60 L(0 - 100000) 
sigma =  31.715 +/- 0.87866 L(0 - 40)

But how could it be solved for slimmer peak?

4

@StephanH Do you know how the fit could be brought to conversion? It does not converges if the sigma is below 17 which is a rather strange behavior in my opinion. I think there is a problem with the RooCBShape class I use here. Because it functions well when trying for example RooGaussian or RooBreitWigner.

5

Hi @JLammering,

I can make a simple fix: Instead of an integral which is exactly zero, I can make RooCBShape return an extremely small number.
RooFit should stop complaining, and the fit should proceed.

I just merged a pull request for this change:

We now have to wait for a new ROOT to be build over night, and then you can test one of the nightlies on lxplus.
Check https://root.cern/nightlies tomorrow to obtain a nightly build.

6

Hi @StephanH,

I finally tested the fix. It seems like the problem is still present.
The current terminal output is here: terminal_output_2.txt (326.8 KB) .

7

Hello @JLammering,

Apparently, I wasn’t generous enough, and the integral was returning a number that’s too small. With 1.E-300, which is still pretty small, it works:
blind_test_ranges.pdf (23.9 KB)

I will commit this soon:

8

I tested the nightly today and it works! Thank you.

9

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.