I am fitting sidebands of the dimuon mass with polynomial background function in the mass intervals of 12-25 and 25-70 GeV. Now I want to know a number of expected background events in the signal region 25-30 GeV and its error.
I take integral of the fitted function in the mass range of 25-30 GeV and it works.
My question is: how I can get an error of the integral ?
thanks ! It work, but gave me a strange result. I do
// define model
RooPolynomial bkg("bkg","background",dimuonmass1sideband,RooArgList(a0,a1,a2));
RooRealVar norm("norm","bkg normalization",1100.,0.0,3000.0);
RooAddPdf model("model","model",RooArgList(bkg),RooArgList(norm));
// fit
RooFitResult* r = model.fitTo(data,Range("R1,R2"),Save());
// Integral of bkg function in the signal region
dimuonmass1sideband.setRange("signalrange",25,30);
RooAbsReal* intPolynX = model.createIntegral(dimuonmass1sideband,NormSet(dimuonmass1sideband),Range("signalrange"));
Double_t integral_undersignal = intPolynX->getVal();
Double_t err_integral_undersignal = intPolynX->getPropagatedError(*r);
cout <<" Integral under the signal region = " << integral_undersignal <<" +/-" << err_integral_undersignal << endl;
and print gives:
Integral under the signal region = 0.133652 +/-1.88998
So, I do not understand why error is larger than integral. Is it in % ?
Thanks ! I do not know yet how ROOT Forum works, but I put my replay to your answer through Forum.
I got a strange thing: integral error is larger than value. Please, check my replay for details.
Thanks for sharing the code, you are doing things correctly. Just one thing that you need to change (my fault for not mentioning it before): if you want the get the expected number of background events, you have to multiply the integral with the yield, and then do the error propagation. You have to do this, because the norm parameter and the pdf shape parameters are generally correlated. If you are not propagating the uncertainty up to the final yield, your result will be wrong.
std::unique_ptr<RooAbsReal> intPolynX{model.createIntegral(dimuonmass1sideband,NormSet(dimuonmass1sideband),Range("signalrange"))};
RooProduct yieldPolynX{"yieldPolynX", "yieldPolynX", {norm, *intPolynX}};
double integralUndersignal = yieldPolynX.getVal();
double errIntegralUndersignal = yieldPolynX.getPropagatedError(*r);
std::cout << " Integral under the signal region = " << integralUndersignal
<< " +/-" << errIntegralUndersignal << std::endl;
Concerning your surprisingly large integral value: apparently the uncertainty on your PDF parameters is so large that your pdf can become negative within the parameter uncertainties. Have you inspected your fit result? My guess is that for some of the polynomial coefficients, there will be a huge uncertainty. And this is probably because your are trying to fit a 3rd order polynomial, which is not very stable because the parameter are very correlated (remember that the coefficients you pass to the RooPolynomial constructor already start with the first order, because the zeroth coefficient is fixed by the normalization condition).
To stabilize the fit, I would use the more robust parametrization of the Chebychev polynomials implemented in RooChebychev. You can also consider the Bernstein polynomials implemented in RooBernstein. They are not as orthogonal as the Chebychev polynomials, but sometimes better for fitting because they are positive definite and the PDF can’t become negative.
If you stabilized your fit model, the problem with the large uncertainties should also be gone.
I hope I’m on the right track here and my answer helps!
Indeed, your problem is not the one I thought, and I think I was wrong with my initial suggestion of doing the error propagation.
For a mathematical reason that I can’t exactly pin-point my finger on just yet, the error propagation is not valid in this case. An safer alternative that just came to mind now that I see your script without error propagation would be to reformulate your model such that you don’t need it, and one of the parameters is already the background yield in the signal region.
In the RooAbsPdf::fitTo() documentation, you see explained the SumCoefRange() option to define for which range the yield coefficients are defined. You can pass SumCoefRange("signalrange") to your fitTo call, and the norm parameter will be directly the normalization you are looking for. Then you can drop all the other integration stuff.
I hope it’s fine for you now, otherwise feel free to follow up again!