Custom constraint function in HistFactory

[JanEric]

Hi,

we are doing a fit using a model constructed to HistFactory.

We are using sample.AddHistoSys() for our coherent shape uncertainties.

The gaussian constraint this uses work perfect for almost all of our uncertainties.

However we have one modelling uncertainty which is really an unconstrained template choice between 2 options where we feel it would be best to have a uniform distribution between 0 and 1 and thats it.

This could be parametrized with a double fermi distribution.

Now my question is how i could best/easiest add this systematic to our measurement given this constraint.

Cheers
Jan-Eric

[mczurylo]

Hi @JanEric,

thank you for your question, maybe @Jonas could help here?

Cheers,
Marta

[jonas]

Hi @JanEric,

thanks for the question and sorry for the late reply!

Before I can give a clear recommendation here, I need to understand a little bit better what you want to achieve mathematically. For sure, what you are trying to achieve here can’t be done with vanilla HistFactory, we will need to dive a bit deeper into RooFit.

So for each histogram systematic, HistFactory creates an associated nuisance parameter that is constrained by a Gaussian. So what you want is a constraint for two nuisance parameters \alpha_1 and \alpha_2 for two different variations that is constraining them like this:

p(\alpha_1,\alpha_2)= 0 if \alpha_1 != 0 and \alpha_2 != 0 , else uniform.

But to avoid the discrete step such that there are non-vanishing gradients to guide the minimization you want to “smooth out” the transition with a sigmoid/Fermi function? Can you confirm this and let me know the mathematical formula for this constraint?

It’s a very interesting idea, and once I understand it better I can see how to best implement this in RooFit and give you a recommendation.

Cheers,
Jonas

[JanEric]

Hi Jonas,

so if i am not mistaken normally for each systematic variation we have a Nominal templates (usually shared between all systematics) and a template for an up and a down variation. The down, nominal and up variation correspond to values of -1, 0 and 1 of the assigned nuisance parameter.

Then for other values the templates are inter/extrapolated.

Now one systematic effect we are considering are higher order corrections. For these we effectively have the nominal case without them and then the variation case where we fully apply them (with a complicated procedure due to which we do not trust them 100% which is why we are not using them directly as nominal).

Currently i have set this up via HistFactory by setting the down variation to the same as the nominal and then the up variation corresponds to the NLO template.

However with the default constraint the fit punishes moving the NP to the NLO template which we would like to avoid. We also do not want to extrapolate to negative NP values or beyond the NLO template as that does not make much sense physically.

So our idea was to swap the gaussian constraint with a double double fermi distribution which is effectively flat between 0 (nominal) and 1 (NLO) but falls to 0 very quickly for everything else.

Cheers
Jan-Eric

[jonas]

Hi, thank you very much for the clarification!

In that case, why not restrict the range of the nuisance parameter to be between zero and one? Then you don’t need the double Fermi distribution, which I find not super elegant because it requires an arbitrary decay parameter. You can then use just a uniform constraint, which is also supported by HistFactory.

So if you have a HistFactory workspace with a HistoSys called “SignalShape”, you could do something like this (the workspace is called ws):

   // The name of the constrained parameter
   const std::string paramName = "alpha_SignalShape";

   // As we want to restrict the variations between the nominal and the "up"
   // variation, we restrict the parameter to that range
   ws->var(paramName)->setRange(0.0, 1.0);

And before creating your workspace, when you have your RooStats::HistFactory::Measurement object, you set the constraint for the systematic to uniform like this:

   meas.AddUniformSyst("SignalShape");

Here is a full ROOT macro that demonstrates this procedure.
custom_constraint_example.C (2.6 KB)

If you print out the workspace, the constraint will appear like this:

RooGaussian::alpha_SignalShapeConstraint[ x=alpha_SignalShape mean=nom_alpha_SignalShape sigma=100 ] = 1

Don’t be surprised it’s still a Gaussian. RooFit encodes uniform constraints with Gaussians that have an “infinite” standard deviation.

I hope this helps you to set up your fit, don’t hesitate to follow up with more questions! In particular, if the uniform constraint is not good for you, I can also explain you how to use arbitrary constraints, although that is hacky and I rather not share the code to do that if not absolutely necessary

Cheers,
Jonas