Dear experts,
I am confused about how to configure constraints in ModelConfig. I am trying to implement a gaussian (yield_Lb2pmunu) constraint in ModelConfig to determine the uppper limit using LEP method. But I am puzzled how to configure the ModelCofig. Here is a summary what I am doing:
TFile f("bpmu_wofirstMLPbin.root");
RooWorkspace* ws = (RooWorkspace*)f.Get("w");
RooRealVar *yield_Lb2pmunu = (RooRealVar*)ws->var("yield_Lb2pmunu");
yield_Lb2pmunu->setVal(yield_Lb2pmunu->getVal()); yield_Lb2pmunu->setConstant();
RooRealVar *gyield_Lb2pmunu = new RooRealVar("gyield_Lb2pmunu","gyield_Lb2pmunu",yield_Lb2pmunu->getMin(),yield_Lb2pmunu->getMax());
RooRealVar *sigma1 = new RooRealVar("sigma1","sigma1",yield_Lb2pmunu->getError());
RooGaussian* gau1 = new RooGaussian("gau1","gau1",*gyield_Lb2pmunu,*yield_Lb2pmunu,*sigma1);
RooAbsPdf* pdf = (RooAbsPdf*)ws->pdf("pdf");
RooProdPdf* pdfconst= new RooProdPdf("pdfconst","pdf with constraints",RooArgSet(*gau1, *pdf));
//construct ModelConfig
ModelConfig sbModel("s+b model",ws);
sbModel.SetWorkspace(*ws);
sbModel.SetObservables(*observables);
sbModel.SetPdf(*pdfconst);
sbModel.SetGlobalObservables(RooArgSet(*gyield_Lb2pmunu));
sbModel.SetNuisanceParameters(RooArgSet(*yield_Lb2pmunu));
sbModel.SetName("S+B Model");
so I fix yield_Lb2pmunu to the value from fit and set it to be the mean value of gaussian function (so mean is a fixed value), and the sigma in the gaussian function is the error of yield_Lb2pmunu from fit. and the variable of gaussian function is a new variable gyield_Lb2pmunu with the same range of yield_Lb2pmunu. And in the ModelConfig, I set the fixed yield_Lb2pmunu to be the SetNuisanceParameters and gyield_Lb2pmunu to be SetGlobalObservables.
Could anyone tell me whether I am correct? I am confused about whether I am right with the mean and variable of the gaussian function and which one I need to fix, variable or mean?
Hi Lingzhu!
For the gaussian, it doesn’t matter which one you use, as mean and variable are interchangable for gaussians. In fact, RooFit/HistFactory is internally inconsistent in how this is handled, but that has no real-world effect on your fit.
I think you should add the floating one to your nuisance parameters, but whether that makes any difference whatsoever depends on if you ever use your modelconfig for any of the RooStats tools, or if you configure them manually anyway.
Thanks! Could you please point out where I am wrong in this code?
I tried to add gaussian constraint on parameter Bdpmu_eff_noPID in this code (the value and uncertainty for it are 0.0387692360337 and 0.000640520336317), so I make this paramater as the second parameter of gaussian function and also set it to be the nuisance parameter. and the first parameter of gaussian function is fixed parameter gBdpmu_eff_noPID, and I set it to be the global variable. Am I right with it?
{
using namespace RooStats;
gROOT->ProcessLine(".L ../B2pmu/RooDoubleCBandGauss.cxx+");
R__LOAD_LIBRARY(../B2pmu/RooDoubleCBandGauss_cxx.so);
TFile *file = TFile::Open("../B2pmu/bpmu_wofirstMLPbin.root");
RooWorkspace *w = (RooWorkspace*)file->Get("w");
//construct function
RooRealVar* m_pmu = (RooRealVar*)w->var("m_pmu");
RooExtendPdf* ext_doublecbandgauss_bd = (RooExtendPdf*)w->pdf("ext_doublecbandgauss_bd_1");
RooExtendPdf* expo = (RooExtendPdf*)w->pdf("ext_expo0_1");
RooAddPdf* pdf = new RooAddPdf("pdf","pdf",RooArgList(*ext_doublecbandgauss_bd, *expo));
//generate samples
RooDataSet* data = pdf->generate(*m_pmu,RooFit::Extended());
RooRealVar *Bdpmu_eff_noPID = (RooRealVar*)w->var("Bdpmu_eff_noPID");
RooRealVar *PIDMLP1 = (RooRealVar*)w->var("PIDMLP1"); PIDMLP1->setConstant(); PIDMLP1->setVal(PIDMLP1->getVal());
RooRealVar *frac_MLP1 = (RooRealVar*)w->var("frac_MLP1"); frac_MLP1->setConstant(); frac_MLP1->setVal(frac_MLP1->getVal());
RooRealVar *tau0_1 = (RooRealVar*)w->var("tau0_1"); //tau0_1->setConstant();
RooRealVar *N_bkg_1 = (RooRealVar*)w->var("N_bkg_1"); //N_bkg_1->setConstant();
RooRealVar * BR_Bd = (RooRealVar*)w->var("BR_Bd");
//construct gaussian construct
RooRealVar *gBdpmu_eff_noPID = new RooRealVar("gBdpmu_eff_noPID","gBdpmu_eff_noPID",0.0387692360337,Bdpmu_eff_noPID->getMin(),Bdpmu_eff_noPID->getMax());
gBdpmu_eff_noPID->setConstant();
RooRealVar *sigma2 = new RooRealVar("sigma2","sigma2",0.000640520336317);
RooGaussian* gau2 = new RooGaussian("gau2","gau2",*gBdpmu_eff_noPID,*Bdpmu_eff_noPID,*sigma2);
RooProdPdf* pdfconst= new RooProdPdf("pdfconst","pdf with constraints",RooArgList(*gau2,*pdf));
//configure ModelConfig
ModelConfig sbModel("s+b model",w);
sbModel.SetWorkspace(*w);
sbModel.SetObservables("m_pmu");
sbModel.SetGlobalObservables(RooArgSet(*gBdpmu_eff_noPID));
sbModel.SetNuisanceParameters(RooArgSet(*tau0_1, *N_bkg_1, *Bdpmu_eff_noPID));
sbModel.SetPdf(*pdfconst);
// sbModel.SetPdf(*pdf);
sbModel.SetName("S+B Model");
sbModel.SetParametersOfInterest(RooArgSet(*BR_Bd));
w->writeToFile("Bd2pmu.root");
RooStats::ProfileLikelihoodCalculator plc(*data, sbModel);
plc.SetTestSize(.05);
RooStats::ConfInterval* lrinterval = plc.GetInterval();
// Let's make a plot
TCanvas *dataCanvas = new TCanvas("dataCanvas");
RooStats::LikelihoodIntervalPlot plotInt((RooStats::LikelihoodInterval *)lrinterval);
plotInt.SetTitle("Profile Likelihood Ratio and Posterior for BR_Bd");
plotInt.Draw();
dataCanvas->SaveAs("Profile_LH_Bd.pdf");
dataCanvas->SaveAs("Profile_LH_Bd.png");
// Get Lower and Upper limits from Profile Calculator
cout << "Profile lower limit = " << ((RooStats::LikelihoodInterval *)lrinterval)->LowerLimit(*BR_Bd) << endl;
cout << "Profile upper limit = " << ((RooStats::LikelihoodInterval *)lrinterval)->UpperLimit(*BR_Bd) << endl;
}
but I got this problem when I run it
Processing add_gaussian_constraints_3.C...
^[[1mRooFit v3.60 -- Developed by Wouter Verkerke and David Kirkby^[[0m
Copyright (C) 2000-2013 NIKHEF, University of California & Stanford University
All rights reserved, please read http://roofit.sourceforge.net/license.txt
[#1] INFO:NumericIntegration -- RooRealIntegral::init(doublecbandgauss_bd_1_Int[m_pmu]) using numeric integrator RooIntegrator1D to calculate Int(m_pmu)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(doublecbandgauss_bd_1_Int[m_pmu]) using numeric integrator RooIntegrator1D to calculate Int(m_pmu)
[#1] INFO:Minization -- p.d.f. provides expected number of events, including extended term in likelihood.
[#1] INFO:NumericIntegration -- RooRealIntegral::init(doublecbandgauss_bd_1_Int[m_pmu]) using numeric integrator RooIntegrator1D to calculate Int(m_pmu)
[#1] INFO:Minization -- createNLL: caching constraint set under name CONSTR_OF_PDF_pdfconst_FOR_OBS_m_pmu with 1 entries
[#1] INFO:Minization -- Including the following contraint terms in minimization: (gau2)
[#1] INFO:Minization -- The following global observables have been defined: (gBdpmu_eff_noPID)
[#0] PROGRESS:Minization -- ProfileLikelihoodCalcultor::DoGLobalFit - find MLE
[#0] PROGRESS:Minization -- ProfileLikelihoodCalcultor::DoMinimizeNLL - using Minuit / Migrad with strategy 1
[#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_pdfconst_pdfData_with_constr) Summation contains a RooNLLVar, using its error level
[#1] INFO:Minization -- RooMinimizer::optimizeConst: activating const optimization
[#1] INFO:NumericIntegration -- RooRealIntegral::init(doublecbandgauss_bd_1_Int[m_pmu]) using numeric integrator RooIntegrator1D to calculate Int(m_pmu)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(doublecbandgauss_bd_1_Int[m_pmu]) using numeric integrator RooIntegrator1D to calculate Int(m_pmu)
[#1] INFO:Minization -- The following expressions have been identified as constant and will be precalculated and cached: (doublecbandgauss_bd_1)
[#1] INFO:Minization -- The following expressions will be evaluated in cache-and-track mode: (ext_doublecbandgauss_bd_1,ext_expo0_1)
[#1] INFO:Minization --
RooFitResult: minimized FCN value: -762.283, estimated distance to minimum: 3.55842e-06
covariance matrix quality: Full, accurate covariance matrix
Status : MINIMIZE=0
Floating Parameter FinalValue +/- Error
-------------------- --------------------------
BR_Bd 1.4852e-08 +/- 9.38e-09
Bdpmu_eff_noPID 3.8818e-02 +/- 6.80e-04
N_bkg_1 2.8494e+03 +/- 5.57e+01
tau0_1 -2.3890e-03 +/- 7.45e-05
[#1] INFO:Minization -- RooProfileLL::evaluate(nll_pdfconst_pdfData_with_constr_Profile[BR_Bd]) Creating instance of MINUIT
[#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_pdfconst_pdfData_with_constr) Summation contains a RooNLLVar, using its error level
[#1] INFO:Minization -- RooProfileLL::evaluate(nll_pdfconst_pdfData_with_constr_Profile[BR_Bd]) determining minimum likelihood for current configurations w.r.t all observable
[#1] INFO:Minization -- RooProfileLL::evaluate(nll_pdfconst_pdfData_with_constr_Profile[BR_Bd]) minimum found at (BR_Bd=1.48658e-08)
.Warning: lower value for BR_Bd is at limit 0
[#1] INFO:NumericIntegration -- RooRealIntegral::init(doublecbandgauss_bd_1_Int[m_pmu]) using numeric integrator RooIntegrator1D to calculate Int(m_pmu)
[#1] INFO:Minization -- RooProfileLL::evaluate(nll_pdfconst_pdfData_with_constr_Profile[BR_Bd]) Creating instance of MINUIT
[#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_pdfconst_pdfData_with_constr) Summation contains a RooNLLVar, using its error level
[#1] INFO:Minization -- RooProfileLL::evaluate(nll_pdfconst_pdfData_with_constr_Profile[BR_Bd]) determining minimum likelihood for current configurations w.r.t all observable
[#0] WARNING:Eval -- RooAbsPdf::getLogVal(pdfconst) WARNING: large likelihood value: 7.70965e+68
[#0] WARNING:Eval -- RooAbsPdf::getLogVal(pdfconst) WARNING: large likelihood value: 1.61098e+67
[#0] WARNING:Eval -- RooAbsPdf::getLogVal(pdfconst) WARNING: large likelihood value: 1.78966e+72
[#0] WARNING:Eval -- RooAbsPdf::getLogVal(pdfconst) WARNING: large likelihood value: 9.86313e+64
[#0] WARNING:Eval -- RooAbsPdf::getLogVal(pdfconst) WARNING: large likelihood value: 7.1038e+65
[#0] WARNING:Eval -- RooAbsPdf::getLogVal(pdfconst) WARNING: large likelihood value: 3.33003e+62
Unfortunately, I cannot run your code because you are using some custom libraries.
At first sight, I cannot see anything inherently wrong about this code, but also the warnings you are seeing are not too worrying - it could very well be that your fit converges just fine despite the warnings. Nevertheless, they can probably be resolved by passing ´RooFit::Offset(true)´ to the fit.
Is there any concrete issue you are seeing with the results of your fit?
I can not get any results, it’s just the warnings and it never ends
about the ´RooFit::Offset(true)´, where should I add it? since in my code, these is no fit.
Hi, I found I could use LEP method in $ROOTSYS/tutorials/roostats/StandardHypoTestInvDemo.C+ to get the upper limit, by passing the ModelConfig I get from the above script. The result seems fine. I do not know why I have problem with RooStats::ProfileLikelihoodCalculator
Thanks! Maybe there are some bugs… but I think I get the reasonable result using LEP method.
The thing I would like to make sure is if I add gaussian constraint on parameters correctly. Please correct me if I am wrong (for example, we would like to constaint parameter a to 1, the uncertainty is 0.2):
1, construct gaussian function, the second parameter is the parameter we want to constraint , so it is a, (it needs floating and needs to be nuisance parameter). the first parameter is a fixed value, to which we want to constraint a ,so it is 1. and it needs to be the global observable. (actually, the first and second parameter can swap position and it does not affact the result). the third parameter is the uncertainty, so it it 0.2.
2, multiply the pdf without uncertainty with gaussian function, the final pdf is as the pdf in ModelConfig
