Bernstein Polynomial + Gaussian (or Crystal Ball) fit to a histogram

Newbie
1

OS: Windows 10
Root Version : ROOT 6.15/01

Hello,

I have a histogram that I am trying to fit with a combination of Bernstein Pol (background) and a Gaussian (or a Crystal Ball) for the signal region. Previously, I had tried doing this with expo and a gaussian (crystal ball) but the fit was poor.

Capture
Capture.JPG1086×512 35.5 KB
[This fit was with expo + gaussian]

RooRealVar x("x","mass",18,500);
RooDataHist dh("dh","e-e+ peak histo",x,Import(*h1)) ;
RooRealVar mass("mass","Central value of Gaussian",90,80,100);
RooRealVar sigma("sigma","Width of Gaussian",20,0,100);
RooGaussian signal("gaus", "The signal distribution", x, mass, sigma); 
RooBernstein bg_bern("bg_bern","background",x, <! here it says it needs a coefficient list>);

How do I go about this?

2

Hi,

you can choose the order of the polynomial by giving the number of coefficients. See e.g. here:
https://root.cern.ch/doc/master/classRooBernstein.html

That would look something like this:

RooRealVar a1("a1","a1",....);
RooRealVar a2("a2","a2",....);
RooRealVar a3("a3","a3",....);
RooBernstein bg_bern("bg_bern","background",x, RooArgList(a1, a2, a3));
3

Hello @StephanH,

I was studying the Bernstein Pol. (here), and saw that an order N polynomial is defined as
B_i,N(t) = (N,i)* t^i (1-t)^(N-i)
Could you tell me what a1,a2,a3 are in the given code ( I am new to fitting with polynomials).
Thanks,
Debajyoti

4

These are the coefficients of the single Bernstein polynomials.

For each order N you have N+1 polynomials that you can sum to form a basis of the space of polynomials of order N or less. These coefficients represent an element of the space of polynomials.
So it’s
a0 * B_{0,N} + a1 * B_{1,N} …

5

Given my histogram, I tried a BP(2)+Gaussian fit (it didn’t work too well). Is there a way by which I can figure out what order of BP I actually need?

Screenshot%20(41)
Screenshot (41).png1366×768 51.5 KB 

RooRealVar x("x","mass",18,500);
         x.setRange("Range1",20.,52.3) ;
         x.setRange("Range2",106.9,500.) ;
         RooRealVar a1("a1","a1",0,10);
         RooRealVar a2("a2","a2",0,10);
         RooRealVar a3("a3","a3",0,10);
         RooDataHist dh("dh","e-e+ peak histo",x,Import(*h1)) ;
         
         

         RooRealVar mass("mass","Central value of Gaussian",90,80,100);
         RooRealVar sigma("sigma","Width of Gaussian",20,0,100);
         RooGaussian signal("gaus", "The signal distribution", x, mass, sigma); 

         RooRealVar b("b", "Number of background events", 0, 1000);
         RooRealVar s("s", "Number of signal events", 0, 1000);
         RooBernstein bg_bern("bg_bern","background",x,RooArgList(a1,a2,a3));

         RooAddPdf fullModel("fullModel", "Signal + bg_bern", RooArgList(signal, bg_bern), RooArgList(s, b));
         RooFitResult* r = fullModel.fitTo(dh) ;
         RooPlot* frame = x.frame();
         dh.plotOn(frame);
         fullModel.plotOn(frame);
         // expo.plotOn(frame);
         
         frame -> Draw();
6

The model looks good. The degree of the polynomials can be hard to determine … Trial and error is the only way if you absolutely don’t know the physics process that’s causing the background distribution.

However
Don’t you think that the number of background events is way too low? In the first 10 bins you already have > 1000 data events. You shouldn’t limit b to 1000.
Look at the fit output. It should say somewhere that b is at limit.

7

It was at the limit. I changed background # to 4000 and tried again. Should I alter the signal events too accordingly?

Also, in order to “guess” the correct order, do I need to keep adding a roorealvar everytime I need to increase the order? Or is there a smarter way of going about it?

8

Try for yourself with the number of events. You can also not give limits at all.

For the order:
It only depends on the length of the RooArgList.

