How to configure constraints in ModelConfig

[aaronsw]

Dear Experts,

I’m sorry if this has been answered somewhere, but I haven’t been able to find it. I’m trying to implement a gaussian (luminosity) constraint in a confidence interval (frequentistCalculator). I’m puzzled as to how frequentistCalculator knows how much to vary the luminosity when it throws toys? Here’s a summary what I’m doing:

w = RooWorkspace("w","w")

# Set up constant+Gaussian S+B model
w.factory("x[0,10]")
w.factory('EXPR::bkgModel("1",x)')
w.factory('Gaussian::sigModel(x,mean[1],sigma[1])')
w.factory('SUM::model(mu[0,10]*sigModel,bkgModel)')

# Set up luminosity constraint
w.factory('Gaussian::lumiConstr(lumiObs[0,2],lumiNp[0,2],0.1)')

# Make model with constraint
w.factory('PROD::modelConstr(model,lumiConstr)')

# But how to set up the modelConfig so that it has all this information?
w.defineSet("poi","mu");
w.defineSet("nuis","lumiNp");
w.defineSet("globObs","lumiObs");

mc = ModelConfig("hypo",w)
mc.SetPdf(w.pdf("modelConstr"))
mc.SetObservables(RooArgSet(w.var("x")))
mc.SetParametersOfInterest(w.set("poi"))
mc.SetNuisanceParameters(w.set("nuis")) # do I need to do this, or is SetGlobalObservables enough?
mc.SetGlobalObservables(w.set("globObs"))
mc.SetSnapshot(w.set("poi"))

How do I provide the lumiConstr gaussian constraint to modelConfig? Is there the equivalent of ExternalConstraints that works with pdf.fitTo()? Thanks very much for any help with this.

(I’m following a procedure from these:
https://twiki.cern.ch/twiki/bin/view/RooStats/RooStatsTutorialsAugust2012#Create_Poisson_Counting_model
https://root.cern.ch/root/html/tutorials/roofit/rf604_constraints.C.html)

[bellenot]

I’m sure @StephanH can help

[StephanH]

Hi @aaronsw,

Some answers:

• You wonder if you have to set the nuisance parameters. The answer is: maybe. The ModelConfig has some auto guessing of parameters/observables. The PDF understands what are parameters and what are observables by looking at the data. Everything that has an entry in the dataset is an observable, the rest are parameters. All parameters that are found to be constant are removed from this set, and this forms the automatic set of nuisance parameters.
  In your case the automatic guessing would probably miss it, because the lumiObs is not in the dataset.
• Generating data points for the constraints: The constraint term has a global observable, the target so to say, lumiObs. For every toy run, one can generate a new value for lumiObs, where lumiObs comes from a Gaussian distribution with mean lumiNp and sigma 0.1. The ToyMCSampler does exactly that with the GlobalObservables. It generates a single value, and goes on to generate multiple values for the main observable(s) afterwards.
  To summarise: As long as you use the full constrained model as PDF and add the constraint observable to the global observables, it should work.

[aaronsw]

Hi @StephanH,

Thanks very much for the information! I have one clarification question. When I make PROD::modelConstr(model,lumiConstr), how does RooFit know that lumiConstr is the Gaussian to be used by ToyMCSampler? For example, what if there were two Gaussians in the product?

Thanks,
Aaron

[StephanH]

With N Gaussians, you would create N global observables. Your likelihood should look like:

L = L(x | a, b, c, ...) * Gauss(obs_a | a, 0.1) * Gauss(obs_b | b, 0.2) * ...

The ToyMCSampler will go through all global observables and generate a data point. Since obs_a only shows up in the Gaussian distribution, its generated distribution will also only be Gaussian.

[aaronsw]

Okay, I see. Thanks very much for the answers, they are very helpful.