RooStats: Change of Constraint does not change Profile Likelihood Plot

Hello everyone,
I have been trying to set up a Frequentist Analysis of a signal+background (Gaussian+Polynomial) PDF where some parameters should be constrained in the likelihood due to Systematic Uncertainties. However, these parameters should not be fitting parameters.
In the code below you find the initialization of the variables and PDFs, followed by a fitting routine and after that the PLL analysis. However, I found that increasing the uncertainty of one of my parameters (sigma_sigma) did not change the PLL plot the way it should have. I would expect that if I increase sigma_sigma from its original value of 0.2 to 100, the PLL plot should get a lot wider.

I would be happy about any suggestions what I might have done wrong. Thank you!

generated_data.root (220.7 KB)

using namespace std;
using namespace RooFit;
using namespace RooStats;

void simplified_testmodel(){

//Create model
//Energy, Energy offset and offset uncertainty/constraint
RooRealVar x("x","Energy",400,500);

RooRealVar offset("offset","energy offset",0);
RooRealVar offset_mean("offset_mean","mean of offset",0);
RooRealVar offset_sigma("offset_sigma","uncertainty of offset",0.3);
RooGaussian err_offset("err_offset","constraint on offset",offset,offset_mean,offset_sigma);

RooConstVar mean("mean","mean",450);
RooFormulaVar mean_offset("mean_offset","mean including energy offset","mean+offset",RooArgList(mean,offset));

//Energy resolution and resolution uncertainty/constraint
RooRealVar sigma("sigma","sigma",4.0);
RooRealVar sigma_mean("sigma_mean","mean of energy resolution",4.0);
RooRealVar sigma_sigma("sigma_sigma","uncertainty of sigma",0.2);
RooGaussian err_sigma("err_sigma","constraint on sigma",sigma,sigma_mean,sigma_sigma);

//Signal PDF
RooGaussian signal("signal","signal",x,mean_offset,sigma);

//Background Variables and PDF
RooRealVar a0("a0","a0",50,0,7000);
RooRealVar a1("a1","a1",1,-100,5);
RooPolynomial bkg("bkg","Background",x,RooArgList(a0,a1));

//Additional Variables for complete model
RooRealVar inv_half_life("inv_half_life","inverse half life in unit of 10^-22 yr",0.005,0,0.15);
RooRealVar live_time("live_time","GERDA Phase II live time in yr",3);
RooRealVar conv_factor("conv_factor","conversion factor",8608427.6);

//efficiency with uncertainties
RooRealVar eff("eff","eff",2.04e-4);
RooRealVar eff_mean("eff_mean","eff_mean",2.04e-4); //setting a random value here to enable importing into workspace
RooRealVar eff_sigma("eff_sigma","uncertainty of efficiency",2.0e-5);
RooGaussian err_eff("err_eff","constraint of efficiency",eff,eff_mean,eff_sigma);

//number of signal and background events
RooFormulaVar nsig("nsig","number of signal events","eff*conv_factor*live_time*inv_half_life",RooArgList(eff,conv_factor,live_time,inv_half_life));
RooRealVar nbkg("nbkg","number of background events",500,0,10000);

//combining PDFs
RooAddPdf add_model("add_model","model",RooArgList(bkg,signal),RooArgList(nbkg,nsig));
RooProdPdf model("model","model with constraints",RooArgList(add_model,err_eff,err_offset,err_sigma));

//Define catecogy to distinguish det1 and det2 samples  
RooCategory sample("sample","sample");

//create Workspace (necessary to set up ModelConfig)
RooWorkspace w("w");
w.import(RooArgSet(model, sample));

//build workspace for simultaneous model
RooSimWSTool simw(w);
RooSimultaneous* tot_model ="tot_model", "model", SplitParam("sigma,a0,a1,eff,offset,nbkg,live_time,sigma_mean,sigma_sigma,eff_mean,eff_sigma,offset_mean,offset_sigma","sample"));
//nsig is not in the SplitParam() function because RooFormulaVars are automatically split if one of their parameters is split //

TFile *f = new TFile("generated_data.root","READ") ;

RooDataSet* data_det1 = (RooDataSet*)f->Get("model_det1Data"); 
RooDataSet* data_det2 = (RooDataSet*)f->Get("model_det2Data"); 
RooDataSet* data_det3 = (RooDataSet*)f->Get("model_det3Data"); 

//bin data sets
RooDataHist* h_data_det1 = new RooDataHist("h_data_det1", "binned version of data_det1", RooArgSet(x),*data_det1);
RooDataHist* h_data_det2 = new RooDataHist("h_data_det2", "binned version of data_det2", RooArgSet(x),*data_det2);
RooDataHist* h_data_det3 = new RooDataHist("h_inv_data_det2", "binned version of data_det3", RooArgSet(x),*data_det3);  
//define combined data set
RooDataHist combData("combData","combined data sample", x, Index(sample), Import("det1",*h_data_det1), Import("det2",*h_data_det2), Import("det3",*h_data_det3));

tot_model->fitTo(combData,Constrain(RooArgSet(eff,offset,sigma))); //


//create ModelConfig
RooStats::ModelConfig mc("mc","ModelConfig",&w);


//Set up Profile Likelihood Calculator
ProfileLikelihoodCalculator plc(combData,mc);
plc.SetConfidenceLevel(0.9); // 90% interval
LikelihoodInterval* interval = plc.GetInterval();

RooRealVar* firstPOI = (RooRealVar*) mc.GetParametersOfInterest()->first();
double lowerLimit = interval->LowerLimit(*firstPOI);
double upperLimit = interval->UpperLimit(*firstPOI);

cout << "\n90% interval on " <<firstPOI->GetName()<<" is : ["<<     lowerLimit << ", "<<     upperLimit <<"] "<<endl; 

LikelihoodIntervalPlot * plot = new LikelihoodIntervalPlot(interval);


Please read tips for efficient and successful posting and posting code

ROOT Version: 6.06/08
Platform: Not Provided
Compiler: Not Provided

(I’m unsure about Platform and Compiler - If someone can tell me how to I can find them, I will add them.)

Hi @bibsession ,
and welcome to the ROOT forum! I think we need our RooStats expert @moneta here, let’s ping him.

1 Like


I see in your model that the parameter sigma is fixed. In this case the constraint is just a constant parameter without any effect. This is true also for the other Gaussian constraints.
You should make the sigma floating, for example:

RooRealVar sigma("sigma","sigma",4.0, 0., 10);



1 Like

Thank you very much, that did the trick!