Dear experts,
I obtained different expected upper limits depending on which parameters were treated as nuisance parameters.
In the following test code, the model consists of a Gaussian plus a first-order polynomial with six parameters, Ns, Nb, sigma, mean, c1, c2. A toy dataset is produced from this model and provided to the ModelConfig. When only Nb is treated as a nuisance parameter, the expected upper limit is 85. In contrast, when Nb, sigma, mean, c1, c2 are all treated as nuisance parameters, the expected upper limit improves to 80.
Why are these two results different? Why does setting the nuisance parameters correctly lead to a better expected upper limit?
from ROOT import RooWorkspace, RooArgSet, RooStats
w = RooWorkspace()
# gaus + poly1 as model
w.factory("SUM::model(Ns[50,0,1000]*Gaussian::sig(m[500,400,600], mean[480, 200, 600], sigma[3, 0, 5]), Nb[20000,200,1000000]*Bernstein::bkg(m, {c1[1],c2[1,0,10]}))")
m = w["m"]
m.setBins(100)
model = w["model"]
# generate a test dataset
data1 = model.generate(w["m"], Extended=True)
# configure (s+b) model
sbModel = RooStats.ModelConfig("sbModel","sbModel", w)
sbModel.SetPdf('model')
sbModel.SetObservables('m')
sbModel.SetNuisanceParameters('Nb') # only Nb as nuisance
#sbModel.SetNuisanceParameters('mean,sigma,Nb,c1,c2') # all the parameters as nuisance except Ns
sbModel.SetParametersOfInterest('Ns') # Ns as POI
poi = sbModel.GetParametersOfInterest().first()
poi.setVal(50)
sbModel.SetSnapshot(poi)
# Clone b-only model with POI=0
bModel = sbModel.Clone('bModel')
poi_old = poi.getVal()
poi.setVal(0)
poi.setConstant()
bModel.SetSnapshot({poi})
poi.setVal(poi_old)
poi.setConstant(False)
# Configure HypoTestInverter according to tutorial StandardHypoTestInvDemo.C
# using FrequentistCalculator and within LEP case
Log_DLL = RooStats.SimpleLikelihoodRatioTestStat(sbModel.GetPdf(), bModel.GetPdf())
nllPars = RooArgSet(sbModel.GetNuisanceParameters())
nllPars.add(sbModel.GetSnapshot())
altPars = RooArgSet(bModel.GetNuisanceParameters())
altPars.add(bModel.GetSnapshot())
Log_DLL.SetNullParameters(nllPars)
Log_DLL.SetAltParameters(altPars)
Log_DLL.SetReuseNLL(True)
FreqCalc = RooStats.FrequentistCalculator(data1, bModel, sbModel)
toyMCs = FreqCalc.GetTestStatSampler()
toyMCs.SetTestStatistic(Log_DLL)
toyMCs.SetGenerateBinned(False)
toyMCs.SetUseMultiGen(True)
FreqCalc.SetToys(1000, 1000)
runHypoTest = RooStats.HypoTestInverter(FreqCalc)
#runHypoTest.SetConfidenceLevel(0.90)
runHypoTest.SetConfidenceLevel(0.95)
runHypoTest.UseCLs(True)
runHypoTest.SetVerbose(True)
runHypoTest.SetFixedScan(20, 0, 1000)
r = runHypoTest.GetInterval()
I have one more question regarding external constraints in ModelConfig. I don’t find any examples or tutorials, or discussions of SetExternalConstraints. In my code, Gaussian constraints are added to ModelConfig using sbModel.SetExternalConstraints("myGaussianConstraints") — is this correct? If so, does that mean I no longer need to set the global observables?
Thank you, and I look forward to your reply.