Global Observable and Non-Gaussian constraint function

[xxiang4]

Hi,

My experiment has code based on RooFit framework to handle profile likelihood ratio to generate limit. I need to build a likelihood model with nuisance parameters to pass in. The snippet below is a simplified model, which contains one signal PDF, one background PDF, and one constraint function.

{
  RooWorkspace* w = new RooWorkspace("w");

  // observable                                                                                                      
  RooRealVar x("x","x",-10,10);
  w->import(x);
  w->defineSet("obs", "x");

  // Parameter Of Interest                                                                                           
  RooRealVar mu_sig("mu_sig","Signal expectation", -100., 1000.);
  w->import(mu_sig);
  w->defineSet("poi","mu_sig");

  //Nuisance Parameters                                                                                              
  RooRealVar mu_bkg("mu_bkg","background expectation", 0, 1000);
  w->import(mu_bkg);
  w->defineSet("nuis", "mu_bkg");

  // signal + background model                                                                                       
  w->factory("Gaussian::gauss_sig(x, m_sig[0], s_sig[1])");
  w->factory("Gaussian::gauss_bkg(x, m_bkg[0], s_bkg[2])");
  TString EventModel_str = "SUM::EventModel(mu_sig*gauss_sig, mu_bkg*gauss_bkg)";
  w->factory(EventModel_str);

  //constraint functions to limit the range of nuisance                                                              
  w->factory("Gaussian::Constraint_bkg(mu_bkg, m_mu_bkg[10], s_mu_bkg[2])");
  //w->factory("CBShape::Constraint_bkg(mu_bkg, m_mu_bkg[10], s_mu_bkg[2], alpha[2], n[5]))");                       

  // Model with Constraint                                                                                           
  TString TotalEventModel_str("PROD::TotalEventModel(EventModel,Constraint_bkg)");
  w->factory(TotalEventModel_str);

  //global observable                                                                                                
  w->defineSet("gobs", "m_mu_bkg");

 // Create Model Config                                                                                             
  //----------------------                                                                                           

  // SIGNAL plus BACKGROUND model                                                                                    
  ModelConfig* SBmodel = new ModelConfig("SBmodel", w);
  SBmodel->SetName("SBmodel");
  SBmodel->SetPdf(*w->pdf("TotalEventModel"));
  SBmodel->SetObservables(*w->set("obs"));
  SBmodel->SetParametersOfInterest(*w->set("poi"));
  SBmodel->SetNuisanceParameters(*w->set("nuis"));
  SBmodel->SetGlobalObservables(*w->set("gobs"));
  w->import(*SBmodel);

  // BACKGROUND-only model                                                                                           
  ModelConfig* Bmodel = (ModelConfig*) SBmodel->Clone();
  Bmodel->SetName("Bmodel");
  RooRealVar* firstPOI = (RooRealVar*) SBmodel->GetParametersOfInterest()->first();
  firstPOI->setVal(0.);
  Bmodel->SetSnapshot(*firstPOI);
  w->import(*Bmodel);

}

Two questions:

1. When to define global observable? The example above set Gaussian mean m_mu_bkg as a global observable, but not the width s_mu_bkg. I’ve also seen other examples online that only define the mean in the constraint function as global observable. Why?

2. Can I use other built-in function (ex. Crystal Ball) for the constrain function instead of Gaussian? If yes, do I need to define all four parameters of Crystal Ball as the global observable? Is it sufficient to set the constraint parameters to constants (ex. use setConstant(kTREU)) instead of defining them as global observable?

Many Thanks,

Xin

[StephanH]

Hi @xxiang4,

@xxiang4

When to define global observable? The example above set Gaussian mean m_mu_bkg as a global observable, but not the width s_mu_bkg . I’ve also seen other examples online that only define the mean in the constraint function as global observable. Why?

When you have a likelihood model that has a constraint term, then you need a “global observable”. A likelihood model that looks like this:

LH(data | parameters) = FitLikelihood(observable | paramA, paramB, paramLumi) * Gauss(lumi | paramLumi, sigmaLumi)

Has lumi as a “global observable”, that’s the observable that someone else measured for you. Think e.g. that the luminosity parameter should be close to 10 / fb, so you introduce the “observable” lumi, and you add a luminosity parameter, which you want to constrain to be close to 10. So you set the Gaussian to:

Gauss(10 | paramLumi, 1)

to get a 10% uncertainty for the luminosity parameter, centred around the external measurement of 10 that you got from the lumi group of your experiment.
Some people interchange the role of x and mu in the gaussian distribution, but that doesn’t matter because it’s symmetric. Sigma on the other hand is a parameter of the constraint term, it cannot be an observable. An observable is something that you can measure.

@xxiang4

Can I use other built-in function (ex. Crystal Ball) for the constrain function instead of Gaussian? If yes, do I need to define all four parameters of Crystal Ball as the global observable? Is it sufficient to set the constraint parameters to constants (ex. use setConstant(kTREU) ) instead of defining them as global observable?

Of course you can use a Crystall ball. Just do

LH_constr = LH_fit(x | param) * CrystalBall(globObs | param, ...)

Note that only the thing that has the role of x in the crystal ball is an observable. The parameters are parameters. And yes, most of the time you will set all parameters of the constraint term constant, because they should not be fitted. Only the parameter that is shared between the fit likelihood and the constraint likelihood will be floating in the fit.

[xxiang4]

Hi @StephanH

Thanks for the reply. You mentioned “Some people interchange the role of x and mu in the gaussian distribution”. Was it a matter of interpretation? For instance, the someone else may claim the luminosity parameter is an expected value from a simulation, and lumi is the one that actually has impact on an experimental outcome. In this case, we would have lumi in the fit likelihood and paramLumi becomes global observable:

FitLikelihood(observable | paramA, paramB, lumi)

In the case of an asymmetric constraint function, it does matter which interpretation it is, right?

Another question. The shared parameters between fit likelihood and constraint likelihood are the nuisance parameters. In your example, paramLumi is the nuisance parameter, right? How does RooFit knows what to float? Does it float all parameters labeled as nuisance?

[StephanH]

Yes, in my opinion a parameter of the fit likelihood should also be a parameter of the constraint likelihood. The “observed” value (the official lumi measurement of your experiment) then has the role of the global observable. In the Gaussian distribution there is no difference if you interchange the two, but you are right that it matters for other constraint functions.
Just read the constraint likelihood as “How likely is it that the luminosity is at the value lumi_par if we measured lumi_obs before”:

Gauss(lumi_obs | lumi_par, sigma)

RooFit floats all variables that are not constant, and not in the dataset (because these are observables).