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");
sample.defineType("det1");
sample.defineType("det2");
sample.defineType("det3");
//create Workspace (necessary to set up ModelConfig)
RooWorkspace w("w");
w.import(RooArgSet(model, sample));
w.Print();
//build workspace for simultaneous model
RooSimWSTool simw(w);
RooSimultaneous* tot_model = simw.build("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 //
tot_model->Print("t");
TFile *f = new TFile("generated_data.root","READ") ;
//f->GetListOfKeys()->Print();
RooDataSet* data_det1 = (RooDataSet*)f->Get("model_det1Data");
RooDataSet* data_det2 = (RooDataSet*)f->Get("model_det2Data");
RooDataSet* data_det3 = (RooDataSet*)f->Get("model_det3Data");
f->Close();
//bin data sets
x.setBins(100);
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->Print("t");
w.Print();
//fit
tot_model->fitTo(combData,Constrain(RooArgSet(eff,offset,sigma))); //
tot_model->Print("t");
//FREQUENTIST ANALYSIS
//create ModelConfig
RooStats::ModelConfig mc("mc","ModelConfig",&w);
mc.SetPdf(*w.pdf("tot_model"));
w.defineSet("POI","inv_half_life");
w.defineSet("nuisP","a0,a1,nbkg,offset,sigma,eff");
mc.SetParametersOfInterest(*w.set("POI"));
mc.SetNuisanceParameters(*w.set("nuisP"));
mc.SetObservables("x");
w.import(mc);
w.Print();
//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;
//plot
LikelihoodIntervalPlot * plot = new LikelihoodIntervalPlot(interval);
plot->SetNPoints(500);
plot->SetMaximum(1.5);
//plot->SetRange(0.008,0.012);
plot->Draw("");
}
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.)