Problems with 1- and 2- sigma bands in the CLs scan (high signal and bkg systematic)

Hello,

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.

Thanks a lot in advance!

Regards,
Cristiano

UpperLimit.C (5.4 KB) CLsScan_ROOT5p30.pdf (20.4 KB) ProfileLikelihoodPlots_ROOT5p30.pdf (117.4 KB) CLsScan_ROOT6p18.pdf (20.4 KB) ProfileLikelihoodPlots_ROOT6p18.pdf (115.9 KB)

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)

Hi,

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.

Lorenzo

Hi @StephanH and @moneta!

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?

Thanks a lot!

Cheers,
Cristiano

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.

Hi @StephanH,

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

// Add penalty                                                                                                                                                                              
RooRealVar *s     = w->var("s");
RooRealVar *sSyst = w->var("sSyst");
RooRealVar *b     = w->var("b");
RooRealVar *bSyst = w->var("bSyst");
RooFormulaVar *penalty = new RooFormulaVar("penalty","(pow(s*sSyst+b*bSyst,2))*((s*sSyst+b*bSyst)<0)",RooArgList(*s,*sSyst,*b,*bSyst));
RooAddition *fNllP = new RooAddition("fNllP","fNllP",RooArgSet(*fNll,penalty));

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.

Thanks a lot!

Cheers,
Cristiano

Yeah, when searching for “penalty”, you wouldn’t find it. The term is “constraint”.

Here is an example:
https://root.cern.ch/doc/master/rf604__constraints_8C.html

Hi,

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

Lorenzo

Hi @StephanH,

of course, RooProdPdf, sorry for the stupid question!

Thanks a lot for the useful discussion!

Cheers,
Cristiano

1 Like

Hi @moneta,

that’s indeed a good point! L’uovo di Colombo…
I agree that with a 40% systematic level the 2 sigma upper band has not constraints.

Thanks a lot for the discussion!

Cheers,
Cristiano

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