9

so can I just repeat elements in the rooarglist to increase the order?
For example, with roorealvar a1,a2,a3 already defined can I do rooarglist(a1,a2,a3,a1,a2,a3…)?

10

That’s a bad idea, unless you really want to have the polynomials 0 and 3 to have the same coefficients. Instead, you should add more RooRealVars to hold the coefficients.

11

Hi,

I tried increasing the order of the polynomial (9 coefficients) and ended up this fit:

Screenshot%20(42)
Screenshot (42).png1366×768 58.3 KB

What I want to do is limit the fit range to x ~170. I tried

x.setRange("Range1",10.,170.) ;
RooAddPdf fullModel("fullModel", "Signal + bg_bern", RooArgList(signal, bg_bern), RooArgList(s, b));
RooFitResult* r = fullModel.fitTo(dh,Range("Range1"),Save(kTRUE)) ;

But that doesn’t work

Screenshot%20(43)
Screenshot (43).png1366×768 48.1 KB

Additionally, I cannot seem to add arbitrarily many coefficients to the RooArgList. The moment I add the tenth coefficient it spits out an error :

RooRealVar a1("a1","a1",-50,400);
         RooRealVar a2("a2","a2",-50,400);
         RooRealVar a3("a3","a3",-50,400);
         RooRealVar a4("a4","a4",-500,1000);
         RooRealVar a5("a5","a5",-50,1000);
         RooRealVar a6("a6","a6",-50,100);
         RooRealVar a7("a7","a7",-50,100);
         RooRealVar a8("a8","a8",-50,100);
         RooRealVar a9("a9","a9",-50,100);
         RooRealVar aa("aa","aa",-50,100);
         RooDataHist dh("dh","e-e+ peak histo",x,Import(*h1)) ;
         RooRealVar lambda("lambda", "slope", -0.01, -1., 1.);
         RooExponential expo("expo", "expo", x, lambda);
         

         RooRealVar mass("mass","Central value of Gaussian",90.,80.,100.);
         RooRealVar sigma("sigma","Width of Gaussian",20,0,100);
         RooGaussian signal("gaus", "The signal distribution", x, mass, sigma); 

         RooRealVar b("b", "Number of background events", 0, 40000);
         RooRealVar s("s", "Number of signal events", 0, 40000);
         RooBernstein bg_bern("bg_bern","background",x, RooArgList(a1,a2,a3,a4,a5,a6,a7,a8,a9,aa));

In file included from input_line_10:1:
D:\debajyoti_project\zdecay.C:166:57: error: no matching constructor for initialization of 'RooArgList'
         RooBernstein bg_bern("bg_bern","background",x, RooArgList(a1,a2,a3,a4,a5,a6,a7,a8,a9,aa));
                                                        ^          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
D:\buildroot\include\RooArgList.h:57:3: note: candidate constructor not viable: no known conversion from 'RooRealVar' to 'const char *' for 10th argument
  RooArgList(const RooAbsArg& var1, const RooAbsArg& var2

But this works for 9 coefficients

12

You could try Range(10., 170.), see RooAbsPdf::fitTo.
“Range1” already exists as a range in your first script. Does this create a problem?

The RooArgList does not take an arbitrary number of arguments in the constructor, see RooArgList. However, you can always use argList.add(myVariable).

13

Range(10., 170.)also gives the same fit. I tried both of them.

14

Please post the logs that RooFit is writing to the console.

15

Hello, the output was rather large. I am attaching a txt.

console.txt (506.7 KB)

16

If you look at the outputs, you see that you get tons of warnings and ultimately a failed fit. It could be that you have so many parameters in your model that fitting with so many coefficients fails because the fit is overdetermined.

If you make the range larger, or reduce the number of coefficients, does it converge?

Successful fits have status 0, the status codes for failed fits can be found here:
Minuit2::Status.

17

I managed to get a fit with Crystal ball and bernstein polynomial. Turns out the issue was the number of signal events and background events. bernstein_CB.pdf (29.4 KB)

This is fairly close to what I want, but I wonder if I can squeeze it a little further.

18

Ok, the technicalities are taken care of. The interpretation of results I leave to you and your group.