I’m running the CLs for a counting experiment with almost zero background (0.03), 0 observed events, and external constraints on both signal and background. The systematics on signal and background are very high (~40%).
The observed and expected limits seem reasonable, while the 1- and 2-sigma bands have (apparently) problems for s > 1. Increasing the printout at level 3, I don’t see any evident issue in the minimization, but maybe I’m missing something here. Do you have any suggestions on how to help the fit converging? I have also tried to play a little bit with the minimisation parameter without a big success. I’m using ROOT 6.18, but I have also tried an older version (5.34/30) which seems to have fewer troubles in two points of the scan, results are attached.
I attach the macro, the CLs scan plots, and the related profile likelihood plots.
Maybe @moneta has an idea. What I see from the logs is that Minuit has problems since with b=0.03 and data=0, the parameter b hits its lower limit of zero. That’s of course expected, but Minuit might have trouble getting proper errors with a parameter at limit.
Do unconditional fit
**********
** 1 **SET PRINT 0
**********
**********
** 2 **SET NOGRAD
**********
PARAMETER DEFINITIONS:
NO. NAME VALUE STEP SIZE LIMITS
1 b 0.00000e+00 4.95321e-01 0.00000e+00 1.00000e+01
MINUIT WARNING IN PARAM DEF
============== STARTING VALUE IS AT LIMIT.
MINUIT WARNING IN PARAMETR
============== VARIABLE1 IS AT ITS LOWER ALLOWED LIMIT.
MINUIT WARNING IN PARAMETR
============== VARIABLE1 BROUGHT BACK INSIDE LIMITS.
2 bSyst 9.92611e-01 4.03849e-01 0.00000e+00 3.00000e+00
3 s 2.06897e-01 1.03448e-01 0.00000e+00 1.00000e+01
4 sSyst 9.92611e-01 4.03849e-01 0.00000e+00 3.00000e+00
**********
** 3 **SET ERR 0.5
**********
**********
** 4 **SET PRINT 0
**********
**********
** 5 **SET STR 1
**********
**********
** 6 **MINIMIZE 2000 1
**********
MINUIT WARNING IN MINImize
============== VARIABLE1 IS AT ITS LOWER ALLOWED LIMIT.
MIGRAD MINIMIZATION HAS CONVERGED.
MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.
FCN=7.50963e-06 FROM MIGRAD STATUS=CONVERGED 59 CALLS 60 TOTAL
EDM=1.50153e-05 STRATEGY= 1 ERROR MATRIX ACCURATE
EXT PARAMETER STEP FIRST
NO. NAME VALUE ERROR SIZE DERIVATIVE
1 b 3.26189e-10 4.46049e-01 2.07458e-04** at limit **
2 bSyst 1.10413e+00 3.94890e-01 1.34839e-04 8.36164e-03
3 s 3.37334e-06 4.16217e-01 2.02701e-04 -6.87957e-03
4 sSyst 1.18449e+00 3.95044e-01 1.33925e-04 4.76113e-03
ERR DEF= 0.5
Do conditional fit
**********
** 1 **SET PRINT 0
**********
**********
** 2 **SET NOGRAD
**********
PARAMETER DEFINITIONS:
NO. NAME VALUE STEP SIZE LIMITS
1 b 3.26189e-10 4.46049e-01 0.00000e+00 1.00000e+01
MINUIT WARNING IN PARAMETR
============== VARIABLE1 IS AT ITS LOWER ALLOWED LIMIT.
MINUIT WARNING IN PARAMETR
============== VARIABLE1 BROUGHT BACK INSIDE LIMITS.
2 bSyst 1.10413e+00 3.94890e-01 0.00000e+00 3.00000e+00
3 sSyst 1.18449e+00 3.95044e-01 0.00000e+00 3.00000e+00
**********
** 3 **SET ERR 0.5
**********
**********
** 4 **SET PRINT 0
**********
**********
** 5 **SET STR 1
**********
**********
** 6 **MINIMIZE 1500 1
**********
MINUIT WARNING IN MINImize
============== VARIABLE1 IS AT ITS LOWER ALLOWED LIMIT.
MIGRAD MINIMIZATION HAS CONVERGED.
MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.
FCN=0.241535 FROM MIGRAD STATUS=CONVERGED 35 CALLS 36 TOTAL
EDM=3.90777e-07 STRATEGY= 1 ERROR MATRIX ACCURATE
EXT PARAMETER STEP FIRST
NO. NAME VALUE ERROR SIZE DERIVATIVE
1 b 0.00000e+00 4.46418e-01 2.53079e-04** at limit **
2 bSyst 1.10321e+00 3.94922e-01 1.64383e-04 -1.17310e-05
3 sSyst 1.15112e+00 3.94999e-01 1.61800e-04 2.27993e-03
ERR DEF= 0.5
EvaluateProfileLikelihood - mu hat = 3.37334e-06, uncond ML = 7.50963e-06, cond ML = 0.241535 pll = 0.241528 time (create/fit1/2) 0 , 0 , 0
POIs:
1) 0x7ff8bd7ad680 RooRealVar:: s = 0.206897 L(0 - 10) "s"
reusing NLL 0x7ff8bfa3a400 new data = 0x7ff8c80ce440
Data set used is: ( nobs = 0)
The problem you observed is probably due to too many fits failing for the background toys when computing the likelihood ratio at large s value. The Blue curve in the test statistic distribution plots.
This causes to be difficult to estimate the +2 sigma.
Strange there is this difference with 5.30, we would need to investigate why.
One think one can try is also to remove the lower bound to zero and put a lower value .The fitted parameters s and b can be negative, it is only their sum that must be >= 0.
Thanks a lot for your prompt replies!
I have a naive question, can I fix the minimum of “nexp” in my expression “sum:nexp(sig,bkg)”? Since “nexp” is a sum of variables, I was naively expecting to be able to set a range on its possible values, but I haven’t found any working solution for it. My second option was to try to constraint “nexp” using a Heaviside function. Does this make sense?
You can set limits for the inputs of the expression a+b, but not for the result (think about this: If the result is not allowed, who is at fault? Should you change a or b?).
Constraining is certainly possible, but a heavyside is brutal. The fitter might not understand what’s wrong if it steps over the boundary.
So, instead, use something like a quadratic penalty term. Penalty = (a+b) < 0 ? (a+b)^2 : 0
If it doesn’t kick in soon enough, shift it to kick in at slightly positive values of a+b. If it doesn’t kick in strong enough, scale it.
thanks a lot for the suggestion!
I have tried to add a penalty (leaving the signal and background floating to negative values). First, I have tried to add it via RooFormulaVar, but this is not possible as far as I understood (right?). I mean, it should be added to the NLL, but, if I want to do it in my code, I have to create my own test statistic from ProfileLikelihoodTestStat so that it uses the new likelihood shifted by a factor (s+b)^2 when (s+b)<0, correct? For example, after line 92 in ProfileLikelihoodTestStat.cxx I can add
I have also tried to add an additional pdf to constrain (s+b) to be positive (that is basically a Heaviside), but as you anticipated, the fit doesn’t like it too much. I was thinking to add an additional pdf (as penalty) to my model (using RooAddPdf) to scale the likelihood when (s+b)<0, but it doesn’t seem a good idea.
Do you know if there is some example? I haven’t found anything around to add a penalty with ProfileLikelihoodTestStat.
Looking better at your original model, you will never be able to compute a 2 sigma bands. This is due to your weak systematic constraint. Basically you have a probability too large that the constraint on your signal parameter(sSyst) is around zero. Now if you write the likelihood you can see that if sSyst is approximatly zero your likelihood is not anymore dependent on s, so you don’t have any way to make an inference on the s parameter. In conclusion, with that model the +2sigma limit is +infinity.
So your original scan plots after all looks correct